样本分析

首先，我们创建一些随机变量。我们设置一个种子所以每次我们都可以得到相同的结果以便观察。
作为一个例子，我们从t分布中抽一个样本。

  1. >>> np.random.seed(282629734)
  2. >>> x = stats.t.rvs(10, size=1000)

这里，我们设置了t分布的形态参数，在这里就是自由度，设为10.使用size=1000表示
我们的样本由1000个抽样是独立的（伪）。当我们不指派loc和scale时，它们具有默认值0和1.

描述统计

X是一个numpy数组。我们可以直接调用它的方法。

  1. >>> print x.max(), x.min()  # equivalent to np.max(x), np.min(x)
  2. 5.26327732981 -3.78975572422
  3. >>> print x.mean(), x.var() # equivalent to np.mean(x), np.var(x)
  4. 0.0140610663985 1.28899386208

如何比较分布本身和它的样本的指标？

  1. >>> m, v, s, k = stats.t.stats(10, moments='mvsk')
  2. >>> n, (smin, smax), sm, sv, ss, sk = stats.describe(x)
  1. >>> print 'distribution:',
  2. distribution:
  3. >>> sstr = 'mean = %6.4f, variance = %6.4f, skew = %6.4f, kurtosis = %6.4f'
  4. >>> print sstr %(m, v, s ,k)
  5. mean = 0.0000, variance = 1.2500, skew = 0.0000, kurtosis = 1.0000
  6. >>> print 'sample:      ',
  7. sample:
  8. >>> print sstr %(sm, sv, ss, sk)
  9. mean = 0.0141, variance = 1.2903, skew = 0.2165, kurtosis = 1.0556

注意：stats.describe用的是无偏的方差估计量，而np.var却用的是有偏的估计量。

T检验和KS检验

我们可以使用t检验是否样本与给定均值（这里是理论均值）存在统计显著差异。

  1. >>> print 't-statistic = %6.3f pvalue = %6.4f' %  stats.ttest_1samp(x, m)
  2. t-statistic =  0.391 pvalue = 0.6955

P值为0.7，这代表第一类错误的概率，在例子中，为10%。我们不能拒绝“该样本均值为0”这个假设，
0是标准t分布的理论均值。

  1. >>> tt = (sm-m)/np.sqrt(sv/float(n))  # t-statistic for mean
  2. >>> pval = stats.t.sf(np.abs(tt), n-1)*2  # two-sided pvalue = Prob(abs(t)>tt)
  3. >>> print 't-statistic = %6.3f pvalue = %6.4f' % (tt, pval)
  4. t-statistic =  0.391 pvalue = 0.6955

这里Kolmogorov-Smirnov检验（KS检验）被被用来检验样本是否来自一个标准的t分布。

  1. >>> print 'KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x, 't', (10,))
  2. KS-statistic D =  0.016 pvalue = 0.9606

又一次得到了很高的P值。所以我们不能拒绝样本是来自t分布的假设。在实际应用中，
我们不能知道潜在的分布到底是什么。如果我们使用KS检验我们的样本对照正态分布，
那么我们将也不能拒绝我们的样本是来自正态分布的，在这种情况下P值为0.4左右。

  1. >>> print 'KS-statistic D = %6.3f pvalue = %6.4f' % stats.kstest(x,'norm')
  2. KS-statistic D =  0.028 pvalue = 0.3949

无论如何，标准正态分布有1的方差，当我们的样本有1.29时。如果我们标准化我们的样本并且
测试它比照正态分布，那么P值将又一次很高我们将还是不能拒绝假设是来自正态分布的。

  1. >>> d, pval = stats.kstest((x-x.mean())/x.std(), 'norm')
  2. >>> print 'KS-statistic D = %6.3f pvalue = %6.4f' % (d, pval)
  3. KS-statistic D =  0.032 pvalue = 0.2402

注释：KS检验假设我们比照的分布就是以给定的参数确定的，但我们在最后估计了均值和方差，
这个假设就被违反了，故而这个测试统计量的P值是含偏的，这个用法是错误的。

分布尾部

最后，我们可以检查分布的右尾，我们可以使用分位点函数ppf，其为cdf函数的逆，来获得临界值，
或者更直接的，我们可以使用残存函数的逆来办。

  1. >>> crit01, crit05, crit10 = stats.t.ppf([1-0.01, 1-0.05, 1-0.10], 10)
  2. >>> print 'critical values from ppf at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f'% (crit01, crit05, crit10)
  3. critical values from ppf at 1%, 5% and 10%   2.7638   1.8125   1.3722
  4. >>> 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))
  5. critical values from isf at 1%, 5% and 10%   2.7638   1.8125   1.3722
  7. >>> freq01 = np.sum(x>crit01) / float(n) * 100
  8. >>> freq05 = np.sum(x>crit05) / float(n) * 100
  9. >>> freq10 = np.sum(x>crit10) / float(n) * 100
  10. >>> print 'sample %%-frequency at 1%%, 5%% and 10%% tail %8.4f %8.4f %8.4f'% (freq01, freq05, freq10)
  11. sample %-frequency at 1%, 5% and 10% tail   1.4000   5.8000  10.5000

