用Python对正态分布样本进行区间估计

介绍

在有参估计中，我们有两种常见的估计参数的方法：

• 点估计
• 区间估计

本次估计方法，就是介绍一下用Python对正态分布样本进行区间估计的方法。

数学方法概述

对呈正态分布的样本的参数进行区间估计的方法，我们大致可以分成两类：去估计$\mu$和去估计$\sigma$，在估计$\mu$ 或$\sigma$又分别有2种情况，则一共有以下四种情况：

1. 已知$\sigma$求$\mu$
2. 未知$\sigma$求$\mu$
3. 已知$\mu$求$\sigma$
4. 未知$\mu$求$\sigma$

面对4种不同的情况我们来分别简要讨论一下计算方法。

1. 已知$\sigma$求$\mu$
   由于$\sigma$已知，则我们设$D=\frac{\overline{X}-\mu }{\sigma /\sqrt{n}}\sim N(0,1)$，设某小概率事件为$\alpha$，则置信度设置为$1-\alpha$，为了尽可能缩小横坐标的取值范围，我们选择使用两边分别取$\frac{\alpha }{2}$的方法。
   这样 P { Φ − α 2 < D < Φ α 2 } = 1 − α P \{\Phi_{-\frac {\alpha}{2}} \lt D\lt \Phi_{\frac {\alpha}{2}} \} = 1-\alpha P{Φ−2α​​<D<Φ2α​​}=1−α
   化简求得，置信区间为：
   $$(\overline{X}-\frac{\sigma }{\sqrt{n}}{\Phi }_{-\frac{\alpha }{2}},\overline{X}+\frac{\sigma }{\sqrt{n}}{\Phi }_{\frac{\alpha }{2}})$$
   也就是说在这个区间里面取的值，由百分之$1-\alpha$的可信度。
2. 未知$\sigma$求$\mu$
   由于未知$\sigma$，那么设$D=\frac{\overline{X}-\mu }{S/\sqrt{n}}\sim t(n-1)$因为这里要进行无偏估计，t分布得自由度为$n-1$，其中$S$为样本方差。
   根据上面的方法同理求得，置信区间为：$$(\overline{X}-\frac{\sigma }{\sqrt{n}}{t}_{-\frac{\alpha }{2}},\overline{X}+\frac{\sigma }{\sqrt{n}}{t}_{\frac{\alpha }{2}})$$
3. 已知$\mu$求$\sigma$
   已知$\mu$，则设$D=\frac{{\sum }_{i=1}^n(x_i-\mu {)}^2}{{\sigma }^2}\sim {\chi }^2(n)$
   求解置信区间：
   $$(\frac{{\sum }_{i=1}^n(x_i-\mu {)}^2}{{\chi }_{1-\frac{\alpha }{2}}^2},\frac{{\sum }_{i=1}^n(x_i-\mu {)}^2}{{\chi }_{\frac{\alpha }{2}}^2})$$
4. 未知$\mu$求$\sigma$
   设$D=\frac{(n-1)S^2}{{\sigma }^2}\sim {\chi }^2(n-1)$，同样是为了保证无偏估计。
   求解置信区间为：$$(\frac{(n-1)S^2}{{\chi }_{\frac{\alpha }{2}}^2},\frac{(n-1)S^2}{{\chi }_{1-\frac{\alpha }{2}}^2})$$

估计期望的代码：

import numpy as np
import scipy.stats  as st
from scipy.stats import t
from scipy.stats import norm
import matplotlib.pyplot as plt

#在进行区间估计时，我们首先要知道样本总体的分布情况和一些参数
#x表示数据取值的列表
#sig表示方差
#alpha一个足够小的概率
#confidence表置信度
def mean_ext(x,sig=None,confidence=0.95):
    alpha = 1 - confidence#求alpha
    if sig != None:#当方差已知时
        phi_alpha_2 = norm.isf(-alpha/2)#isf指的是残存函数的逆
        lower_limit = np.mean(x) - phi_alpha_2 * sig / np.sqrt(len(x)) #边界，最小值
        upper_limit = np.mean(x) + phi_alpha_2 * sig / np.sqrt(len(x))#边界，最大值
    else:#当方差未知时
        t_alpha_2 = t.isf(-alpha/2,df = len(x)-1)
        t_test = 1-t.isf(alpha/2,df = len(x)-1)
        std = np.std(x,ddof=1)
        lower_limit = np.mean(x) - t_alpha_2 * std / np.sqrt(len(x)-1)
        upper_limit = np.mean(x) + t_alpha_2 * std / np.sqrt(len(x)-1)
   
   
    return np.array(round(lower_limit, 6), round(upper_limit, 6))#返回一个四舍五入的值
