分类数据的检验

必要的包

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from IPython.display import display
1
2
3
4
5

---

分类型数据（频数数据）

频数：落入某个特定分组的样本数量
当分组的维度只有1时，数据为单因素频数
维度为2时，数据为列联表数据

---

4.1 分类数据的常用检验

分类型数据的检验有：卡方检验，Fisher检验，McNermr检验，Cochran’s Q检验

卡方检验与Fisher检验是分类数据检验中最常见的检验

4.1.1 卡方检验和Fisher检验

• 两种检验的思想与卡方检验原理

卡方检验与Fisher检验都是分析频数表中绝对频数与期望频数的偏差程度

理论频数和期望频数是一样的，期望频数可以人为规定

检验的两个假设为:
$$H_0:\text{绝对频数与期望频数没有差别}\leftrightarrow H_1:\text{绝对频数与期望频数存在差别}$$

绝对频数：
在一个给定的频数表中，我们称其中的频数为绝对频数，记为$o_i,o_{ij}$,i,j为频数位置

期望频数：
根据进行卡方/Fisher检验的用途变化而变化，记为$e_i$

例如:
在列联表独立性检验中，每个格子的期望频数就是这个格子在表格中的“期望值”；
在特定分布的拟合优度检验中，期望频数就是特定分布在每个取值区间下的频数。

---

如果绝对频数与期望频数的差值越小，则两者越接近，我们越不能拒绝原假设

统计量
基于此卡方检验构造了下述检验统计量：
$$Teststatistics=\sum _i\frac{{\left(o_i-e_i\right)}^2}{e_i}$$
而该统计量近似服从卡方分布（在大样本下）：
$$Teststatistics\sim {\chi }_f^2$$

自由度f
其中，$f$为自由度。
对于有n个格子的单因素频数表而言，$f=n-1$；
对于$r\times c$的双因素列联表而言，$f=\left(r-1\right)\left(c-1\right)$

p值计算规则为：
$$pvalue=P\left({\chi }_f^2>Teststatistics\right)$$

卡方检验中的检验统计量是近似服从而非精确服从卡方分布
只有在大样本下（绝对频数与期望频数都很大）的情况下，卡方检验的精确度才高，而在小样本下，卡方检验的效用不及Fisher检验。相比于卡方检验这种“近似的”检验，Fisher检验是一种精确的检验，但是它的计算要比卡方检验复杂。

卡方检验的适用范围(注意绝对频数与理论频数是否过低)

卡方检验适用于单因素频数表、双因素频数表中的$2\times 2$与$r\times c$列联表。

单因素频数表

1. 在单因素频数表中，每一类的绝对频数$o_i$不能小于5。

$2\times 2$列联表
2. 样本总量$\sum o_i>40$，且所有期望频数$e_i>5$，可使用Pearson卡方检验
3. 样本总量$\sum o_i>40$，但存在期望频数$1<e_i<5$，可使用连续型校正的卡方检验
4. 若样本总量$\sum o_i<40$，或存在期望频数$1<e_i$，建议使用Fisher检验⚪

$r\times c$列联表
5. 表中期望频数$e_i<5$的格子不能超过1/5。
6. 不得出现期望频数$1<e_i$的情况。

Fisher检验的适用范围(2x2)

Fisher检验仅仅适用于双因素频数表中的$2\times 2$列联表⚪
Fisher检验在$2\times 2$列联表中适用范围内很广，弥补了卡方检验的缺点。

总结:
卡方检验在多种频数表中都可以应用，但要注意绝对频数与理论频数是否过低；
Fisher检验只适用于$2\times 2$列联表，但是在该表中的精确度与适用性都优于卡方检验。因此在$2\times 2$列联表中，我推荐大家使用Fisher检验；在其他表格中使用卡方检验。

4.2 单因素频数表的检验

接下来，我们开始使用卡方检验与Fisher检验，对特定的问题进行假设检验。
首先看单因素频数表的检验，常见的问题有：适度性检验与分布的拟合优度检验。

4.2.1 适度性检验

适度性检验的目的是：检验多项式实验得到的绝对频数与“期望值”是否存在显著差异。
在单因素频数表中，这个期望值就是表格的观测总数除以格子数：
$$expectedFrequency=\frac{\sum o_i}{i}$$
在前面的120次掷骰子的实验中，每个数字的期望值就是120/6=20。

Exp10

你和你的5个朋友们租了一个公寓一起生活。大家每天晚上都通过抽签决定某一个人去洗碗。经过几十天下来，你们统计了这段时间每个人的洗碗次数，如下图所示

你：20，小红：12，小明：10，李华：8，二狗：10，三猫：6

你发现有点不对劲，怎么你洗碗的次数一骑绝尘呢？你怀疑有人在作弊，想检验一下这个洗碗次数的分布是不是不太符合常理，你会怎么检验呢？

事实上，你想检验的本质上就是“每个人的实际洗碗次数与期望洗碗次数是否存在显著差异”，我们便可以使用适度性检验。

