概率统计Python计算：中心极限定理的验证

中心极限定理告诉我们，独立同分布的随机变量序列$X_1,X_2,\cdots ,X_n,\cdots$，若$E(X_i)=\mu$，$D(X_i)={\sigma }^2$，$i=1,2,\cdots$。则
$$\underset{n\to \infty }{\lim }P\left(\left|\frac{\frac{{\sum }_{i=1}^n{X}_i}{n}-\mu }{\frac{\sigma }{\sqrt{n}}}\right|\le x\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x{e}^{-\frac{t^2}{2}}dt.$$
即无论诸$X_i$的分布是什么，其前$n$项均值$\overline{X}=\frac{1}{n}\sum _{i=1}^n{X}_i$的标准化变量的绝对值，收敛于标准正态分布。或等价地，$\overline{X}$收敛于$N(\mu ,{\sigma }^2/n)$。
此处，为验证该定理，设$\Theta$~$U\left(-\frac{\pi }{2},\frac{\pi }{2}\right)$，$V=380\sin \Theta$。先计算$V$的分布。
由于$\Theta$的分布函数为
$$F_{\Theta }(\theta )=\{\begin{array}{ll}0& \theta \le -\frac{\pi }{2}\\ \frac{\theta +\frac{\pi }{2}}{\pi }& -\frac{\pi }{2}<\theta <\frac{\pi }{2}\\ 1& \theta \ge \frac{\pi }{2}\end{array}=\{\begin{array}{ll}0& \theta \le -\frac{\pi }{2}\\ \frac{\theta }{\pi }+\frac{1}{2}& -\frac{\pi }{2}<\theta <\frac{\pi }{2}\\ 1& \theta \ge \frac{\pi }{2}\end{array}$$
且函数$g(\theta )=380\sin \theta$，$\theta \in \left(-\frac{\pi }{2},\frac{\pi }{2}\right)$的反函数为$g^{-1}(v)=\arcsin \left(\frac{v}{380}\right)$，$v\in (-380,380)$。于是，$V$的分布函数为
$$F_V(v)=\{\begin{array}{ll}0& v\le -380\\ \frac{1}{\pi }\arcsin \frac{v}{380}+\frac{1}{2}& -380<v<380\\ 1& v\ge 380\end{array}.$$
我们用scipy.stats的rv_continuous类，构造一个分布函数为$F_V(v)$的连续型概率分布vdist（详见博文《连续型随机变量函数分布》）。

from scipy.stats import rv_continuous
import numpy as np
from matplotlib import pyplot as plt
class vdist(rv_continuous):
    def _cdf(self, v):                      #定义累积分布函数
        y=np.zeros(v.size)                  #返回值初始化为0
        d=np.where((v>-380)&(v<380))        #非0非1区域
        y[d]=np.arcsin(v[d]/380)/np.pi+1/2  #非0非1函数值
        d=np.where(v>=380)                  #值为1的区域
        y[d]=1                              #值为1的函数值
dist=vdist()
v=np.linspace(-500, 500, 256)
plt.plot(v, dist.pdf(v), color='black')
plt.show()
Ev=dist.expect()
Dv=dist.var()
print('E(V)=%.4f, D(V)=%.4f'%(Ev, Dv))

程序的第4~11行定义rv_continuous的子类vdist。这只需要做一件事情：在第5~11行定义该分布类的累积分布函数cdf。这是一个属于对象的函数（含有参数self）该函数只含有负责传递外部数据的参数v。该参数是数组类对象，表示分布函数$F_V(v)$的自变量。第6行将函数值初始化为0。根据上述$F_V(v)$的定义式，第7行将介于-380~380的部分区间设为d，第8行将对应区间d的函数值设为$\frac{1}{\pi }\arcsin \frac{v}{380}+\frac{1}{2}$，第9行将不小于380的区间设为d，第10行将对应该区间的函数值置为1。
第12行创建一个服从vdist分布的随机变量对象dist。第13行设置绘图用的横坐标范围v，即区间$(-500,500)$。第14行调用绘图对象plt的plot函数绘制vdist的概率密度函数pdf的图形，第15行显示图形（见下图）。

显然dist分布与正态分布大不相同（其密度函数呈现双峰形态）。
第16、17行分别调用vdist的expect函数和var函数计算其期望$E(V)$与方差$D(V)$。运行程序，输出

E(V)=0.0000, D(V)=72200.0000

为vdist的期望值和方差值。
接下来，利用上列程序创建的dist对象，产生服从该分布的随机变量$V$的10000个值。

V=dist.rvs(size=10000)                                      #产生10000个V的取值
plt.hist(V, density=True, histtype='stepfilled', alpha=0.2) #画直方图
plt.show()

程序的第1行调用dist的rvs函数产生10000个随机数模拟服从dist分布的随机变量$V$，存于数组V。第2行调用绘图对象plt的绘制直方图的函数hist，绘制V的直方图。运行程序，显示如下图形。

可以看到，上图中V数据直方图与前图的密度函数拟合得很好。现在从存储在V中的10000个数据中随机抽取500组，每组100个数据，并计算出每组数据的平均值，然后考察这500个均值的分布。

from scipy.stats import norm                        #导入norm 
x=np.array([np.random.choice(a=V, size=100).mean()  #V中100个随机数据的均值
            for k in range(500)])                   #共生成500个数据存于x
plt.hist(x, density=True, histtype='stepfilled')    #绘制x的直方图
v=np.linspace(-80, 80,256)                          #横坐标区间(-80, 80)
sigma=np.sqrt(Dv)/10                                #计算标准差
plt.plot(v, norm.pdf(v, scale=sigma))               #绘制正态密度曲线
plt.show()

程序的2~3行调用numpy的random对象的choice函数从V中随机取得100个数据（形成一个匿名数组），调用该数组对象的mean函数计算这100个数据的均值，加入数组x中。反复做500次（第3行的for循环决定），故x中存有500个这样的数据。第4行绘制x中数据的直方图。第5行设置绘制正态分布密度函数图形的横坐标范围，区间v=(-80, 80)。第6行计算$V$的均方差$\sqrt{D(V)}=\sigma$。第7行绘制正态分布$N(0,{\sigma }^2/100)$的密度函数norm.pdf图形。运行该程序，展示图形如下图所示。

从上图中我们看到，500个$\frac{1}{100}\left(\sum _{i=1}^{100}{V}_i\right)$数据值作成的直方图与正态分布$N(0,{\sigma }^2/100)$的密度函数图形高度拟合。其中$V_1,V_2,\cdots V_{100}$独立（从V中随机抽取）同分布（与V的分布相同），${\sigma }^2=D(V)$，$0=E(V)$。这就验证了中心极限定理：尽管随机变量$V$所服从的分布与正态分布形态迥异，但是100个独立且与$V$同分布的随机变量$V_1,V_2,\cdots V_{100}$之均值$\frac{V_1+V_2+\cdots +V_{100}}{100}$却近似服从$N(0,{\sigma }^2/100)$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》