Python进行假设检验

网上关于假设检验的内容太多了，说法也是各式各样。在整理这篇博客时参考的资料都在末尾中给出了，大神们可以自行去看。
这篇博客就是个小心求证的过程！！！

1 中心极限定理

  中心极限定理指出，如果样本量足够大，则变量均值的采样分布将近似于正态分布，而与该变量在总体中的分布无关。
  举个例子，以$X$表示10次抛硬币试验中正面出现的次数，则很明显随机变量$X$服从二项分布。现在我们从$X$取值的总体中抽样出60次试验的结果，并根据这60次试验的结果计算随机变量$X$的均值(即采样均值$\overline{X}$)，中心极限定理的意思就是说：采样均值$\overline{X}$近似服从正态分布。为了验证这个定理，下面我们用python模拟了1000次$\overline{X}$的均值计算，每次抽取60个样本计算$\overline{X}$，其采样均值$\overline{X}$分布结果如下图。

import pandas as pd
import numpy as np
from scipy import stats
import random
from scipy.special import comb
import seaborn as sns
##10次抛硬币的二项分布正面出现的次数取值为0,1,……,10,计算每个取值出现的概率
#prob为k_prob中index所指示的取值出现的概率
#prob_sum为累计概率，主要是为了方便后面的轮盘赌算法
k_prob=[comb(10,i)*(0.5**10) for i in range(0,11)]
k_prob=pd.DataFrame(k_prob,columns=['prob'])
k_prob['prob_sum']=k_prob['prob'].cumsum()
#轮盘赌
x_mean=[]
for _ in range(0,1000):
    sample=[] #从总体中取60个样本
    for i in range(60):
        p=random.random()
        sample.append(k_prob[k_prob['prob_sum']>=p]['prob_sum'].idxmin())
    x_mean.append(sum(sample)/len(sample))
#画图        
sns.distplot(x_mean,kde=False,fit=stats.gamma)#曲线做了fit拟合

再假设随机变量$X$服从0-1均匀分布，其1000次采样均值$\overline{X}$的分布结果如下图。

##0-1均匀分布计算的均值采样分布
num=list(range(1,1001))
x_mean=[]
for k in num:
    x=[random.random() for _ in range(60)]
    x_mean.append(sum(x)/len(x))
sns.distplot(x_mean,kde=False,fit=stats.gamma)

2 假设检验

2.1 基本思想

  假设检验，是用来判断样本与样本、样本与总体的差异是由抽样误差引起的还是本质差别造成的统计推断方法。其基本思想是小概率和反证法。首先提出一个假设，用适当的统计方法确定当这个假设成立时，获得现在样本的概率大小。如果获得当前样本的概率极小，则认为在当前这个假设下获得当前采样值是一件不可能的事。据此则认为该假设不成立，否则接受该假设。

2.2 原假设和备择假设

  在假设检验中，首先提出两种假设，即原假设和备择假设。

• 原假设：记为$H_0$。该假设认为样本与总体或样本与样本之间的差异是由抽样误差引起的。原假设一般是统计者想要拒绝的假设。
• 备择假设：记为$H_1$。该假设认为样本与总体或样本与样本之间存在本质差异。备择假设一般是统计者想要接受的假设。

  设总体参数$\mu$的假设值为${\mu }_0$。对总体参数$\mu$的假设检验一定要采取下面三种形式之一：
  (1) $H_0:\mu \ge {\mu }_0,H_1:\mu <{\mu }_0$ (单侧检验)
  (2) $H_0:\mu \le {\mu }_0,H_1:\mu >{\mu }_0$ (单侧检验)
  (3) $H_0:\mu ={\mu }_0,H_1:\mu \ne {\mu }_0$ (双侧检验)