obs1=[20,12,10,8,10,6]
chi=stats.chisquare(obs1)
#stats.chisquare函数默认期望频数等于期望
# 若理论频数是期望，则我们只需要输入实际频数即可。
chi
#Power_divergenceResult(statistic=10.727272727272728, pvalue=0.057063677925693244)
#p值大于0.05,不能拒绝原假设
1
2
3
4
5
6
7

4.2.2 分布的拟合优度检验

单因素频数表检验的另外一个重要的应用场景就是分布的拟合优度检验

什么是分布拟合优度检验

设样本$x_i$是来自一个潜在总体$F(x)$,分布拟合优度检验所检验的问题是：
通过这些样本$x_i$判断潜在总体$F(x)$是否等同于一个假设的理论分布$F_0(x)$。原假设为
$$H_0:F\left(x\right)=F_0\left(x\right)$$

理论分布$F_0(x)$必须是一个分布函数已知的分布。
由于分布可以分为离散分布与连续分布，我们将分开讨论这两类分布的检验。

---

样本$x_i$来自离散分布

既然样本来自一个离散分布，那么理论分布$F_0(x)$也应当是一个离散分布。

设总体$X$为取有限个值$a_1$，$a_2$，……的离散随机变量。我们既可以将单独的$a_j$看成是一个类$A$，也可以将多个$a_j$合并成一个类$A$。例如，我们可以将“取值为1”的样本看成是一个类，也可以将“取值为1/2/3”的样本看成是一个类。我们将这有限个类记为$A_1$,……,$A_r$。

分类的唯一准则是：落入类A中的样本观测数不得少于5 ,否则就要与其他类合并。例如：取值为1的样本个数为3，取值为2的样本个数为4，那么我们就需要将它们合并为一个“取值为1和2”的类，因为这样才能保证落入该类的样本个数高于5。

记$P\left(X\in A_i\right)=p_i(i=1,2,\cdots ,r)$为理论分布$F_0(x)$在该分类下的“比例”，显然，我们用这个比例乘样本个数$n$，所得到的$n{p}_i$就是样本在这个分布下的期望频数$e_i$

如果样本的潜在总体$F(x)$服从理论总体$F_0(x)$，样本实际落入类$A_i$的实际频数$o_i$应当接近于期望频数$n{p}_i$，原假设便转换为了

$$H_0:o_i=n{p}_i$$

这恰巧就是卡方分布所检验的内容。于是我们便可以是用卡方分布进行分布拟合优度检验。

---

卡方分布进行拟合优度检验的步骤

总结一下使用卡方分布进行拟合优度检验的步骤：

1. 分类。根据实际频数$o_i$确定类别，若某个变量取值下样本的个数大于5，则直接将该取值作为一个单独的类；若小于5，则与相邻的取值合并为一个类。

2. 计算理论分布$F_0(x)$在该分类规则下，每个分类中的理论频数$n{p}_i$。

3. 进行卡方检验。

Exp11

我们验证卢瑟福实验的数据是否服从一个泊松分布。以下数据观测的是一枚放射性$\alpha$物质在单位时间内放射的质点数，实验重复2608次。

我们对该数据进行泊松分布的拟合优度检验。提示:我们用极大似然估计法计算出泊松分布参数$\lambda =3.87$

import numpy as np
import pandas as pd
from scipy import stats
#原始数据

