Python 数学建模——独立性检验

前言

  独立性检验是概率论中比较重要的内容，常常用于检验两个变量是否具有关联关系。下面结合代码说说数学建模中如何使用 Python 进行独立性检验。

卡方独立性检验

原理

  ${\chi }^2$ （卡方）独立性检验是一种基于列联表的检验，用于分析两种属性之间的独立性。假设一共有 $N$ 个样本，每个样本都有 $X,Y$ 两个属性，属性的取值集合为 { X 1 , X 2 , ⋯ , X n } × { Y 1 , Y 2 , ⋯ , Y m } \{{X_1},{X_2},\cdots ,{X_n}\}\times \{{Y_1},{Y_2},\cdots ,{Y_m}\} {X1​,X2​,⋯,Xn​}×{Y1​,Y2​,⋯,Ym​}。其中：

• 属性组合为 $(X_i,Y_j)$ 的样本的频数记为 $f_{ij}$。
• $X$ 属性为 $X_i$ 的样本总数记为 $N_{X_i}$，$Y$ 属性为 $Y_j$ 的样本总数记为 $N_{Y_j}$。

  那么独立性检验的步骤就是：

1. 建立原假设 $H_0$：两变量相互独立。备择假设 $H_1$：两变量不相互独立。
2. 计算理论频数 $e_{ij}=N_{X_i}\times N_{Y_j}/N$。
3. 计算自由度 $df=(n-1)(m-1)$。
4. 计算卡方统计量 ${\chi }^2=\sum _{i=1}^n\sum _{j=1}^m\frac{(f_{ij}-e_{ij}{)}^2}{e_{ij}}\sim {\chi }^2(df)$。

实际观察次数与理论次数之差的平方再除以理论次数得到的统计量近似服从卡方分布。

查 ${\chi }^2$ 分布表，确定满足式 ${\chi }_{1-\alpha /2}^2(df)<{\chi }^2<{\chi }_{\alpha /2}^2(df)$ 的最大 $\alpha$ 值（这个值记作 $p$ 值）。如果 $p>0.05$，就可以接受原假设；否则拒绝原假设，接受备择假设。
参考文献：列联表分析——独立性检验（卡方检验）_列联表卡方检验-CSDN博客

注意使用条件：

• 总样本数大于 $40$ 且每个类别的理论频数大于等于 $5$。
• 自由度 $df>1$。

  不符合上面的第一条使用条件时，采用 Fisher 精确检验。
  不符合上面的第二条使用条件时，需要使用 Yates 矫正的卡方检验（实际代码会根据 $df$ 自动判断是否矫正，但是论文中应说明使用的是矫正的卡方检验）。

代码实例

  验证“性别”和“信来世”是否独立。联合分布频数表如下：

	信来世$(Y_1)$	不信来世$(Y_2)$	合计
女性$(X_1)$	$509$	$116$	$625$
男性$(X_2)$	$398$	$104$	$502$
合计	$907$	$220$	$1127$

import numpy as np
from scipy.stats import chi2_contingency# 输入列联表 table
table = np.array([[509,116],[398,104]])res = chi2_contingency(table,correction = False)
print(res)
# 结果 卡方值 = res[0] = 0.8245764245493725,p_value = res[1] = 0.3638455155911725

  这个 $p=0.36>0.05$，所以是接受原假设，认为“性别”和“信来世”是相互独立的。设置参数correction = False的目的是关闭 Yates 矫正，使得计算结果与手算相同。

注：使用chi2_contingency进行卡方检验时，如果自由度 $df=(n-1)(m-1)=1$，也即 $m=n=2$，则底层代码会自动进行 Yates 校正。需要手动关闭 correction 才可以使得结果和手算相同。
参考文献：python用chi2_contingency做卡方检验结果和手算以及spss不同的问题_spss和python 2x2卡方值不同-CSDN博客

Fisher 精确检验

  对于不满足 ${\chi }^2$ 独立检验条件的，可以采用 Fisher 独立性检验，即适用范围为：总样本数小于$40$ 或者有一个类别的理论频数小于 $5$。
  Python 中内置有 Fisher 精确检验库（scipy.stats.fisher_exact），但是该库只能进行 $2\times 2$ 列联表的检验，应用场景极其狭窄。通过 Python 引用 R 语言包（需要pip install rpy2）可以打破 $m=n=2$ 的限制：

import rpy2.robjects as ro
import numpy as np# 创建一个列联表
matrix_data = np.array([[11,11],[6,0],[17,13]])# 转换为 R 矩阵格式
r_matrix = ro.r.matrix(ro.IntVector(matrix_data.flatten()),nrow = matrix_data.shape[0],byrow=True)# 执行多维 Fisher 精确检验
fisher_test = ro.r['fisher.test']
result = fisher_test(r_matrix)
print(result)
"""Fisher's Exact Test for Count Datadata:  structure(c(11L, 6L, 17L, 11L, 0L, 13L), dim = 3:2)
p-value = 0.08361
alternative hypothesis: two.sided
"""