怎么用python进行grs检验因子模型 python q检验

假设检验步骤

1. 给出原假设$H_0$，通常为积极肯定的一面，例如原数据集符合某类分布$F(X)$。
2. 挑选统计量(该统计量服从分布$F(X)$)，根据样本计算统计量的值。
3. 根据预先设定的显著性程度或者置信度，计算临界值，若统计值超出临界值则否定原假设；或者根据统计值计算p值(符合原假设的概率)，若p值小于显著性程度则否定原假设。

基于上述思想，可以先假设没有异常值且样本服从某种分布，若检验结果否定原假设，则存在异常值。

异常值检验

单变量检验

箱型图

上下边界之外的点就是异常点。
四分位距(IQR)=上四分位(75%)-下四分位(25%)
上边界=上四分位+1.5*IQR
下边界=下四分位-1.5*IQR

3$\sigma$准则

适合于符合正态分布的数据集的检验。
数据位于$(\mu-3\sigma,\mu+3\sigma)$中的概率为0.9973，若某样本超出该范围则有理由相信为异常值。

格拉布斯检验(Grubbs)

3$\sigma$准则的加强版，适合于近似正态分布数据的检验。

统计量为$G=max\frac{|Y_i-\bar Y|}{s}$式中$\bar Y$为样本均值，s为样本标准差。若满足：

式中，t表示t分布，$\alpha$为显著性程度。

#添加包 pip install outlier_utils==0.0.3
from outliers import smirnov_grubbs as grubbs
grubbs.test(data, alpha) # 默认双边检验

ESD检验(generalized Extreme Studentized Deviate)

适合多个异常值检验。预先设定异常值个数上限r，执行r轮单独的检验，每轮检验的变量为$R_i=max\frac{|Y_i-\bar Y|}{s}$相应的临界值为

式中，$t_{p,n-i-1}$指自由度为$n-i-1$的t分布的100p($p=1-\frac{\alpha}{2(n-i+1)}$)百分位点对应的值。

若离群点为满足$R_i>\lambda_i$的点。

import math
import numpy as np
from scipy import stats


def get_r(arr):
    m_arr = np.mean(arr)
    d_arr = abs(arr - m_arr)
    s = np.std(arr, ddof=1)  # 样本标准差，注意分母n-1
    out_ind = np.argmax(d_arr)
    return np.max(d_arr) / s, out_ind


def esd(data, alpha=0.05, max_anoms=0.10):
    n = len(data)
    if isinstance(max_anoms, float):
        r = math.ceil(max_anoms * n)
    else:
        r = max_anoms
    outliers = []
    for i in range(1, r + 1):
        p = 1 - alpha / (n - i + 1) / 2
        t = stats.t.ppf(p, n - i - 1)  # p分位点
        _lambda = (n - i) * t / math.sqrt((n - i - 1 + t ** 2) * (n - i + 1))
        arr = np.delete(data, outliers)
        _r, out_ind = get_r(arr)
        if _r > _lambda:  # 超出临界值，视为异常点
            outliers.append(out_ind)
        else:
            break
    return np.delete(data, outliers), data[outliers]

时间序列中异常值检验

STL算法将时间序列分解为长期趋势项T+季节变动项S+余项R，然后提取余项，并假设其服从正态分布，那么接下来就可以采用ESD检验了。

1. S-ESD检验(Seasonal ESD)
   S-ESD将趋势项T替换成中位数。
2. S-H-ESD检验(Seasonal Hybrid ESD)
   S-H-ESD将趋势项T替换成绝对中位差(MAD,Median Absolute Deviation)。
   建议参考Twitter工程师原论文

from pyculiarity import detect_ts
res = detect_ts(df, max_anoms=0.10, direction='pos',alpha=0.05) # 仅检测正向序列
# rstl包中还含有stl分解时间序列的函数

高维数据异常值检测

局部异常因子(LOF,Local Outlier Factor)算法

参考博客异常检测—LOF算法简介以及Python实现。

from sklearn.neighbors import LocalOutlierFactor
lof = LocalOutlierFactor(n_neighbors=5)
y_pred = lof.fit_predict(data)
score = -lof.negative_outlier_factor_