在这三种情况中，我们的样本有有一个更重的尾部，即实际在理论分界值右边的概率要高于理论值。
我们可以通过使用更大的样本来获得更好的拟合。在以下情况经验频率已经很接近理论概率了，
但即使我们重复这个过程若干次，波动依然会保持在这个程度。

  1. >>> freq05l = np.sum(stats.t.rvs(10, size=10000) > crit05) / 10000.0 * 100
  2. >>> print 'larger sample %%-frequency at 5%% tail %8.4f'% freq05l
  3. larger sample %-frequency at 5% tail   4.8000

我们也可以比较它与正态分布的尾部，其有一个轻的多的尾部：

  1. >>> print 'tail prob. of normal at 1%%, 5%% and 10%% %8.4f %8.4f %8.4f'% \
  2. ...       tuple(stats.norm.sf([crit01, crit05, crit10])*100)
  3. tail prob. of normal at 1%, 5% and 10%   0.2857   3.4957   8.5003

卡方检验可以被用来测试，是否一个有限的分类观测值频率与假定的理论概率分布具有显著差异。

  1. >>> quantiles = [0.0, 0.01, 0.05, 0.1, 1-0.10, 1-0.05, 1-0.01, 1.0]
  2. >>> crit = stats.t.ppf(quantiles, 10)
  3. >>> print crit
  4. [       -Inf -2.76376946 -1.81246112 -1.37218364  1.37218364  1.81246112
  5.   2.76376946         Inf]
  6. >>> n_sample = x.size
  7. >>> freqcount = np.histogram(x, bins=crit)[0]
  8. >>> tprob = np.diff(quantiles)
  9. >>> nprob = np.diff(stats.norm.cdf(crit))
  10. >>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
  11. >>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
  12. >>> print 'chisquare for t:      chi2 = %6.3f pvalue = %6.4f' % (tch, tpval)
  13. chisquare for t:      chi2 =  2.300 pvalue = 0.8901
  14. >>> print 'chisquare for normal: chi2 = %6.3f pvalue = %6.4f' % (nch, npval)
  15. chisquare for normal: chi2 = 64.605 pvalue = 0.0000

我们看到当t分布检验没被拒绝时标准正态分布却被完全拒绝。在我们的样本区分出这两个分布后，
我们可以先进行拟合确定scale与location再检查拟合后的分布的差异性。

我们可以先进行拟合，再用拟合分布而不是默认（起码location和scale是默认的）分布去进行检验。

  1. >>> tdof, tloc, tscale = stats.t.fit(x)
  2. >>> nloc, nscale = stats.norm.fit(x)
  3. >>> tprob = np.diff(stats.t.cdf(crit, tdof, loc=tloc, scale=tscale))
  4. >>> nprob = np.diff(stats.norm.cdf(crit, loc=nloc, scale=nscale))
  5. >>> tch, tpval = stats.chisquare(freqcount, tprob*n_sample)
  6. >>> nch, npval = stats.chisquare(freqcount, nprob*n_sample)
  7. >>> print 'chisquare for t:      chi2 = %6.3f pvalue = %6.4f' % (tch, tpval)
  8. chisquare for t:      chi2 =  1.577 pvalue = 0.9542
  9. >>> print 'chisquare for normal: chi2 = %6.3f pvalue = %6.4f' % (nch, npval)
  10. chisquare for normal: chi2 = 11.084 pvalue = 0.0858

在经过参数调整之后，我们仍然可以以5%水平拒绝正态分布假设。然而却以95%的p值显然的不能拒绝t分布。

正态分布的特殊检验

自从正态分布变为统计学中最常见的分布，就出现了大量的方法用来检验一个样本
是否可以被看成是来自正态分布的。

首先我们检验分布的峰度和偏度是否显著地与正态分布的对应值相差异。

  1. >>> print 'normal skewtest teststat = %6.3f pvalue = %6.4f' % stats.skewtest(x)
  2. normal skewtest teststat =  2.785 pvalue = 0.0054
  3. >>> print 'normal kurtosistest teststat = %6.3f pvalue = %6.4f' % stats.kurtosistest(x)
  4. normal kurtosistest teststat =  4.757 pvalue = 0.0000

将这两个检验组合起来的正态性检验

  1. >>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(x)
  2. normaltest teststat = 30.379 pvalue = 0.0000

在所有三个测试中，P值是非常低的，所以我们可以拒绝我们的样本的峰度与偏度与正态分布相同的假设。

当我们的样本标准化之后，我们依旧得到相同的结果。

  1. >>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % \
  2. ...                      stats.normaltest((x-x.mean())/x.std())
  3. normaltest teststat = 30.379 pvalue = 0.0000

因为正态性被很强的拒绝了，所以我们可以检查这种检验方式是否可以有效地作用到其他情况中。

  1. >>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(stats.t.rvs(10, size=100))
  2. normaltest teststat =  4.698 pvalue = 0.0955
  3. >>> print 'normaltest teststat = %6.3f pvalue = %6.4f' % stats.normaltest(stats.norm.rvs(size=1000))
  4. normaltest teststat =  0.613 pvalue = 0.7361

我们检验了小样本的t分布样本的观测值以及一个大样本的正态分布观测值，在这两种情况中我们
都不能拒绝其来自正态分布的空假设。得到这样的结果是因为前者是因为无法区分小样本下的t分布，
后者是因为它本来就来自正态分布。