Анализ одной выборки#
Сначала мы создаем несколько случайных величин. Мы устанавливаем seed, чтобы при каждом запуске получать идентичные результаты для анализа. В качестве примера мы берем выборку из распределения Стьюдента:
>>> import numpy as np
>>> import scipy.stats as stats
>>> x = stats.t.rvs(10, size=1000)
Здесь мы устанавливаем требуемый параметр формы t-распределения, который в статистике соответствует степеням свободы, равным 10. Использование size=1000 означает, что наша выборка состоит из 1000 независимо сгенерированных (псевдо)случайных чисел. Поскольку мы не указали ключевые аргументы loc и scale#14990
Описательная статистика#
x является массивом numpy, и у нас есть прямой доступ ко всем методам массива, например,
>>> print(x.min()) # equivalent to np.min(x)
-3.78975572422 # random
>>> print(x.max()) # equivalent to np.max(x)
5.26327732981 # random
>>> print(x.mean()) # equivalent to np.mean(x)
0.0140610663985 # random
>>> print(x.var()) # equivalent to np.var(x))
1.28899386208 # random
Как свойства выборки соотносятся с их теоретическими аналогами?
>>> m, v, s, k = stats.t.stats(10, moments='mvsk')
>>> n, (smin, smax), sm, sv, ss, sk = stats.describe(x)
>>> sstr = '%-14s mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
>>> print(sstr % ('distribution:', m, v, s ,k))
distribution: mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000 # random
>>> print(sstr % ('sample:', sm, sv, ss, sk))
sample: mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556 # random
Примечание: stats.describe использует несмещенную оценку для дисперсии, в то время как np.var является смещенной оценкой.
Для нашей выборки выборочные статистики немного отличаются от их теоретических аналогов.
T-тест и KS-тест#
Мы можем использовать t-тест, чтобы проверить, отличается ли среднее значение нашей выборки статистически значимым образом от теоретического ожидания.
>>> print('t-statistic = %6.3f pvalue = %6.4f' % stats.ttest_1samp(x, m))
t-statistic = 0.391 pvalue = 0.6955 # random
P-значение равно 0,7, это означает, что при ошибке альфа, например, 10%, мы не можем отвергнуть гипотезу о том, что выборочное среднее равно нулю, ожиданию стандартного t-распределения.
В качестве упражнения мы можем рассчитать наш t-тест также напрямую, без использования предоставленной функции, что должно дать нам тот же ответ, и так оно и есть:
>>> tt = (sm-m)/np.sqrt(sv/float(n)) # t-statistic for mean
>>> pval = stats.t.sf(np.abs(tt), n-1)*2 # two-sided pvalue = Prob(abs(t)>tt)
>>> print('t-statistic = %6.3f pvalue = %6.4f' % (tt, pval))
t-statistic = 0.391 pvalue = 0.6955 # random
Критерий Колмогорова-Смирнова может использоваться для проверки гипотезы о том, что выборка происходит из стандартного t-распределения
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,)))
KS-statistic D = 0.016 pvalue = 0.9571 # random
Опять же, p-значение достаточно велико, чтобы мы не могли отвергнуть гипотезу о том, что случайная выборка действительно распределена согласно t-распределению. В реальных приложениях мы не знаем, каково базовое распределение. Если мы выполним тест Колмогорова-Смирнова для нашей выборки против стандартного нормального распределения, то мы также не сможем отвергнуть гипотезу о том, что наша выборка была сгенерирована нормальным распределением, учитывая, что в этом примере p-значение составляет почти 40%.
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 'norm'))
KS-statistic D = 0.028 pvalue = 0.3918 # random
Однако стандартное нормальное распределение имеет дисперсию 1, в то время как наша выборка имеет дисперсию 1.29. Если мы стандартизируем нашу выборку и проверим её против нормального распределения, то p-значение снова достаточно велико, чтобы мы не могли отвергнуть гипотезу о том, что выборка получена из нормального распределения.
>>> d, pval = stats.kstest((x-x.mean())/x.std(), 'norm')
>>> print('KS-statistic D = %6.3f pvalue = %6.4f' % (d, pval))
KS-statistic D = 0.032 pvalue = 0.2397 # random
Примечание: Критерий Колмогорова-Смирнова предполагает, что мы тестируем против распределения с заданными параметрами, поскольку в последнем случае мы оценили среднее и дисперсию, это предположение нарушено, и распределение тестовой статистики, на которой основано p-значение, не является корректным.
Хвосты распределения#
Наконец, мы можем проверить верхний хвост распределения. Мы можем использовать функцию процентной точки ppf, которая является обратной к функции cdf, чтобы получить критические значения, или, более прямо, мы можем использовать обратную функцию выживания
>>> crit01, crit05, crit10 = stats.t.ppf([1-0.01, 1-0.05, 1-0.10], 10)
>>> print('critical values from ppf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % (crit01, crit05, crit10))
critical values from ppf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> print('critical values from isf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' % tuple(stats.t.isf([0.01,0.05,0.10],10)))
critical values from isf at 1%, 5% and 10% 2.7638 1.8125 1.3722
>>> freq01 = np.sum(x>crit01) / float(n) * 100
>>> freq05 = np.sum(x>crit05) / float(n) * 100
>>> freq10 = np.sum(x>crit10) / float(n) * 100
>>> print('sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f' % (freq01, freq05, freq10))
sample %-frequency at 1%, 5% and 10% tail 1.4000 5.8000 10.5000 # random
Во всех трех случаях наша выборка имеет больший вес в верхнем хвосте, чем исходное распределение. Мы можем кратко проверить большую выборку, чтобы увидеть, получим ли мы более близкое соответствие. В этом случае эмпирическая частота довольно близка к теоретической вероятности, но если мы повторим это несколько раз, колебания все еще довольно велики.
>>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100
>>> print('larger sample %%-frequency at 5%% tail %8.4f' % freq05l)
larger sample %-frequency at 5% tail 4.8000 # random
Мы также можем сравнить его с хвостом нормального распределения, который имеет меньший вес в хвостах:
>>> print('tail prob. of normal at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f' %
... tuple(stats.norm.sf([crit01, crit05, crit10])*100))
tail prob. of normal at 1%, 5% and 10% 0.2857 3.4957 8.5003
Критерий хи-квадрат может использоваться для проверки того, отличаются ли наблюдаемые частоты для конечного числа интервалов значимо от вероятностей предполагаемого распределения.
>>> quantiles = [0.0, 0.01, 0.05, 0.1, 1-0.10, 1-0.05, 1-0.01, 1.0]
>>> crit = stats.t.ppf(quantiles, 10)
>>> crit
array([ -inf, -2.76376946, -1.81246112, -1.37218364, 1.37218364,
1.81246112, 2.76376946, inf])
>>> n_sample = x.size
>>> freqcount = np.histogram(x, bins=crit)[0]
>>> tprob = np.diff(quantiles)
>>> nprob = np.diff(stats.norm.cdf(crit))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 2.30 pvalue = 0.8901 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 64.60 pvalue = 0.0000 # random
Мы видим, что стандартное нормальное распределение явно отвергается, в то время как стандартное t-распределение не может быть отвергнуто. Поскольку дисперсия нашей выборки отличается от обоих стандартных распределений, мы можем снова повторить тест, учитывая оценку масштаба и местоположения.
Метод fit распределений может использоваться для оценки параметров распределения, и тест повторяется с использованием вероятностей оцененного распределения.
>>> tdof, tloc, tscale = stats.t.fit(x)
>>> nloc, nscale = stats.norm.fit(x)
>>> tprob = np.diff(stats.t.cdf(crit, tdof, loc=tloc, scale=tscale))
>>> nprob = np.diff(stats.norm.cdf(crit, loc=nloc, scale=nscale))
>>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
>>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
>>> print('chisquare for t: chi2 = %6.2f pvalue = %6.4f' % (tch, tpval))
chisquare for t: chi2 = 1.58 pvalue = 0.9542 # random
>>> print('chisquare for normal: chi2 = %6.2f pvalue = %6.4f' % (nch, npval))
chisquare for normal: chi2 = 11.08 pvalue = 0.0858 # random
Учитывая оцененные параметры, мы всё ещё можем отвергнуть гипотезу, что наша выборка взята из нормального распределения (на уровне 5%), но снова, с p-значением 0.95, мы не можем отвергнуть t-распределение.
Специальные тесты для нормальных распределений#
Поскольку нормальное распределение является наиболее распространенным распределением в статистике, существует несколько дополнительных функций для проверки, могла ли выборка быть взята из нормального распределения.
Сначала мы можем проверить, отличаются ли асимметрия и эксцесс нашей выборки значимо от таковых у нормального распределения:
>>> print('normal skewtest teststat = %6.3f pvalue = %6.4f' % stats.skewtest(x))
normal skewtest teststat = 2.785 pvalue = 0.0054 # random
>>> print('normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x))
normal kurtosistest teststat = 4.757 pvalue = 0.0000 # random
Эти два теста объединены в тест на нормальность
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x))
normaltest teststat = 30.379 pvalue = 0.0000 # random
Во всех трёх тестах p-значения очень низкие, и мы можем отвергнуть гипотезу о том, что наша выборка имеет асимметрию и эксцесс нормального распределения.
Поскольку асимметрия и эксцесс нашей выборки основаны на центральных моментах, мы получаем точно такие же результаты, если тестируем стандартизированную выборку:
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest((x-x.mean())/x.std()))
normaltest teststat = 30.379 pvalue = 0.0000 # random
Поскольку нормальность отвергается столь сильно, мы можем проверить, дает ли normaltest разумные результаты для других случаев:
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.t.rvs(10, size=100)))
normaltest teststat = 4.698 pvalue = 0.0955 # random
>>> print('normaltest teststat = %6.3f pvalue = %6.4f' %
... stats.normaltest(stats.norm.rvs(size=1000)))
normaltest teststat = 0.613 pvalue = 0.7361 # random
При проверке нормальности небольшой выборки t-распределенных наблюдений и большой выборки нормально распределенных наблюдений, в обоих случаях мы не можем отвергнуть нулевую гипотезу о том, что выборка происходит из нормального распределения. В первом случае это происходит потому, что тест недостаточно мощный, чтобы различить t-распределенную и нормально распределенную случайную переменную в небольшой выборке.