data={‘counts': list(range(15)),'observe':[57,203,383,525,532,408,273,139,45,27,10,4,2,0,0]}
df=pd.DataFrame(data)
#将实际频数小于5的类别合并
df.iloc[11,1]=6
df=df[:12]
#Attention!!!
1
2
3
4
5
6
7
8
9
10
11

注意：这里的counts=11实际上是大于等于11的类，
因此在计算它的$p_i$时，不应当计算为$p_i=P\left(X=11\right)$，而应当是$p_i=P\left(X⩾11\right)=1-P\left(X<11\right)$

#根据自变量count的值计算每个自变量对应的理论频率
Poiss=stats.poisson(mu=3.87)
df['prop']=Poiss.pmf(df['counts'])
#pmf函数可以根据输入的自变量，输出对应的频率，也就是 理论概率 ）

#上述注意的修正
df.iloc[11,2]=1-Poiss.cdf(10)
#修正：由于数据框中counts=11实际上是大于等于11
#因此在这里修正counts大于11对应的概率
# cdf函数为左侧累积概率函数

# 用理论频率乘样本数，就可以得到理论频数
df['T_counts']=2608*df['prop']

#用卡方检验，比较实际频数与理论频数的差别
#就可以检验出数据是否服从泊松分布
chi=stats.chisquare(df['observe'],df['T_counts',ddof=1])
#ddof调整pvalue的自由度k - 1 - ddof
#其中k是观察到的频率数。默认值为ddof为 0。
#Power_divergenceResult(statistic=12.97449924325357, pvalue=0.2251010902807956)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

样本来自连续分布

样本为连续分布的情况要复杂一些，因为变量$X$的取值是不可数的，因此我们不能直接将某个取值作为一个类。一般采取的方法是：选$r-1$个实数$a_1<\cdots <a_{r-1}$，将实数组分为$r$个区间
$$(-\infty ,a_1],(a_1,a_2],\cdots ,(a_{r-1},\infty )$$
以这些区间为类，记录样本落入这些区间的个数，并计算这些区间在分布$F_0(x)$下的累积概率，往后的步骤就与前面离散分布检验的步骤一样了。

4.3 双因素列联表的检验

列联表检验常见的应用是两个特征的独立性检验与一致性检验

这两种检验虽然问题不同，但是检验的原理是完全一样的：
即比较绝对频数与表格中每个格子的**“期望值”**的差异。

在列联表中，每个单元格个期望值计算非常简单：
单元格所在行的行样本数乘单元格所在列的列样本数再除以样本总数
$$\text{ expectedFrequency }=\frac{\text{ RowTotal }\times \text{ ColumnTotal }}{n}$$

单元格男性右利手的期望值为52*87/100=45.2。

4.3.1 两个特征维度的独立性检验

独立性检验的是两个特征之间是否独立，一个变量作行，一个变量为列。

以上述的男女性/左右利手问题为例，我们可以检验性别特征与左右利手特征是否独立。

卡方检验

原假设独立，备择假设不独立

obs1=np.array([[43,9],[44,4]])
print(obs1)

#定义一个函数，使它可以输出连续性校正/未校正卡方检验的p值，以及期望频数表
def chiSquare(data):
	from scipy.stats import chi2_contingency
	chi2_corrected=stats.chi2_contingency(data, correction=True)
	#连续性校正
	chi_uncorrected=stats.chi2_contingency(data, correction=False)
	#非连续性校正
	
	print('Pearson卡方检验的p值为：{}'.format(chi2_uncorrected[1]))
    print('连续性校正的卡方检验的p值为：{}'.format(chi2_corrected[1]))
    print('-----------------------------------------------')
    print('期望频数表：')
    print(chi2_corrected[3])

#use
chiSquare(obs1)
Pearson卡方检验的p值为：0.1824670652605479
连续性校正的卡方检验的p值为：0.300384770390566
-----------------------------------------------
期望频数表：
[[45.24  6.76]
 [41.76  6.24]]
期望频数都大于5，应该使用非校正的卡方检验
p值约为0.3，无法拒绝原假设，即两种特征独立。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Fisher检验

使用Fisher检验解决上面的例子

def fisherExact(data,confidence=0.05):
	from scipy.stats import fisher_exact
	fisher_pvalue=fisher_exact(data)[1]

	if fisher_pvalue<confidence:
        print('在显著性水平{0:}下，可以拒绝原假设(p={1:.4f})'.format(confidence,fisher_pvalue))
	else:
		print('在显著性水平{0:}下，不可以拒绝原假设(p={1:.4f})'.format(confidence,fisher_pvalue))

fisherExact(obs1)
1
2
3
4
5
6
7
8
9
10

4.2 统一性检验

统一性检验的问题是：
两个或两个以上总体的某一特征分布，也就是各类别的比例是否统一或相近。

Example.12 某公司想了解南京和北京的市民对最低生活保障的满意程度是否相同，调查结果如下所示

在这个例子中，我们分析的是南京与北京这两个总体的一个特征“最低生活保障的满意程度”分布是否统一。

卡方检验

obs2=np.array([[100,110],[150,160],[180,170],[170,160]])
print(obs2)
chiSquare(obs2)
1
2
3

作业

一个淘宝网购商家搜集了一年中每天的订单数$X$，除去春节期间及双十一前后外，按330天记，数据如下

请用卡方分布验证订单数是否泊松分布。已知：通过极大似然估计得知泊松分布参数$\lambda =5.3$

data=pd.DataFrame({'订单数': list(range(14)) + [16],
    '天数': [3, 6, 21, 46, 48, 61, 52, 42, 27, 11, 6, 4, 1, 1, 1]})
data
1
2
3

#类别数量小于5的合并
#
data.loc[1,'天数']=9
data.loc[11,'天数']=7
data=data[1:12]
data

1
2
3
4
5
6
7
8

#计算理论频率,据题意mu=5.3
Poiss=stats.poisson(mu=5.3)
data['prop']=Poiss.pmf(data['订单数'])

#修正
data.loc[11,'prop']=1-Poiss.cdf(10)
#cdf函数为左侧累积概率函数
data.loc[1,'prop']=Poiss.cdf(1)
#计算理论概数
data['T_counts']=data['天数'].sum()*data['prop']
data
1
2
3
4
5
6
7
8
9
10
11

#卡方分布检验
chi = stats.chisquare(data['天数'], data['T_counts'], ddof=1)
chi
Power_divergenceResult(statistic=3.9705897417232916, pvalue=0.9133375422858901)
#p>confidence因此可以认为样本总体服从泊松分布
1
2
3
4
5

教程来源:https://github.com/Git-Model/Modeling-Universe/tree/main/Data-Story