Python中的卡方检验

Question

我在R 中使用了以下代码来确定观察值（例如 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

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 值。

Answer 2

只是想指出，虽然答案在句法上似乎是正确的，但您不应该在示例中使用卡方分布，因为您观察到频率太小而无法进行准确的卡方检验。

“当每个类别中观察到的或预期的频率太小时，此测试无效。一个典型的规则是所有观察到的和预期的频率至少应为 5。”看： http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chisquare.html#scipy.stats.chisquare

Answer 3

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

• 通过使 R 脚本作为命令行工具运行。有关使用Rscript 从命令行运行 R 脚本的更多信息，请参阅this link。然后，您可以在 python 中通过使用subprocess 或os.system 执行系统调用来运行R 脚本。任何数据交换都是通过文本或二进制文件完成的。我喜欢这种方法，因为它非常简单，并且很容易将 R 脚本与 python 代码分开调试。缺点是所有数据都通过硬盘驱动器，这可能会非常缓慢。
• 通过使用rpy 或rpy2 直接从python 中运行R 代码。这样集成更加紧密，但是这个链接也引入了自己的小怪癖。例如，根据我的经验，调试通过 rpy 调用的 R 代码有点难调试。