python中独立性检验 独立性检验函数

在单样本问题中, 人们想要检验的是总体的中心是否等于一个已知的值. 但在实际问题中, 更受注意的往往是比较两个总体的位置参数; 比如, 两种训练方法中哪一种更出成绩, 两种汽油中哪一种污染更少, 两种市场营销策略中哪种更有效等等.

1. $\chi^2$独立性检验的原理

若随机变量$X,Y$的分布函数分别为$F_1(x) 和 F_2(y)$, 且联合分布为$F(x, y)$, 则X与Y的独立性归结为假设检验问题：

$H_0: F(x,y)=F_1(x)F_2(y) \quad H_1: F(x,y)\neq F_1(x)F_2(y)$

若X与Y为分类变量，其中X的取值为$X_1, X_2,...,X_r$, Y的取值为$Y_1,Y_2,...,Y_s$, 将X与Y的各种情况组合用一张$r \times s$列联表表示，称为$r\times s$二维列联表，如下图所示：

表中$n_{ij}$表示n个随机试验的结果中X取$X_i$及Y取$Y_j$的频数，$\sum_{i=1}^r\sum_{j=1}^sn_{ij}=n$.

$n_{i.}=\sum_{j=1}^sn_{ij}, i=1,2,...,r, 表示各行之和$

$n_{.j}=\sum_{i=1}^rn_{ij}, j=1,2,...s, 表示各列之和$

令$p_{ij}=P(X=X_i, Y=Y_j), p_{i.}=P(X=X_i), p_{.j}=P(Y=Y_j), i=1,2,...,r; j=1,2,...,s,$ 则X与Y的独立性检验等价于下述检验：

$H_0: p_{ij} = p_{i.}p_{.j}, \forall 1\le i \le r, 1\le j \le s \quad H_1: \exist (i,j), p_{ij}\neq p_{i.}p_{.j}$

注: 若X与Y 为连续型随机变量, 这时将它们的取值范围分成r个及s个互不相交的小区间, 用$n_{ij}$ 表示n个随机试验的结果中“X属于第i个小区间, Y 属于第k个小区间”的频数 $(i=1,2,...,r; k=1,2,...,s)$ 这时可将X与Y 的独立性转化为列联表的独立性检验问题.

由于$p_{i.}$的极大似然估计为$\hat p_{i.} = n_{i.}/n$, $p_{.j}$的极大似然估计为$\hat p_{.j}/n$, 因此若$H_0$成立，则$p_{ij}$的极大似然估计为$\hat p_{i.} \hat p_{.j}=n_{i.}n_{.j}/n^2$. 从而X取代$X_i$, Y取代$Y_j$(试验数据落入第$(i,j)$个类)的理论频数为$n \times n_{i.}n_{.j}/n^2 = n_{i.}n_{.j}/n$. 由此构造检验统计量为：

$\chi^2=\sum_{i=1}^r\sum_{k=1}^s[n_{ij}-\frac{n_{i.}n_{.j}}{n}]^2/\frac{n_{i.}n_{.j}}{n}$

可以证明在原假设成立时， $\chi^2$近似服从$\chi^2((r-1)(s-1))$

2. 使用Python展示和计算

在python中的pandas库中有一个DataFrame的数据结构，可以用来保存二维数据，并进行计算。

例1： 对表1所示频数分布表，以59%显著水平，检验色觉与性别是否相关。

2. 1. 初始化数据

import numpy as np
import pandas as pd
from scipy import stats
d=np.array([[442, 514],[38,6]])
r,s = d.shape
df = pd.DataFrame(d, columns=['male', 'female'], index=['normal', 'blindness'])

上述代码已将初始化数据存入df，显示如下：

2.2. 计算边际频数

首先在df基础上建立一个数据框（数据一致）。并在新建的数据框基础继续计算。

df1 = df[:]
df1= pd.DataFrame(df1)
df1['r_tot']=np.sum(df1,axis=1)
df1.loc['c_tot']=np.sum(df1, axis=0)

经过上面的处理，df1就是带有编辑频数的数据框。显示如下：

2.3. 计算理论频数分布表

r_tot = df1['r_tot'][:-1]
c_tot=df1.loc['c_tot'][:-1]
total = np.sum(r_tot)
data = np.zeros((r,s))
for i in range(len(r_tot)):
    for j in range(len(c_tot)):
        data[i,j]=r_tot[i]*c_tot[j]/total
df2=pd.DataFrame(data, index=['normal', 'blindness'], columns=['male', 'female'])

至此，我们将计算得到的理论频数分布表保存在了df2中，而原始数据在df中，如下图所示：

2.4. 统计量计算

为了清晰整个计算过程，这里添加一步（对计算结果来说是多余的）产生统计量数据框的代码：

(df-df2)**2/df2

其产生的统计量数据框如下图：

其实，下面一行代码即可获得卡方统计量的值并获取pval：

chi_square=np.sum((df-df2)**2/df2).sum()
stats.chi2.sf(chi_squre, df=(r-1)*(s-1))

得到$\chi^2=27.13874340443378, pval=1.8936459120177876e-07$
因为$pval \lt \alpha=0.05$, 所以拒绝原假设，可以认为色盲和性别有关系