Python统计学04——参数估计

参考书目：贾俊平. 统计学——Python实现. 北京: 高等教育出版社，2021.

---

（最近很多同学找我要这个数据，可以参考：贾俊平统计学，整本书的实验数据都在里面）

---

推断统计中最重要的核心手段就是参数估计和假设检验。本章介绍用Python实现参数估计，主要是代码实现，依赖scipy包，原理比较基础简单也不多介绍了。

先导入包

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

plt.rcParams ['font.sans-serif'] ='SimHei'               #显示中文
plt.rcParams ['axes.unicode_minus']=False               #显示负号

首先模拟一下大数中心极限定理。统计量的三个性质

无偏性

x,m,v = [],[],[]
n=10
for i in range(10000):
    d= np.random.normal(loc=50,scale=10,size=n)
    x.append(np.mean(d))
    m.append(np.median(d))
    v.append(np.var(d,ddof=1))
print(f"样本均值的均值：{np.mean(x): .4f}\n样本中位数的均值：{np.mean(m): .4f}\n 样本方差的均值：{np.mean(v):.4f}")

 有效性

x,m=[],[]
n=10
for i in range(10000):
    d = np.random.normal(size=n)
    x.append(np.mean(d))
    m.append(np.median(d))
print("样本的均值的方差:%.4f\n样本中位数的方差: %.4f" %(np.var(x,ddof=1), np.var(m,ddof=1)))

画图看

#绘制样本均值和样本中位数分布的直方图
def plot_dis(df,ax,xlabel):
    df.plot(bins=20,kind="hist",density=True,ax=ax,legend=False)
    df.plot(kind="density",linewidth=2,ax=ax,legend=False)
    ax.set_ylim(0,1.3)
    ax.set_xlim(-1.5,1.5)
    ax.set_xlabel(xlabel)
    ax.set_title(xlabel+"的分布")
    
plt.subplots(1,2,figsize=(12,5))
ax1=plt.subplot(121)
plot_dis(pd.DataFrame(x),ax1,"样本均值")
ax2=plt.subplot(122)
plot_dis(pd.DataFrame(m),ax2,"样本中位数")
plt.show()

 

 一致性

np.random.seed(2020)
N=np.random.normal(loc=50,scale=10,size=1000)
mu = np.mean(N)

xbar10=np.mean(np.random.choice(N,10,replace=False))
xbar100 = np.mean(np.random.choice(N,100,replace=False))
xbar500 = np.mean(np.random.choice(N,500,replace=False))
xbar900 = np.mean(np.random.choice(N,900,replace=False))

pd.DataFrame([mu,xbar10,xbar100,xbar500,xbar900],["总体均值","xbar10","xbar100","xbar500","xbar900"]).T

---

总体均值区间的估计 

读取案例数据，然后计算这条数据的90%置信区间

import scipy.stats as st
example5_1 = pd.read_csv("example5_1.csv",encoding="gbk")
int=st.norm.interval(0.90,loc=np.mean(example5_1),scale=st.sem(example5_1))
                #样本均值=mmean(example5_1) 均值的标准误差st.sem(example5_1)
np.round(int,4) #以数组形式返回置信区间

 使用t分布的情况，要传入自由度

example5_2 = pd.read_csv("example5_2.csv",encoding="gbk")
int=st.t.interval(0.95,len(example5_2)-1,loc=np.mean(example5_2),scale=st.sem(example5_2))
                #样本均值=mmean(example5_2) 均值的标准误差st.sem(example5_2)
np.round(int,4) #以数组形式返回置信区间

---

两个总体均值差的估计

读取案例数据

from scipy.stats  import norm
example5_3 = pd.read_csv("example5_3.csv",encoding="gbk")
example5_3 

 计算两个总体均值之差的置信区间

x1=example5_3 ["男性工资"];x2=example5_3 ["女性工资"]

conf_level=0.95
xbar1=x1.mean();s1=x1.std();n1=len(x1)
xbar2=x2.mean();s2=x2.std();n2=len(x2)

interval=norm.interval(alpha=conf_level,loc=(xbar1-xbar2),scale=np.sqrt(s1**2/n1+s2*2/n2))
print(f"男女平均工资之差的95%置信区间：{np.round(interval,2)}")

 当然两个总体还分为总体方差相等和不相等的情况，自由度相同不同的情况，具体的联合方差或者联合自由度会有区别。

---

总体方差区间估计

方差检验不用z或者t分布了，使用卡方分布

读取案例数据

#总体方差区间估计
from scipy.stats  import chi2
example5_2 = pd.read_csv("example5_2.csv",encoding="gbk")
example5_2

x=example5_2["食品重量"]

conf_level = 0.95
sigma2 = x.var();n=len(x)

LCI=(n-1)*sigma2/chi2.ppf(q=(1+conf_level)/2,df=n-1)
UCI=(n-1)*sigma2/chi2.ppf(q=(1-conf_level)/2,df=n-1)
print(f"置信区间95为:[{LCI:.5f}   {UCI:.5f}]")

---

两个总体方差比的估计 

方差比用F分布

from scipy.stats  import f
example5_4 = pd.read_csv("example5_4.csv",encoding="gbk")
example5_4

x1=example5_4 ["方法一"];x2=example5_4 ["方法二"]

conf_level=0.95
var1=x1.var();var2=x2.var()
n1=len(x1);n2=len(x2)

LCI=(var1/var2)/f.ppf(q=(1+conf_level)/2,dfn=n1-1,dfd=n2-1)
UCI=(var1/var2)/f.ppf(q=(1-conf_level)/2,dfn=n1-1,dfd=n2-1)
print(f"置信区间95为:[{LCI:.5f}  {UCI:.5f}]")