3.4.6 卡方分布和卡方检验
卡方分布（χ2）是概率论与统计学中常用的一种概率分布。k个独立的标准正态分布变量的平方和服从自由度为k的卡方分布。下面通过随机数验证该分布，结果如图3-21所示：
图3-21 使用随机数验证卡方分布
a = np.random.normal(size=(300000, 4))
cs = np.sum(a**
2, axis=1)
sample_dist, bins = np.histogram(cs, bins=100, range=(0, 20), density=True)
x = 0.5 *
(bins[:-1] + bins[1:])
chi2_dist = stats.chi2.pdf(x, 4)
print "max error:", np.max(np.abs(sample_dist - chi2_dist))
max error: 0.00340194486328
卡方分布可以用来描述这样的概率现象：袋子里有5种颜色的球，抽到每种球的概率相同，从中选N次，并统计每种颜色的次数Oi。则下面的χ2符合自由度为4的卡方分布，其中E=N/5为每种球被抽选的期望次数：
下面用程序模拟这个过程，结果如图3-22所示。❶首先调用randint( )创建从0到5的随机数，其结果ball_ids的第0轴表示实验次数，第1轴为每次实验抽取的100个球的编号。❷使用bincount( )统计每次实验中每个编号出现的次数，由于它不支持多维数组，因此这里使用apply_along_axis( )对第1轴上的数据循环调用bincount( )。为了保证每行的统计结果的长度相同，设置minlength参数为5，apply_along_axis( )会将所有关键字参数传递给进行实际运算的函数。❸使用上面的公式计算χ2统计量cs2，❹并用gaussian_kde( )计算cs2的分布情况。
图3-22 模拟卡方分布
repeat_count = 60000
n, k = 100, 5
np.random.seed(42)
ball_ids = np.random.randint(0, k, size=(repeat_count, n)) ❶
counts = np.apply_along_axis(np.bincount, 1, ball_ids, minlength=k) ❷
cs2 = np.sum((counts - n/k)**
2.0/(n/k), axis=1) ❸
k = stats.kde.gaussian_kde(cs2) ❹
x = np.linspace(0, 10, 200)
pl.plot(x, stats.chi2.pdf(x, 4), lw=2, label=u"$\chi ^{2}$分布")
pl.plot(x, k(x), lw=2, color="red", alpha=0.6, label=u"样本分布")
pl.legend(loc="best")
pl.xlim(0, 10)
卡方检验可以用来评估观测值与理论值的差异是否只是因为随机误差造成的。在下面的例子中，袋子中各种颜色的球按照probabilities参数指定的概率分布，choose_balls(probabilities, size)从袋中选择size次并返回每种球被选中的次数。袋子1中的球的概率分布为：0.18、0.24、0.25、0.16、0.17。袋子2中各种颜色的球的个数一样多。通过调用choose_balls( )得到两组数字：80 93 97 64 66和89 76 79 71 85。现在需要判断袋子中的球是否是平均分布的。
def choose_balls(probabilities, size):
r = stats.rv_discrete(values=(range(len(probabilities)), probabilities))
s = r.rvs(size=size)
counts = np.bincount(s)
return counts
np.random.seed(42)
counts1 = choose_balls([0.18, 0.24, 0.25, 0.16, 0.17], 400)
counts2 = choose_balls([0.2]*
5, 400)
counts1 counts2
-------------------- --------------------
[80, 93, 97, 64, 66] [89, 76, 79, 71, 85]
使用chisquare( )进行卡方检验，它的参数为每种球被选中次数的列表，如果没有设置检验的目标概率，就测试它们是否符合平均分布。卡方检验的零假设为样本符合目标概率，由下面的检验结果可知，第一个袋子对应的p值只有0.02，也就是说如果第一个袋子中的球真的符合平均分布，那么得到的观测结果80 93 97 64 66的概率只有2%，因此可以推翻零假设，即袋子中的球不太可能是平均分布的。第二个袋子对应的p值为0.64，无法推翻零假设，即我们的结论是不能否定第二个袋子中的球是平均分布的。注意，和前面介绍的t检验一样，零假设只能用来否定，因此不能根据观测结果89 76 79 71 85得出袋子中的球是符合平均分布的结论。
chi1, p1 = stats.chisquare(counts1)
chi2, p2 = stats.chisquare(counts2)
print "chi1 =", chi1, "p1 =", p1
print "chi2 =", chi2, "p2 =", p2
chi1 = 11.375 p1 = 0.0226576012398
chi2 = 2.55 p2 = 0.635705452704
卡方检验是通过卡方分布进行计算的。图3-23显示自由度为4的卡方分布的概率密度函数，以及chi1和chi2对应的位置
和
。p1是
右侧部分的面积，而p2是
右侧的面积。
图3-23 卡方检验计算的概率为阴影部分的面积
卡方检验还可以用于二维数据。前面介绍的彩色球的例子中，只是按照球的颜色分组，而二维数据则按照样本的两个属性分组，统计学上称之为列联表。例如下面是性别与惯用手的列联表，我们希望知道这个统计结果能否说明性别与惯用手之间存在某种联系。
下面使用chi2_contingency( )对列联表进行卡方检验，零假设为性别与惯用手之间不存在联系，即男性与女性惯用左右手的概率相同。由于p值为0.3，因此不能推翻零假设，即该实验数据中没有明显证据表明男性和女性在使用左右手习惯上存在区别。
table = [[43, 9], [44, 4]]
chi2, p, dof, expected = stats.chi2_contingency(table)
chi2 p
------------------ -------------------
1.0724852071005921 0.30038477039056899
对于上面的2×2的数值较小的列联表，可以使用fisher_exact( )计算出精确的p值：
stats.fisher_exact(table)
(0.43434343434343436, 0.23915695682225618)