Внедрение теста Колмогорова Смирнова в python scipy

Ответ 1

Ваши данные были сгенерированы с помощью mu = 0,07 и сигма = 0,89. Вы тестируете эти данные против нормального распределения со средним значением 0 и стандартным отклонением 1.

Нулевая гипотеза (H0) заключается в том, что распределение, в котором ваши данные являются образцом, равно стандартным нормальным распределениям со средним 0, std отклонение 1.

Небольшое значение p указывает на то, что ожидаемая статистическая статистика, равная D, с вероятностью p-значения.

Иными словами, (с p-значением ~ 8.9e-22) маловероятно, что H0 истинно.

Это разумно, поскольку средства и отклонения std не совпадают.

Сравните ваш результат с:

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 и std отклонение 1:

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

Предупреждение: (большое спасибо user333700 (aka scipy developer Josef Perktold)) Если вы не знаете mu и sigma, оценка параметров делает неверным значение 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

Ответ 2

Обновление для ответа unutbu:

Для распределений, которые зависят только от местоположения и масштаба, но не имеют параметра формы, распределения нескольких статистических тестов пригодности не зависят от значений местоположения и шкалы. Распределение нестандартное, однако оно может быть сведено в таблицу и использоваться с любым местоположением и масштабом лежащего в основе распределения.

Тест Колмогорова-Смирнова для нормального распределения с оценкой местоположения и масштаба также называется тест Lilliefors.

Теперь он доступен в 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)

Большинство исследований в Монте-Карло показывают, что тест Андерсона-Дарлин более мощный, чем тест Колмогорова-Смирнова. Он доступен в scipy.stats с критическими значениями и в statsmodels с приблизительными значениями p:

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

Ни один из тестов не отклоняет гипотезу Null о том, что образец нормально распределен. В то время как kstest в вопросе отклоняет гипотезу Null, что образец стандартный нормальный распределен.

Ответ 3

Вы также можете рассмотреть возможность использования теста Шапиро-Вилка, который "проверяет нулевую гипотезу о том, что данные были получены из нормального распределения". Он также реализован в 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 < 0,01 (или 0,05, если вы предпочитаете - это не имеет значения), у нас есть все основания отклонить нулевую гипотезу о том, что эти данные были получены из нормального распределения.

Ответ 4

В дополнение к ответу @unutbu вы также можете предоставить параметры распределения для тестового дистрибутива в kstest. Предположим, что у нас были некоторые выборки из переменной (и назвали их datax), и мы хотели проверить, могут ли эти образцы не прийти из логнормального, унифицированного или нормального. Обратите внимание, что для scipy stats способ ввода входных параметров для каждого дистрибутива немного меняется. Теперь, благодаря "args" (кортеж или последовательность) в kstest, возможно предоставить аргументы для дистрибутива 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}, 'uniform': {' KS ': 0.15300000000000002,' p-value ': 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-value': 0.0}}

Примечание. Для равномерного распределения scipy.stats a и b берутся как a = loc и b = loc + scale (см. документация).