孤立森林(isolation forest)算法

参考博客异常检测—IsolationForest算法简介以及Python实现

from sklearn.ensemble import IsolationForest
IF = IsolationForest(max_samples=0.9, contamination=0.1) # max_samples指每颗基本树使用的样本量，contamination指异常值比例
res = IF.fit_predict(data)

一类支持向量机(OneClass SVM)

根据sklearn官方文档的介绍，outlier detetction指原数据集中被污染，目的是找出其中的离群点。novelty detection指原数据集未被污染，检验新的观测样本是否为异常点。OneClass由于本身对离群点敏感，因此仅适合于novelty detection。
OneClass是非监督学习算法，算法旨在寻找训练样本的边界，当新的样本点超出边界时，则视为异常点。

from sklearn.svm import OneClassSVM
import numpy as np

X = 0.3 * np.random.randn(100, 2)
X_train = np.r_[X + 2, X - 2]
X_outliers = np.random.uniform(low=-4, high=4, size=(20, 2))
clf = OneClassSVM(nu=0.1, kernel="rbf", gamma=0.1) # nu指错误预测样本(边界外样本)比例上限，支持向量下限，因为训练数据集通常未被污染，因此该参数值通常较小
clf.fit(X_train)
y_pred_train = clf.predict(X_train)
y_pred_outliers = clf.predict(X_outliers)

相关性检验

KS检验与卡方检验都是采用实际频数与期望频数之差来检验某一样本是否符合某一分布。

KS检验

适合连续值的检验。

1. 单样本
   给定理论分布函数$F(X)$(通常理论分布的参数可由样本计算得到)，检验样本是否来自该分布。计算样本的累积经验分布函数$F_i(X)$。则统计量为$D = max|F_I(X)-F(X)|$查表得到临界值，若D大于临界值则表示不符合。
2. 双样本
   与单样本检验原理相同，只不过将理论分布换成另一个样本的累积经验分布函数。

卡方检验

适合于离散值的检验(也可用于连续值，只需对连续值分组即可，但数据利用率就低了，不建议使用)。

离散变量有个理论分布列(例如抛硬币正反概率为1/2)，可求得期望频数E，根据样本计算得实际频数A，那么统计量为$\chi ^2=\sum \frac{(A-E)^2}{E}$该统计量服从卡方分布，查表可得临界值。

卡方检验的另一个重要应用：检验变量是否相关。

通常假设变量X与Y无关，根据原假设计算Y的理论值，再与Y的实际值比较计算得到$\chi ^2$的值，若超出临界值则表示有关。(参考资料)

检验是否相关时，通常原假设H0是不相关，此时便于计算理论值。例如参考资料中假设“是否线上买生鲜”与“性别”不相关时，理论频率计算如下：

计算卡方值为58.66，对应p值是1.13e-12。显著性水平（原假设为真时犯错的概率）$\alpha$=0.05时，对应分位数为3.84。犯错概率如此低的情况下仍旧出现，表明原假设为假，“是否线上买生鲜”与“性别”是相关的。

from scipy import stats
stats.kstest(x, 'norm') # 检验样本是否符合正态分布
stats.ks_2samp(x1, x2) # 检验两样本是否来自同一分布
stats.chisquare(f_obs,f_exp) # 若不传入f_exp则默认是均匀分布的，值为样本均值

from sklearn.feature_selection import SelectKBest, chi2
SelectKBest(chi2,k=10).fit(X,y)

注意：sklearn中用chi2进行卡方值计算时，要求每个特征只能是相对于类标签的布尔值或频数，不能是变量的数值标签（例如LabelEncoder编码），便于抹去自由度的影响。如果特征有多个取值，不能使用sklearn的chi2进行变量筛选，此时统计的卡方值是不正确的，p值也没有意义。

参考资料

http://pdfs.semanticscholar.org/85c3/099bc202a0b175039c8675cc8e01c8f9d256.pdf
https://www.jianshu.com/p/807b2c2bfd9b
https://scikit-learn.org/stable/auto_examples/svm/plot_oneclass.html#sphx-glr-auto-examples-svm-plot-oneclass-py
https://www.jianshu.com/p/64974c4de9d4