x = np.random.rand(500)生产数据
xnew = np.linspace(0,10,10000) #画图数据
interval_result = mean_ext(x,confidence=0.99) #区间估计结果,为一个区间
u = np.mean(interval_result) #算个平均期望
f = norm(loc = u,scale = np.std(x)).pdf(xnew) #pdf为概率密度函数
f1 = norm(loc = u,scale = np.std(x)).pdf(x)
plt.plot(x,f1,'r+',xnew,f,'g-')
plt.legend(['point','pdf'],loc='best')
plt.show()

画图结果：
估计方差的代码：

import numpy as np
import scipy.stats  as st
from scipy.stats import t
from scipy.stats import norm
import matplotlib.pyplot as plt

#在进行区间估计时，我们首先要知道样本总体的分布情况和一些参数
#x表示参数取值的列表
#mean表期望
#sig表示方差
#alpha上分位
def mean_ext(x,sig=None,confidence=0.95):
    alpha = 1 - confidence
    if sig != None:
        phi_alpha_2 = norm.isf(-alpha/2)
        lower_limit = np.mean(x) - phi_alpha_2 * sig / np.sqrt(len(x))
        upper_limit = np.mean(x) + phi_alpha_2 * sig / np.sqrt(len(x))
    else:
        t_alpha_2 = t.isf(alpha/2,df = len(x)-1)
        t_test = t.isf(1-alpha/2,df = len(x)-1)
        std = np.std(x,ddof=1)
        lower_limit = np.mean(x) - t_alpha_2 * std / np.sqrt(len(x)-1)
        upper_limit = np.mean(x) + t_alpha_2 * std / np.sqrt(len(x)-1)
   
   
    return np.array(round(lower_limit, 3), round(upper_limit, 3))


def sigma_ext(x,mean = None,confidence = 0.99):
    alpha = 1-confidence
    top:float = 0#分子
    if mean != None:
        for xi in x:
            top += (xi-mean)**2
        chi2_alpha_2 = st.chi2.isf(1-alpha/2,df = len(x))
        chi2_alpha_2C = st.chi2.isf(alpha/2,df = len(x))
        chi2_alpha_test = chi2_alpha_2+chi2_alpha_2C
    else:
        d = np.mean(x)
        for xi in x:
            top +=(xi-d)**2
       
        chi2_alpha_2 = st.chi2.isf(1-alpha/2,df = len(x)-1)
        chi2_alpha_2C = st.chi2.isf(alpha/2,df = len(x)-1)
        chi2_alpha_test = chi2_alpha_2+chi2_alpha_2C
    lower_limit:float = top/chi2_alpha_2C
    upper_limit:float = top/chi2_alpha_2

    return np.array(round(lower_limit, 3), round(upper_limit, 3))

#数据准备
x = np.random.normal(0,1,300)
xnew = np.linspace(-3,3,10000)
interval_mean_result1 = mean_ext(x,confidence=0.99)
interval_mean_result2 = mean_ext(x,1,confidence=0.99)
interval_sigma_result1 = sigma_ext(x,confidence=0.99)
interval_sigma_result2 = sigma_ext(x,0,confidence =0.99)
u1 = np.mean(interval_mean_result1)
u2 = np.mean(interval_mean_result2)
sigma1 = np.mean(interval_sigma_result1)
sigma2 = np.mean(interval_sigma_result2)
f = norm(loc = 0,scale = 1).pdf(x)
f1= norm(loc = u1,scale = np.std(x)).pdf(xnew)
f2= norm(loc = u2,scale = 1).pdf(xnew)
f3= norm(loc = np.mean(x),scale = np.sqrt(sigma1)).pdf(xnew)
f4= norm(loc = 0,scale =np.sqrt(sigma2)).pdf(xnew)

#画图模块
plt.figure(figsize=(12,12))
plt.subplot(221)
plt.plot(x,f,'r+',xnew,f1,'g--')
plt.legend(['point','pdf'],loc='best')
plt.title("Mean Estimation(Unknown Sigma)")
plt.subplot(222)
plt.plot(x,f,'r+',xnew,f2,'g--')
plt.legend(['point','pdf'],loc='best')
plt.title("Mean Estimation(Known Sigma)")
plt.subplot(223)
plt.plot(x,f,'r+',xnew,f3,'b--')
plt.legend(['point','pdf'],loc='best')
plt.title("Sigma Estimation(Unknown Mean)")
plt.subplot(224)
plt.plot(x,f,'r+',xnew,f4,'b--')
plt.legend(['point','pdf'],loc='best')
plt.title("Sigma Estimation(Known Mean)")
plt.show()

图像：