2.3 检验统计量

  确定假设之后，一般不直接利用样本观测数据来建立检验准则，因为原始的样本观测数据的分布可能是未知的。所以需要寻找一个抽样分布可以确定或近似确定的统计量来计算，然后将样本数据代入到统计量的公式中。不同的假设检验对应不同的统计量，具体分为以下几种。

2.3.1 总体均值的检验

  这类检验的目的是验证总体均值$\overline{X}$是否等于(或大于，或小于)${\mu }_0$。依据总体方差和样本总数的不同情况，其常用统计量如下表:

情况	检验统计量
总体服从正态分布且方差已知	$Z=\frac{\overline{x}-{\mu }_0}{\sigma /\sqrt{n}}$，其中$n$为抽样样本数量，$\overline{x}$为样本均值(下同)，$\sigma$为总体分布的标准差，统计量$Z$服从标准正态分布(中心极限定理)
总体分布及其方差均未知但为大样本($n>30$)	$Z=\frac{\overline{x}-{\mu }_0}{s/\sqrt{n}}$,其中样本方差$s=\frac{{\sum }_{i=1}^n(x_i-\overline{x}{)}^2}{n-1}$(下同)，统计量$Z$服从标准正态分布(中心极限定理)
总体服从正态分布，方差未知且为小样本	$t=\frac{\overline{x}-{\mu }_0}{s/\sqrt{n}}$，统计量$t$服从自由度为$n-1$的$t$分布

2.3.2 两个均值之差的检验

  设两个总体的均值分别为${\overline{X}}_1$和${\overline{X}}_2$,两个总体的方差分别为$S_1^2$和$S_2^2$，来自两个总体的样本容量分别$n_1$和$n_2$，样本均值分别为${\overline{x}}_1$和${\overline{x}}_2$。检验的目的是两个总体的均值是否相等，或两个总体的均值之差是否为零。依据总体方差和样本总数的不同情况，其常用统计量如下表:

情况	检验统计量
两个总体服从正态分布且方差已知	$Z=\frac{{\overline{x}}_1-{\overline{x}}_2}{\sqrt{\frac{S_1^2}{n_1}+\frac{S_2^2}{n_2}}}$，${\overline{x}}_1$和${\overline{x}}_2$分别为两个总体的样本均值(下同)，$Z$服从标准正态分布
两个总体方差未知且不相等但为大样本	$Z=\frac{{\overline{x}}_1-{\overline{x}}_2}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}$,$s_1^2$和$s_2^2$为样本方差(下同)，$Z$服从标准正态分布
两个总体方差未知且相等但为大样本	$Z=\frac{{\overline{x}}_1-{\overline{x}}_2}{\sqrt{\frac{s_1^2}{n_2}+\frac{s_2^2}{n_1}}}$，$Z$服从标准正态分布
两个总体服从正态分布，方差未知且相等但为小样本	$t=\frac{{\overline{x}}_1-{\overline{x}}_2}{s\sqrt{\frac{1}{n_1}+\frac{1}{n_2}}}$,$s=\sqrt{\frac{(n_1-1)s_1^2+(n_2-1)s_2^2}{n_1+n_2-2}}$，$t$服从自由度为$n_1+n_2-2$的$t$分布
两个总体服从正态分布，方差未知且不相等但为小样本	$t=\frac{({\overline{x}}_1-{\overline{x}}_2)}{\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}}$，$t$服从自由度为$d{f}^{\prime }$的$t$分布，$d{f}^{\prime }=\frac{(\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}{)}^2}{\frac{(\frac{s_1^2}{n_1}{)}^2}{n_1-1}+\frac{(\frac{s_2^2}{n_2}{)}^2}{n_2-1}}$

2.3.3 总体成数(比例)的检验

  检验的目的是判断总体成数$P$是否等于$P_0$。其检验统计量如下：
$$Z=\frac{p-P_0}{\sqrt{\frac{P_0(1-P_0)}{n}}}$$其中$p$为样本成数，该统计量服从标准正态分布。

