参考书目:贾俊平. 统计学——Python实现. 北京: 高等教育出版社,2021.
参数估计和假设检验是统计学的核心。上一节介绍了常见的参数估计,本章介绍Python实现假设检验的流程。
导入包,读取案例数据
- import numpy as np
- import pandas as pd
- from statsmodels.stats.weightstats import ztest
-
- example6_3=pd.read_csv("example6_3.csv",encoding="gbk")
- example6_3.head()
数据前五行
对该组数据检验,原假设均值是否大于81,备择假设是否小于81。
计算如下
- z, p_value=ztest(x1=example6_3["PM2.5值"],value=81,alternative="smaller")
- print(f"样本均值={example6_3['PM2.5值'].mean() :.2f}, z统计量值={z: .4f}, p值={p_value: .4f}")
主要看p值,p大于0.05,说明不能拒绝原假设,现在的pm2.5的均值没有显著小于81,空气质量不一定变好了。
读取案例数据
- example6_4=pd.read_csv("example6_4.csv",encoding="gbk")
- example6_4.head()
对该组数据检验,原假设均值等于5,备择假设是不等于5。
计算如下
- from scipy.stats import ttest_1samp
- t,p_value = ttest_1samp(a=example6_4['厚度'],popmean=5) #假设popmean为总体均值
- print(f"样本均值={example6_4['厚度'].mean():.2f},t统计量值={t:.4f},p值={p_value:.4g}")
p值很小,说明拒绝原假设,均值显著不为5。
既然整体均值不为5,那么每个样本平均下来均值和5差了多少呢,计算如下
- mu = 5
- xbar=example6_4["厚度"].mean()
- sd=example6_4["厚度"].std()
- d=abs(xbar-mu)/sd
- print(f"效应量 d={d:.4f}")
平均均值和标准的均值差了1.2583个标准差
读取案例数据
- from statsmodels.stats.weightstats import ztest
- example6_5=pd.read_csv("example6_5.csv",encoding="gbk")
- example6_5.head()
计算
- z,p_value = ztest(x1=example6_5['男生上网时间'],x2=example6_5['女生上网时间'],alternative='two-sided')
- print(f"男生平均上网时间为{example6_5['男生上网时间'].mean():.4f}, 女生平均上网时间为{example6_5['女生上网时间'].mean():.4f}\
- \n z统计量值={z:.4f},p值={p_value:.4f}")
p大于0.05,没有证据说明男女生上网时间有显著性差异。
当然两个总体的检验还分为总体方差是否相等,样本自由度是否相等很多种不同的情况,计算和上面这种情况有差异的。但是公式过于繁琐,而且用的少就不展示了。
- #一个总体方差的检验
- from scipy.stats import chi2
- example6_11=pd.read_csv("example6_11.csv",encoding="gbk")
- example6_11.head()
计算
- x=example6_11["填装量"]
- sigma0_2=16
- s2=x.var();n=len(x);df=n-1
- chi2_value=(n-1)*s2/sigma0_2
- p_value=1-chi2.cdf(chi2_value,df=df)
- print(f"样本填装量的方差={s2:.4f},卡方统计量值={chi2_value:.4f},自由度df={df},p值={p_value:.4f}")
- from scipy.stats import f
- example6_6=pd.read_csv("example6_6.csv",encoding="gbk")
- example6_6.head()
计算
- x1=example6_6["甲企业"];x2=example6_6["乙企业"]
- s1_square=x1.var();dfn=len(x1)-1
- s2_square=x2.var();dfd=len(x2)-1
- f_value=s1_square/s2_square
- p_value=f.cdf(f_value,dfn=dfn,dfd=dfd)*2
- print(f"F统计量值={f_value:.4f}, p值={p_value:.4f}")
有时候我们要求数据是服从正态分布的是最好的,因此要进行检验,看数据的分布是否为正态分布,直观上可以画QQ图看。还可以做K-S检验。
- #正态概率图
- import pandas as pd
- import seaborn as sns
- import statsmodels.api as sm
- from matplotlib import pyplot as plt
- plt.rcParams["font.sans-serif"]=["SimHei"]
- plt.rcParams["axes.unicode_minus"]=False
-
- example6_3=pd.read_csv("example6_3.csv",encoding="gbk")
- example6_3.head()
画图
- pplot=sm.ProbPlot(example6_3["PM2.5值"],fit=True)
- fig=plt.figure(figsize=(12,5))
-
- #绘制Q-Q图
- ax1 = plt.subplot(141)
- pplot.qqplot(line="r",ax=ax1,xlabel="期望正态值",ylabel="标准化的观测值")
- ax1.set_title("正态性Q-Q图",fontsize=15)
-
- #绘制P-P图
- ax2 = plt.subplot(142)
- pplot.ppplot(line="45",ax=ax2,xlabel="期望的累积概率",ylabel="观测的累积概率")
- ax2.set_title("正态性p-p图",fontsize=15)
-
- #绘制核密度图
- ax3 = plt.subplot(143)
- sns.kdeplot(example6_3["PM2.5值"],shade=True,alpha=0.5,lw=0.6)
- #plt.xlabel("密度")
- #plt.ylabel("PM2.5")
- ax3.set_title("核密度图",fontsize=15)
-
- #绘制箱线图
- ax4 = plt.subplot(144)
- sns.boxplot(data=example6_3["PM2.5值"],width=0.6,fliersize=2)
- fig.tight_layout()
K-S检验
- #KS检验
- from scipy.stats import kstest
- d=example6_4['厚度']
- D,p_value=kstest(d,'norm',alternative='two-sided',mode='asymp',args=(d.mean(),d.std()))
- print(f"统计量D={D:.5f}, p值={p_value:.4g}")
p>0.05,不拒绝原假设,没有证据证明不服从正态分布。