Python中的卡方检验

Question

我在中使用了以下代码来判断观察到的数值（比如20, 20, 0和0）与预期的数值/比例（比如每种情况各占25%）的契合度：

> chisq.test(c(20,20,0,0), p=c(0.25, 0.25, 0.25, 0.25))

    Chi-squared test for given probabilities

data:  c(20, 20, 0, 0)

X-squared = 40, df = 3, p-value = 1.066e-08

我该如何在Python中实现这个功能呢？我尝试使用scipy中的chisquare函数，但得到的结果差别很大；我不确定这是否是正确的函数。我查阅了scipy的文档，但那内容实在太多了，有1000多页；而numpy的文档几乎还要多50%。

Answer 1

另一种方法是从Python中调用你的R代码。你可以这样做：

• 通过将R脚本作为命令行工具运行。想了解更多关于如何使用 Rscript 从命令行运行R脚本的信息，可以查看这个链接。然后在Python中，你可以通过使用 subprocess 或 os.system 来执行系统调用，从而运行R脚本。数据交换是通过文本文件或二进制文件进行的。我喜欢这种方法，因为它非常简单，而且可以很容易地单独调试R脚本，而不影响Python代码。缺点是所有数据都要经过硬盘，这可能会很慢。
• 通过使用rpy 或 rpy2，可以直接在Python中运行R代码。这样集成得更紧密，但这个链接也有自己的一些小问题。例如，根据我的经验，通过rpy调用的R代码调试起来会稍微困难一些。

Answer 2

我想指出的是，虽然这个答案在语法上看起来是正确的，但在你的例子中不应该使用卡方分布，因为你观察到的频率太小，无法进行准确的卡方检验。

“当每个类别中的观察频率或预期频率太小时，这个检验就是无效的。一个常见的规则是，所有观察到的和预期的频率都应该至少为5。” 你可以查看这个链接了解更多信息： http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare

Answer 3

scipy.stats.chisquare 这个函数需要你提供实际观察到的频率和预期的绝对频率，而不是比例。你可以用下面的代码来得到你想要的结果：

>>> observed = np.array([20., 20., 0., 0.])
>>> expected = np.array([.25, .25, .25, .25]) * np.sum(observed)
>>> chisquare(observed, expected)
(40.0, 1.065509033425585e-08)

不过，如果你预期的值在各个类别中是均匀分布的，那你可以不计算预期值：

>>> chisquare(observed)
(40.0, 1.065509033425585e-08)

这个函数返回的第一个值是 χ² 统计量，第二个值是测试的 p 值。