2.3.4 总体成数之间的检验

  设两个总体成数分别为$P_1$和$P_2$,来自两个总体的样本容量分别为$n_1$和$n_2$，样本成数分别为$p_1$和$p_2$。检验两个总体成数是否相等，或两个总体成数之差是否相等。当$n_1$和$n_2$和都足够大时($n_1{P}_1$，$n_1{P}_1(1-P_1)$，$n_2{P}_2$，$n_2{P}_2(1-P_2)$都大于5)，其检验统计量如下：
$$Z=\frac{(p_1-p_2)}{\sqrt{\hat{P}(1-\hat{P})(\frac{1}{n_1}+\frac{1}{n_2}}）}$$其中$$\hat{P}=\frac{n_1{p}_1+n_2{p}_2}{n_1+n_2}$$检验统计量$Z$近似服从正态分布。
  当检验两个总体成数之差为某一非零的常数的假设，其检验统计量如下：
$$Z=\frac{(p_1-p_2)-d_0}{\sqrt{\frac{p_1(1-p_1)}{n_1}+\frac{p_2(1-p_2)}{n_2}}}$$

2.3.5 总体方差的检验

  检验目的是判断正态总体方差$S^2$是否等于$S_0^2$。其检验量为:
$${\chi }^2=\frac{(n-1)s^2}{S_0^2}$$其中$s^2$为样本方差，该检验量服从自由度为$n-1$的${\chi }^2$分布(卡方分布)。

2.3.6 总体方差之比的检验

  设两个总体方差$S_1^2$和$S_2^2$，相应的样本方差分别为$s_1^2$和$s_2^2$，样本容量分别为$n_1$和$n_2$，该检验是为了判断两个总体的方差是否相等。其检验统计量为：
$$F=\frac{s_1^2}{s_2^2}$$该统计量服从分子自由度为$n_1-1$,分母自由度为$n_2-1$的$F$分布。

2.4 显著性水平

  在假设检验之前，需要设定一个概率$\alpha$，这个概率即显著性水平。通常$\alpha$取0.05或0.01，有时也取0.1。将样本代入到检验统计量中，计算出检验统计值，然后求出统计量分布曲线中与检验统计值相对应的，称之为观测到的显著性水平$P$-值。将$P$-值与事先确定好的显著性水平$\alpha$值相比较来判断是否拒绝原假设。原假设类型不同判断标准也不同。

• 在双侧检验中，若$P<\alpha /2$则拒绝原假设，否则不拒绝原假设。
• 在单侧检验中，若$P<\alpha$则拒绝原假设，否则不拒绝原假设。

  假设检验还可以用临界值来判断，但python中的函数可以直接返回$P$值，使用非常方便。(所以就不介绍临界值规则了，网上资料很多)

3 python代码实现

  (1)以参考资料4中作者提供的数据(作者提供了下载链接)为例，构造几个问题并用python进行假设检验。其代码如下：

import pandas as pd
import numpy as np
from scipy import stats
import statsmodels.stats.weightstats as ws

path=r'/Users/sherry/Downloads' #换成自己的文件地址
data=pd.read_csv(path+r'/test.csv')
#问题 判断人体体温是否服从正态分布
#问题1 判断人体体温均值是否为98.6
#问题2 从数据中随机抽取25条数据，使用t检验验证问题1
#问题3 判断男性和女性之间体温均值是否相等
#问题4 分别从男性和女性体温数据中抽取25条数据，验证男性和女性体温是否相等

#问题0 人体体温是否服从正态分布
stats.kstest(data['Temperature'],'norm',(data['Temperature'].mean(),data['Temperature'].std()))
#这里测试的pval=0.62 这说明人体体温是服从正态分布的(这和参考资料中的结论不一致)


#问题1 判断人体体温均值是否为98.6
#原假设 H0：体温均值=98.6
# 三种实现方法
zstat,pval=ws.ztest(x1=data['Temperature'],value=98.6,alternative='two-sided',ddof=1)

ztest=ws.DescrStatsW(data['Temperature'])
zstat,pval=ztest.ztest_mean(value=98.6,alternative='two-sided')

x1_mean=data['Temperature'].mean()
x1_std=np.std(data['Temperature'],ddof=1)
zstat=(x1_mean-98.6)/(x1_std/np.sqrt(data.shape[0]))
pval=stats.norm.sf(np.abs(zstat))*2 #残存函数=1-累计分布函数

#zstat=-5.45 p=4.9e-8<0.05 拒绝原假设
##从ztest的源码中发现,ztest中没有办法处理方差已知情况,ztest函数源码中已经默认使用样本标准差进行计算了
#ztest直接返回依据样本计算的统计量值和在原假设成立时获得现在的值的概率


#问题2 从数据中随机抽取25条数据，使用t检验验证问题1(总体均值t检验)
data_index=[86,85,51,127,16,13,93,24,111,91,115,81,27,122,116,36,104,2,95,49,98,
            1,17,41,20] #抽取的样本索引
data_t=data.loc[data_index]

#原假设H0:体温均值=98.6
#三种实现方法
x1_mean=data_t['Temperature'].mean()
x1_std=np.std(data_t['Temperature'],ddof=1)
tstat=(x1_mean-98.6)/(x1_std/np.sqrt(data_t.shape[0]))
pval=stats.t.sf(np.abs(tstat),df=data_t.shape[0]-1)*2

ttest=ws.DescrStatsW(data_t['Temperature'])
tstat,pval,_=ttest.ttest_mean(value=98.6,alternative='two-sided')

tstat,pval=stats.ttest_1samp(data_t['Temperature'],popmean=98.6)
#tstat=-2.81 pval=0.009<0.05 拒绝原假设


#问题3 判断男性和女性之间的体温均值否相等(总体均值之差的检验)
#原假设H0:男女体温均值相等
data1=data[data['Gender']==1] #男性，65条数据
data2=data[data['Gender']==2] #女性，65条数据

#(1)先验证男性和女性体温方差是否相等(F检验)
#原假设H0：男性体温方差=女性体温方差
d1=np.var(data1['Temperature'],ddof=1)
d2=np.var(data2['Temperature'],ddof=1)

fstat=d1/d2
#F分布不是对称分布
pval=1-2*abs(0.5-stats.f.cdf(fstat,data1.shape[0]-1,data2.shape[0]-1))

#若男性、女性体温方差相等 
d1=ws.DescrStatsW(data1['Temperature'])
d2=ws.DescrStatsW(data2['Temperature'])
com=ws.CompareMeans(d1,d2)
zstat,pval=com.ztest_ind(alternative='two-sided',usevar='pooled',value=0)
#zstat=-2.28 pval=0.022 拒绝原假设
#若方差不相等,usevar='unequal'

x1_mean=np.mean(data1['Temperature'])
x2_mean=np.mean(data2['Temperature'])
x1_var=np.var(data1['Temperature'],ddof=1)
x2_var=np.var(data2['Temperature'],ddof=1)
zstat=(x1_mean-x2_mean)/np.sqrt(x1_var/data1.shape[0]+x2_var/data2.shape[0])
pval=stats.norm.sf(abs(zstat))*2
#当方差不相等时，zstat=(x1_mean-x2_mean)/np.sqrt(x1_var/data2.shape[0]+x2_var/data1.shape[0])


#问题4 分别从男性和女性体温数据中抽取25条数据，验证男性和女性体温均值是否相等
df1_ind=[62,56,44,15,26,37,33,0,43,49,41,60,12,7,52,57,31,50,30,1,39,2,48,32,21]
df2_ind=[129,85,79,82,89,92,117,113,116,96,106,103,114,80,100,76,95,73,101,
          75,118,74,119,81,124]

#默认男性、女性方差相等
#原假设H0：男性体温平均值=女性体温平均值
d1=ws.DescrStatsW(data1.loc[df1_ind,'Temperature'])
d2=ws.DescrStatsW(data2.loc[df2_ind,'Temperature'])
com=ws.CompareMeans(d1,d2)
tstat,pval,_=com.ttest_ind(alternative='two-sided',usevar='unequal',value=0)
#当方差不想等时，usevar='unequal'

x1_mean=np.mean(data1.loc[df1_ind,'Temperature'])
x2_mean=np.mean(data2.loc[df2_ind,'Temperature'])
x1_var=np.var(data1.loc[df1_ind,'Temperature'],ddof=1)
x2_var=np.var(data2.loc[df2_ind,'Temperature'],ddof=1)
s=((len(df1_ind)-1)*x1_var+(len(df2_ind)-1)*x2_var)/(len(df1_ind)+len(df2_ind)-2)
tstat=(x1_mean-x2_mean)/np.sqrt(s*(1/len(df1_ind)+1/len(df2_ind)))
pval=stats.t.sf(abs(tstat),len(df1_ind)+len(df2_ind)-2)*2

"""
若方差不相等时
tstat=(x1_mean-x2_mean)/np.sqrt(x1_var/len(df1_ind)+x2_var/len(df2_ind))
df_1=(x1_var/len(df1_ind)+x2_var/len(df2_ind))**2
df_2=(x1_var/len(df1_ind))**2/(len(df1_ind)-1)+(x2_var/len(df2_ind))**2/(len(df2_ind))
df=df_1/df_2
pval=stats.t.sf(abs(tstat),df)*2
"""
#这个问题的pval=0.08 如果按照0.05的显著性水平，在这里我们无法
#拒绝原假设。这和问题3的结论相反


#其他数据 测得10只电阻的电阻值，验证其方差与60是否有显著差异？(卡方检验,显著水平要求0.1)
val=[114.2,91.9,107.5,89.1,87.2,87.6,95.8,98.4,94.6,85.4] 

#原假设H0:电阻方差=60
val_var=np.var(val,ddof=1)
chi_stat=(len(val)-1)*val_var/60
pval=stats.chi2.sf(np.abs(chi_stat),len(val)-1)*2

这里有以下几点需要说明：

• python包中的检验统计量公式与博客中列出的并不完全一致，差别主要在样本方差、样本标准差以及自由度上。但最终的计算结果与用公式计算出的结果差别不大。
• weightstats中能进行$Z$检验和$t$检验的还有其他方法，这里就不全列出来了。但这些方法中的参数和ztest()方法中的参数大同小异。ztest()方法中的参数及其说用如下。

参数名	作用
x1	接受array-like类型的1维或2维变量。当为2位变量时，默认变量之间相互独立
x2	默认值为空，可接受array-like类型的1维或2维变量。当x2为空时，该函数仅对x1中的变量进行总体均值的假设检验；当变量x2不为空时，做变量x1和变量x2之间的总体均值之差的假设检验
value	当x2为空时，value为原假设下变量x1的均值假设值。当x2不为空时，value为原假设下变量x1和变量x2的均值假设值之差
alternative	备择假设类型 ，取值为“two-sided”时：双边检验，$H_1$中为$\ne$；取值为“larger”时：$H_1$中为$>$；取值为“smaller”时：$H_1$中为$<$
usevar	目前只能取“pooled“”，但是在这个函数中该参数除了抛出异常之外并没有影响实际的计算。但在ztest_ind()函数中除了可以取“pooled"(指示样本方差相等)，还可以取“unequal"(指示样本方差不想等)。
ddof	方差计算时的自由度。

参考资料

1. 中心极限定理
2. 知乎专栏文章: 假设检验
3. 百度文库：假设检验、t检验和Z检验
4. 假设检验的Python实现
5. 假设检验原理及Python实现