python实现抽样分布单个正态总体模拟

数理统计中，我们利用样本数据推断总体的信息。在应用时，往往不是直接使用样本本身，而是针对不同的问题构造样本统计量，利用样本统计量进行推断。因此，清楚样本统计量服从什么分布就至关重要了。样本统计量的分布就称为抽样分布。这里，样本是随机变量，样本统计量作为样本的函数，自然也是随机变量，因此才有分布。它是指当多次抽样，得到多个样本，每个样本都可以计算一个样本统计量的观测值；由于抽样的随机性，这些值一般不相等，但有一定的规律，即服从某种分布。抽样分布就指出了它们所服从的分布。

以单个正态总体的均值为例，设总体?~?(?,${\sigma }^2$), 典型的参数估计问题是，已知一个样本的观测值为$x_1,x_2$,⋯,$x_n$, 求总体均值?的估计。我们知道样本均值? ̅ = 1/? ∑ ?? =1 是总体均值?的无偏估计，给定一个样本的观测值?1,?2,⋯,??, 可以 计算出样本均值的观测值$\bar{X}=\frac{1}{n}{\sum }_{i=0}^n{x}_i=1$ . 此处容易迷惑了，由一个样本得到了一个样本均值的观测值，怎么还有样本均值的分布呢？事实上，这里所指的分布，是指多次抽样，得到多个样本，每个样本都可以计算样本均值。由于抽样的随机性，这些样本均值一般来说不相等，但它们的取值有一定的规律，即服从一定的分布。比如，方差${\sigma }^2$已知时，$\bar{X}\sim (\mu ,\frac{{\sigma }^2}{n})$。
下面的定理给出了单个总体的抽样分布结论：
定理 设$x_1,x_2$,⋯,$x_n$是来自正态总体(?,${\sigma }^2$)的样本，$\bar{X},S^2$分别为样本均值和样本方差，则有：

1. $\bar{X}$~?(?,$\frac{{\sigma }^2}{n}$)，进而$\frac{\bar{X}-\mu }{\sigma /\sqrt{n}}\sim N(0,1)$;
2. $(n-1)\frac{S^2}{{\sigma }^2}\sim {\chi }^2(n-1)$
3. ${\bar{X}}^2$与$S^2$相互独立；
4. $\frac{\bar{X}-\mu }{S/\sqrt{n}}\sim t(n-1)$
   下面，利用 numpy 生成服从正态分布的随机数，每一列表示一个样本，利用得到的多个样本，模拟抽样，并比较由抽样结果得到的概率密度估计值和上面定理给出的理论抽样分布的概率密度曲线。
   （1）方差 sigma2 已知时,样本均值的抽样分布 即$\frac{\bar{X}-\mu }{\sigma /\sqrt{n}}\sim N(0,1)$

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
from numpy import matlib
import math
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
n = 10  # 样本容量
m = 10000  # 抽样次数
# N(0,1)
mu = 0.0
sigma = 1.0
X = sigma * matlib.randn((n, m)) + mu  # 模拟抽样，每一列是一个样本容量为n的样本，共抽样m次
Xmean = np.mean(X, axis=0)  # 样本均值
z = (Xmean - mu) / (sigma / math.sqrt(n))  # （? ̅−? ）/（?/√?）~?(0,1);
Z = np.array(z).flatten()
x = np.linspace(-5, 5, 100)
N, bins, patches = plt.hist(Z, 50, density=1, facecolor="yellow", edgecolor="black", alpha=0.7)
plt.plot(x, stats.norm.pdf(x), "b", label="N(0,1)")
plt.xlabel('x')
plt.ylabel('pdf')
plt.suptitle("方差sigma2已知时,样本均值的抽样分布")
plt.legend()
plt.show()

(2)方差sigma2未知时,样本均值的抽样分布 即$\frac{\bar{X}-\mu }{S/\sqrt{n}}\sim t(n-1)$

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
from numpy import matlib
import math
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
n = 10  # 样本容量
m = 10000  # 抽样次数
#N(0,1)
mu=0.0
sigma=1.0
X=sigma * matlib.randn((n, m)) + mu#模拟抽样，每一列是一个样本容量为n的样本，共抽样m次
Xmean=np.mean(X,axis=0)#样本均值
Xstd=np.std(X,axis=0)#标准差
T=(Xmean-mu)/(Xstd/math.sqrt(n))#（? ̅−?）/（?/√?）~?(?−1).
t=np.array(T).flatten()
x=np.linspace(-5, 5, 100)
N, bins, patches =plt.hist(t, 50, density=1, facecolor="yellow", edgecolor="black", alpha=0.7)
plt.plot(x, stats.norm.pdf(x),"b", label="N(0,1)")
plt.plot(x, stats.t.pdf(x, n-1),"r", label="t(n-1)")
plt.xlabel('x')
plt.ylabel('pdf')
plt.suptitle("方差sigma2未知时,样本均值的抽样分布")
plt.legend()
plt.show()

(3)方差 S^2 的抽样分布 即 $(n-1)\frac{S^2}{{\sigma }^2}\sim {\chi }^2(n-1)$

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
from numpy import matlib
import math
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
n = 10  # 样本容量
m = 10000  # 抽样次数
#N(0,1)
mu=0.0
sigma=1.0
X=sigma * matlib.randn((n, m)) + mu#模拟抽样，每一列是一个样本容量为n的样本，共抽样m次
Xmean=np.mean(X,axis=0)#样本均值
zhen2=np.var(X, axis=0)#样本方差
xchi2=(np.dot((n), np.array(zhen2)))/(sigma*sigma)#(?−1)?2/?2 ~?2(? −1)
Xchi2=np.array(xchi2).flatten()
c2=stats.chi2.isf(0.0000001, n-1)
x=np.linspace(0, c2, 100)
N, bins, patches =plt.hist(Xchi2, 50, density=1, facecolor="yellow", edgecolor="black", alpha=0.7)
plt.plot(x, stats.chi2.pdf(x, n-1),"b", label="$\chi2(n-1)$")
plt.xlabel('x')
plt.ylabel('pdf')
plt.suptitle("方差 S^2 的抽样分布")
plt.legend()
plt.show()