在 python scipy 中实现 Kolmogorov Smirnov 测试

Posted

技术标签:

【中文标题】在 python scipy 中实现 Kolmogorov Smirnov 测试【英文标题】:Implementing a Kolmogorov Smirnov test in python scipy 【发布时间】:2011-12-15 18:44:51 【问题描述】:

我有一个关于 N 个数字的数据集,我想测试它的正态性。 我知道 scipy.stats 有一个 kstest function 但是没有关于如何使用它以及如何解释结果的示例。 有哪位熟悉的可以给我一些建议吗?

根据文档,使用 kstest 返回两个数字,KS 检验统计量 D 和 p 值。 如果 p 值大于显着性水平(例如 5%),那么我们不能拒绝数据来自给定分布的假设。

当我通过从正态分布中抽取 10000 个样本并测试高斯性来进行测试运行时:

import numpy as np
from scipy.stats import kstest

mu,sigma = 0.07, 0.89
kstest(np.random.normal(mu,sigma,10000),'norm')

我得到以下输出:

(0.04957880905196102, 8.9249710700788814e-22)

p 值小于 5%,这意味着我们可以拒绝数据呈正态分布的假设。但是样本是从正态分布中抽取的!

有人可以理解并向我解释这里的差异吗?

(正态性测试是否假设 mu = 0 和 sigma = 1?如果是,我如何测试我的数据是高斯分布但具有不同的 mu 和 sigma?)

【问题讨论】:

【参考方案1】:

您的数据是使用 mu=0.07 和 sigma=0.89 生成的。 您正在根据平均值为 0 且标准差为 1 的正态分布测试此数据。

原假设 (H0) 是您的数据作为样本的分布等于均值为 0,标准差为 1 的标准正态分布。

较小的 p 值表明与 D 一样大的检验统计量预期具有概率 p 值。

换句话说,(p 值约为 8.9e-22)H0 是真的不太可能。

这是合理的,因为均值和标准偏差不匹配。

将您的结果与:

In [22]: import numpy as np
In [23]: import scipy.stats as stats
In [24]: stats.kstest(np.random.normal(0,1,10000),'norm')
Out[24]: (0.007038739782416259, 0.70477679457831155)

要测试您的数据是否为高斯分布,您可以对其进行移位和重新缩放,使其正常,均值为 0,标准偏差为 1:

data=np.random.normal(mu,sigma,10000)
normed_data=(data-mu)/sigma
print(stats.kstest(normed_data,'norm'))
# (0.0085805670733036798, 0.45316245879609179)

警告: (many thanks to user333700 (aka scipy developer Josef Perktold)) 如果你不知道musigma,估计参数会使p 值无效:

import numpy as np
import scipy.stats as stats

mu = 0.3
sigma = 5

num_tests = 10**5
num_rejects = 0
alpha = 0.05
for i in xrange(num_tests):
    data = np.random.normal(mu, sigma, 10000)
    # normed_data = (data - mu) / sigma    # this is okay
    # 4915/100000 = 0.05 rejects at rejection level 0.05 (as expected)
    normed_data = (data - data.mean()) / data.std()    # this is NOT okay
    # 20/100000 = 0.00 rejects at rejection level 0.05 (not expected)
    D, pval = stats.kstest(normed_data, 'norm')
    if pval < alpha:
        num_rejects += 1
ratio = float(num_rejects) / num_tests
print('/ = :.2f rejects at rejection level '.format(
    num_rejects, num_tests, ratio, alpha))     

打印

20/100000 = 0.00 rejects at rejection level 0.05 (not expected)

这表明stats.kstest 可能不会拒绝预期数量的空假设 如果使用样本的均值和标准差对样本进行归一化

normed_data = (data - data.mean()) / data.std()    # this is NOT okay

【讨论】:

嗨,unutbu,这是有道理的,谢谢!如何测试我的数据是否以示例中给出的不同 mu 和 sigma 值高斯分布?如果我可以将高斯定义为 g = lambda x: ( 1.0/(np.sqrt(2.0*np.pisigma2)) ) exp(-((x -mu)**2/(2.0*sigma2) ) ) ,然后我可以将 kstest 作为 kstest(data,g) 运行吗? 您可以提供一个可调用对象作为kstest 的第二个参数,但可调用对象必须表示累积分布函数,而不是概率密度函数。我认为在这种情况下,重新规范化您的数据会更容易(如上所示)。 Kolmogorov-Smirnov 检验假定未估计参数。如果您标准化或使用估计参数,则 kstest 的 pvalues 不正确。如果您只是想测试数据是否对于 mu 和 sigma 是正态分布的,那么我建议您进行其他测试。 嗨,有哪些测试只是为了测试? @user528025:嗨,我不明白你的意思。如果您要问统计问题,最好在Cross Validated 发帖提问。【参考方案2】:

关于 unutbu 答案的更新:

对于仅取决于位置和尺度但没有形状参数的分布,几个拟合优度检验统计量的分布与位置和尺度值无关。分布是非标准的,但是,它可以制成表格并与基础分布的任何位置和规模一起使用。

具有估计位置和规模的正态分布的 Kolmogorov-Smirnov 检验也称为Lilliefors test。

