概率统计Python计算：总体未知参数的矩估计

设样本$(X_1,X_2,\cdots ,X_n)$来自$X$，$X$的分布中含有$m$个未知参数${\theta }_1,{\theta }_2,\cdots ,{\theta }_m$。设$X$存在直到$m$阶的原点矩$E(X),E(X^2),\cdots ,E(X^m)$。显然，$E(X^k)={\mu }_k$，仍然含有$m$个参数${\theta }_1,{\theta }_2,\cdots ,{\theta }_m$，记为${\mu }_k({\theta }_1,{\theta }_2,\cdots ,{\theta }_m)$，$k=1,2,\cdots ,m$。我们得到方程组
$$\{\begin{array}{l}{\mu }_1={\mu }_1({\theta }_1,{\theta }_2,\cdots ,{\theta }_m)\\ {\mu }_2={\mu }_2({\theta }_1,{\theta }_2,\cdots ,{\theta }_m)\\ \cdots \cdots \cdots \\ {\mu }_m={\mu }_m({\theta }_1,{\theta }_2,\cdots ,{\theta }_m)\end{array}$$
解此方程组，并用样本的$k$阶原点矩$A_k=\frac{1}{n}\sum _{i=1}^n{X}_i^k$代替总体$k$阶原点矩${\mu }_k$，得到参数${\theta }_1,{\theta }_2,\cdots ,{\theta }_m$的矩估计量
$$\{\begin{array}{l}{\stackrel{\wedge }{\theta }}_1=g_1(A_1,A_2,\cdots ,A_m)\\ {\stackrel{\wedge }{\theta }}_2=g_2(A_1,A_2,\cdots ,A_m)\\ \cdots \cdots \cdots \\ {\stackrel{\wedge }{\theta }}_m=g_m(A_1,A_2,\cdots ,A_m)\end{array}.$$
矩估计量的观测值称为矩估计值。
为计算参数的矩估计值，则可用样本数据数组x（numpy.array类对象）的mean函数计算样本1阶原点矩；若需要计算$k$（$\ge 2$）阶原点矩，则只需调用x**k的mean函数即可。
例1随机地取8只活塞环，测得它们的直径为（单位：mm）：
$$74.001,74.005,74.003,74.001,74.000,73.998,74.006,74.002$$
计算总体期望$\mu$，方差${\sigma }^2$的矩估计值，并计算样本方差值$s^2$。
解： 设总体为$X$，由${\sigma }^2=D(X)=E(X^2)-[E(X){]}^2$，得方程组
$$\{\begin{array}{l}\mu =E(X)=A_1\\ {\sigma }^2=E(X^2)-[E(X){]}^2=A_2-A_1^2\end{array}.$$
下列代码计算本例中总体期望$\mu$和方差${\sigma }^2$的矩估计值，以及样本方差观测值$s^2$。

import numpy as np                              #导入numpy
x=np.array([74.001, 74.005, 74.003, 74.001,     #样本观测值
            74.000, 73.998, 74.006, 74.002])
a1=x.mean()                                     #1阶样本原点矩
print('mu=%.4f'%a1)                             #总体均值的矩估计
a2=(x**2).mean()                                #2阶样本原点矩
print('sigma^2=%.5f'%(a2-a1**2))                #总体方差的矩估计
s2=x.var(ddof=1)								#样本方差
print('s^2=%.4f'%s2)

程序的第2~3行设置样本观测值数据，记为x。第4行调用x的mean函数，计算样本均值，即样本1阶原点矩$A_1$，记为a1。第5行显示用$A_1$估计总体期望$\mu$。第6行调用x**2（即$x^2$）的mean函数，计算2阶样本原点矩$A_2$，记为a2。第7行显示用样本原点矩表达式$A_2-A_1^2$估计总体方差${\sigma }^2$。第8~9行调用x的var函数计算并输出样本方差。运行程序，输出

mu=74.0020
sigma^2=0.00001
s^2=0.00001

对于连续型总体，若知道其分布类型，scipy.stats提供的所有分布对象，均有函数fit_loc_scale，其调用接口为
$$\text{fit\_loc\_scale(data)}$$
该函数直接对参数data所传递的样本数据的1、2阶原点矩对该分布的参数loc和scale作矩估计。我们可以通过了解loc和scale与待估参数的关系，算出待估参数的估计值。例如，若总体$X$~$U(a,b)$，scipy.stats表示均匀分布的uniform对象的参数loc对应分布参数$a$，scale参数对应分布参数$a$与$b$的差$b-a$。则loc的估计值就是$a$的矩估计值，scale+loc的估计值就是$b$的估计值。
例2 设总体$X$~$U(a,b)$，$a$和$b$未知。来自$X$的，容量$n=20$的样本观测值为
$$\begin{gathered} 1.248,1.664,1.101,1.967,1.468,1.140,1.434,1.063,1.878,1.375 \\ 1.819,1.704,1.328,1.619,1.830,1.764,1.034,1.553,1.878,1.166 \end{gathered}$$
试计算$a$和$b$的矩估计值。
解： 下列代码完成本利计算。

import numpy as np                              #导入numpy
from scipy.stats import uniform                 #导入uniform 
x=np.array([1.248, 1.664 ,1.101 ,1.967 ,1.468,  #设置样本数据数组
            1.140, 1.434, 1.063, 1.878, 1.375,
            1.819, 1.704, 1.328, 1.619, 1.830,
            1.764, 1.034, 1.553, 1.878, 1.166])
l, s=uniform.fit_loc_scale(x)                   #均匀分布参数loc,scale的矩估计
a=l                                             #a的矩估计
b=a+s                                           #b的矩估计
print('用样本矩估计a=%.4f, b=%.4f'%(a, b))

注意程序中第7行调用uniform对象的fit_loc_scale函数，用数据数组计算表示均匀分布的uniform对象的属性参数loc和scale的矩估计值l和s。根据loc和scale与均匀分布的参数$a$和$b$的关系，第8行和第9行分别计算出$a$和$b$的矩估计值。运行程序，输出

用样本矩估计a=0.9832, b=2.0201

写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！
代码诚可贵，原理价更高。若为AI学，读正版书好。
返回《导引》