它现在可以在 statsmodels 中使用,具有相关决策范围的近似 p 值。

>>> import numpy as np
>>> mu,sigma = 0.07, 0.89
>>> x = np.random.normal(mu, sigma, 10000)
>>> import statsmodels.api as sm
>>> sm.stats.lilliefors(x)
(0.0055267411213540951, 0.66190841161592895)

大多数蒙特卡罗研究表明,Anderson-Darling 检验比 Kolmogorov-Smirnov 检验更有效。它在具有临界值的 scipy.stats 和具有近似 p 值的 statsmodels 中可用:

>>> sm.stats.normal_ad(x)
(0.23016468240712129, 0.80657628536145665)

这两项检验都不拒绝样本是正态分布的 Null 假设。 而问题中的 kstest 拒绝了样本是标准正态分布的 Null 假设。

【讨论】:

Anderson(或 Shapiro 或 D'Agostino)都因样本量大而失败 他们失败了是什么意思?您可以单独询问问题的详细信息。【参考方案3】:

您可能还想考虑使用 Shapiro-Wilk 检验,该检验“检验数据来自正态分布的零假设”。它也在scipy中实现:

http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.shapiro.html

您需要将数据直接传递到函数中。

import scipy

W, p = scipy.stats.shapiro(dataset)
print("Shapiro-Wilk test statistic, W:", W, "\n", "p-value:", p)

返回类似:

 Shapiro-Wilk test statistic, W: 0.7761164903640747 
 p-value: 6.317247641091492e-37

对于 p

【讨论】:

链接文档指出,对于 N > 5000,p 值不能保证是正确的。原题以10000为例。【参考方案4】:

作为对@unutbu 答案的补充,您还可以在 kstest 中提供测试分布的分布参数。假设我们有一些来自变量的样本(并将它们命名为 datax),我们想检查这些样本是否可能来自对数正态、均匀或正态。请注意,对于 scipy stats,为每个分布采用输入参数的方式会有所不同。现在,感谢 kstest 中的“args”(元组或序列),可以为您要测试的 scipy.stats 分布提供参数。

:) 我还添加了使用双样本测试的选项,以防您想以任何一种方式进行:

import numpy as np
from math import sqrt
from scipy.stats import kstest, ks_2samp, lognorm
import scipy.stats

def KSSeveralDists(data,dists_and_args,samplesFromDists=100,twosampleKS=True):
    returnable=
    for dist in dists_and_args:
        try:
            if twosampleKS:
                try:
                    loc=dists_and_args[dist][0]
                    scale=dists_and_args[dist][1]
                    expression='scipy.stats.'+dist+'.rvs(loc=loc,scale=scale,size=samplesFromDists)'
                    sampledDist=eval(expression)
                except:
                    sc=dists_and_args[dist][0]
                    loc=dists_and_args[dist][1]
                    scale=dists_and_args[dist][2]
                    expression='scipy.stats.'+dist+'.rvs(sc,loc=loc,scale=scale,size=samplesFromDists)'
                    sampledDist=eval(expression)
                D,p=ks_2samp(data,sampledDist)
            else:
                D,p=kstest(data,dist,N=samplesFromDists,args=dists_and_args[dist])
        except:
            continue
        returnable[dist]='KS':D,'p-value':p
    return returnable

a=lambda m,std: m-std*sqrt(12.)/2.
b=lambda m,std: m+std*sqrt(12.)/2.
sz=2000

sc=0.5 #shape 
datax=lognorm.rvs(sc,loc=0.,scale=1.,size=sz)
normalargs=(datax.mean(),datax.std())

#suppose these are the parameters you wanted to pass for each distribution
dists_and_args='norm':normalargs,
               'uniform':(a(*normalargs),b(*normalargs)),
               'lognorm':[0.5,0.,1.]
              
print "two sample KS:"
print KSSeveralDists(datax,dists_and_args,samplesFromDists=sz,twosampleKS=True)
print "one sample KS:"
print KSSeveralDists(datax,dists_and_args,samplesFromDists=sz,twosampleKS=False)

它给出的输出类似于:

两个样本KS: 'lognorm': 'KS': 0.023499999999999965, 'p-value': 0.63384188886455217, 'norm': 'KS': 0.10600000000000004, 'p-value': 2.918766666723155e-10, KS':0.15300000000000002,'p值':6.443660021191129e-21

一个样本 KS: 'lognorm':'KS':0.01763415915126032,'p-value':0.56275820961065193,'norm':'KS':0.10792612430093562,'p-value':0.0,'uniform':'KS': 0.14910036159697559,“p值”:0.0

注意:对于 scipy.stats 均匀分布,a 和 b 被视为 a=loc 和 b=loc + scale(参见documentation)。

【讨论】:

以上是关于在 python scipy 中实现 Kolmogorov Smirnov 测试的主要内容,如果未能解决你的问题,请参考以下文章

如何在 Scipy 中实现日志均匀分布?

如何在 Python 中实现 KS-Test

python / scipy中的多元样条插值?

结合scipy.linalg在Python中使用线性系统

Python中的约束最小二乘估计

python中的离散优化