一元正态总体参数的极大似然估计_附Python实现

1.极大似然估计思想

    我们获得了样本$x_1,x_2,...,x_m$，这些样本被我们从总体中抽到了，从概率的角度思考，我们认为这些样本出现的概率是极大的，因此，可以构造出这些样本的联合概率函数，使得联合概率函数最大的参数值就是极大似然估计值.
    现假设随机变量$X_1,X_2,...,X_m$的联合概率密度函数或联合概率质量函数为$p(X;\theta \theta )$，$\theta \theta \in \Theta \Theta$，其中$\theta \theta$是一个或多个未知参数组成的参数向量，$\Theta \Theta$是参数空间，$x_1,x_2,...,x_m$是来自总体的样本，将样本的联合概率函数看成是$\theta \theta$的函数，用$L(\theta \theta ;x_1,x_2,...,x_m)$表示，简记为$L(\theta \theta )$.
    如果$X_i$是i.i.d的，那么$L(\theta \theta )$可以写为
$$L(\theta \theta )=L(\theta \theta ;x_1,x_2,...,x_m)=\prod _{i=1}^{m}p(x_i,\theta \theta )$$
$L(\theta \theta )$称为样本的似然函数。如果统计量
$$\widehat{\theta \theta }=\hat{\theta }\hat{\theta }(x_1,x_1,...,x_m)$$
满足
$$L(\widehat{\theta \theta })=\underset{\theta \theta \in \Theta \Theta }{\max }L(\theta \theta )$$
则称$\widehat{\theta \theta }$是$\theta \theta$的极大似然估计，简记为MLE(maximum likelihood estimate).

2.一元正态总体参数MLE推导

    设样本$x_1,x_2,...,x_m$来自正态总体$N(\mu ,{\sigma }^2)$，且相互独立，则$\theta \theta =(\mu ,{\sigma }^2)$.由于样本来自正态总体，则
p ( x i ; μ , σ 2 ) = 1 2 π σ e x p { − ( x i − μ ) 2 2 σ 2 } , i = 1 , 2 , . . . , m p(x_{i};\mu,\sigma^{2})=\frac{1}{\sqrt{2\pi}\sigma}exp\{-\tfrac{(x_{i}-\mu )^{2}}{2\sigma^{2}}\} ,i=1,2,...,m p(xi​;μ,σ2)=2π ​σ1​exp{−2σ2(xi​−μ)2​},i=1,2,...,m
故样本的似然函数为
L ( θ ) = ∏ i = 1 m p ( x i , θ ) = ∏ i = 1 m 1 2 π σ e x p { − ( x i − μ ) 2 2 σ 2 } = ( 2 π σ 2 ) − m / 2 e x p { − 1 2 σ 2 ∑ i = 1 m ( x i − μ ) 2 } \begin{aligned} L(\pmb{\theta})&=\prod_{i=1}^{m} p(x_{i},\pmb{\theta })\\ &=\prod_{i=1}^{m} \frac{1}{\sqrt{2\pi}\sigma}exp\{-\tfrac{(x_{i}-\mu )^{2}}{2\sigma^{2}}\} \\ &=(2\pi\sigma^{2})^{-m/2}exp\{-\frac{1}{2\sigma^{2}}\sum_{i=1}^{m}(x_{i}-\mu )^{2}\} \end{aligned} L(θθ)​=i=1∏m​p(xi​,θθ)=i=1∏m​2π ​σ1​exp{−2σ2(xi​−μ)2​}=(2πσ2)−m/2exp{−2σ21​i=1∑m​(xi​−μ)2}​
对数似然函数为
$$\begin{array}{rl}lnL(\theta \theta )& =-\frac{m}{2}ln(2\pi )-\frac{m}{2}ln({\sigma }^2)-\frac{1}{2{\sigma }^2}\sum _{i=1}^m(x_i-\mu {)}^2\end{array}$$
将$lnL(\theta \theta )$对两个分量分别求偏导得
$$\frac{\partial lnL(\theta \theta )}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^m(x_i-\mu )=0$$
$$\frac{\partial lnL(\theta \theta )}{\partial {\sigma }^2}=\frac{1}{2{\sigma }^4}\sum _{i=1}^m(x_i-\mu {)}^2-\frac{m}{2{\sigma }^2}=0$$
解上述方程组得$\mu$的极大似然估计为
$$\hat{\mu }=\frac{1}{m}\sum _{i=1}^m{x}_i=\bar{x}$$
${\sigma }^2$的极大似然估计为
$$\widehat{{\sigma }^2}=\frac{1}{m}\sum _{i=1}^m(x_i-\bar{x}{)}^2$$

3.代码实现

代码使用python进行实现，本文使用两种方式实现参数估计

• 根据估计量的表达式进行计算.
• 使用scipy的optimize模块进行优化.

3.1.抽取样本

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optim

首先从$\mu =2,\sigma =4$的正态分布中抽取$m=100000$个样本

np.random.seed(10)
mu=2;sigma=4;m=100000
X=np.sort(np.random.normal(mu,sigma,m))
y=(1/np.sqrt(2*np.pi)*sigma)*np.exp(-(X-mu)**2/(2*sigma**2))

使用样本绘制图形

def plot_density(X,y,mu,sigma):
    fig=plt.Figure(figsize=(8,6))
    grid=plt.GridSpec(1,1)
    axes=fig.add_subplot(grid[0,0])
    axes.plot(X,y)
    axes.set_xlabel("x")
    axes.set_ylabel("$p(x)$")
    for dire in ['top','right']:
        axes.spines[dire].set_visible(False)
    axes.vlines(2,ymin=0,ymax=np.max(y),color='red')
    axes.text(8, 0.25, '$\mu={},\sigma={}$'.format(mu,sigma))
    fig.subplots_adjust(left=0.1, bottom=0.12, right=0.96, top=0.95, wspace=0.5, hspace=0.5)
    axes.text(mu-0.5, 0.5*np.max(y), '$x={}$'.format(mu))
    return fig
fig1=plot_density(X,y,mu,sigma)

3.2.根据估计量的表达式进行计算

$\mu$的极大似然估计为
$$\hat{\mu }=\frac{1}{m}\sum _{i=1}^m{x}_i=\bar{x}$$
${\sigma }^2$的极大似然估计为
$$\widehat{{\sigma }^2}=\frac{1}{m}\sum _{i=1}^m(x_i-\bar{x}{)}^2$$
根据上述两式计算：

hat_mu=np.mean(X)
hat_sigma=np.sqrt((X-hat_mu).dot(X-hat_mu)/m)
print("The estimated value of mean:",hat_mu)
print("The estimated value of Standard deviation:",hat_sigma)

结果：
The estimated value of mean: 1.9880863402374647
The estimated value of Standard deviation: 3.9976758117844575

3.3.使用scipy的optimize模块进行优化

极大似然估计要求使得似然函数最大化的参数$\theta \theta$，也即
$$\hat{\theta }\hat{\theta }=agr\underset{\theta \theta \in \Theta \Theta }{\max }L(\theta \theta )$$
上述优化问题等价于
$$\hat{\theta }\hat{\theta }=agr\underset{\theta \theta \in \Theta \Theta }{\min }-L(\theta \theta )$$
即
θ ^ = a g r min ⁡ θ ∈ Θ − L ( θ ) = a g r min ⁡ θ ∈ Θ − { − m l n ( σ ) − m 2 l n ( 2 π ) − 1 2 σ 2 ∑ i = 1 m ( x i − μ ) 2 } = a g r min ⁡ θ ∈ Θ { m l n ( σ ) + m 2 l n ( 2 π ) + 1 2 σ 2 ∑ i = 1 m ( x i − μ ) 2 } \begin{aligned} \pmb{\hat{\theta}} &=agr\min_{\pmb{\theta}\in \pmb{\Theta}}-L(\pmb{\theta})\\ &=agr\min_{\pmb{\theta}\in \pmb{\Theta}}-\{-mln(\sigma)-\frac{m}{2}ln(2\pi)-\frac{1}{2\sigma^{2}}\sum_{i=1}^{m} (x_{i}-\mu)^{2}\}\\ &=agr\min_{\pmb{\theta}\in \pmb{\Theta}}\{mln(\sigma)+\frac{m}{2}ln(2\pi)+\frac{1}{2\sigma^{2}}\sum_{i=1}^{m} (x_{i}-\mu)^{2}\} \end{aligned} θ^θ^​=agrθθ∈ΘΘmin​−L(θθ)=agrθθ∈ΘΘmin​−{−mln(σ)−2m​ln(2π)−2σ21​i=1∑m​(xi​−μ)2}=agrθθ∈ΘΘmin​{mln(σ)+2m​ln(2π)+2σ21​i=1∑m​(xi​−μ)2}​
由于右边第二项$\frac{m}{2}ln(2\pi )$是常数，不影响优化，故可以去掉。
将上述极小化问题用代码实现如下：

def likelihood_func(theta,X):
	#theta表示要优化的参数向量，其中theta[0]表示均值，theta[1]表示标准差。
    m=len(X)
    z=m*np.log(theta[1])+(1/(2*theta[1]**2))*((X-theta[0]).dot(X-theta[0]))
    return z

将${\mu }_0=4,{\sigma }_0=8$作为初始点，执行优化如下：

pt=np.array([4,8])
result = optim.minimize(likelihood_func, pt,args=(X,), method='Nelder-Mead')
print('Status : %s' % result['message'])
print('Total Evaluations: %d' % result['nfev'])
solution = result['x']
print('Solution: f(%s) ' % solution)

结果：
Status : Optimization terminated successfully.
Total Evaluations: 82
Solution: f([1.98805021 3.99765954])

4.总结

    本文根据极大似然估计的思想给出了极大似然估计一般定义，并推导了一元正态总体参数的极大似然估计，同时，使用两种计算方式得出了参数的极大似然估计。从两种计算结果来看，估计差异很小，且估计误差也很小。两种方式中，$\mu$和$\sigma$的估计值均非常接近真实值。

参考文献

[1] 茆诗松,程依明,濮晓龙.概率论与数理统计教程[M]. 北京:高等教育出版社, 2019. 277-279.
[2] Rice J A. Mathematical statistics and data analysis[M]. Cengage Learning, 2006.

所有代码

#导入包
import scipy.optimize as optim
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt

#抽取样本
np.random.seed(10)
mu=2;sigma=4;m=100000
X=np.sort(np.random.normal(mu,sigma,m))
y=(1/np.sqrt(2*np.pi)*sigma)*np.exp(-(X-mu)**2/(2*sigma**2))

#绘制图形
def plot_density(X,y,mu,sigma):
    fig=plt.Figure(figsize=(6,4))
    grid=plt.GridSpec(1,1)
    axes=fig.add_subplot(grid[0,0])
    axes.plot(X,y)
    axes.set_xlabel("x")
    axes.set_ylabel("$p(x)$")
    for dire in ['top','right']:
        axes.spines[dire].set_visible(False)
    axes.vlines(2,ymin=0,ymax=np.max(y),color='red')
    axes.text(8, 0.25, '$\mu={},\sigma={}$'.format(mu,sigma))
    fig.subplots_adjust(left=0.1, bottom=0.12, right=0.96, top=0.95, wspace=0.5, hspace=0.5)
    axes.text(mu-0.5, 0.5*np.max(y), '$x={}$'.format(mu))
    return fig
fig1=plot_density(X,y,mu,sigma)

#根据估计量的表达式进行计算
hat_mu=np.mean(X)
hat_sigma=np.sqrt((X-hat_mu).dot(X-hat_mu)/m)
print("The estimated value of mean:",hat_mu)
print("The estimated value of Standard deviation:",hat_sigma)

#使用scipy的optimize模块进行优化
def likelihood_func(theta,X):
    m=len(X)
    z=m*np.log(theta[1])+(1/(2*theta[1]**2))*((X-theta[0]).dot(X-theta[0]))
    return z
pt=np.array([4,8])
result = optim.minimize(likelihood_func, pt,args=(X,), method='Nelder-Mead')
print('Status : %s' % result['message'])
print('Total Evaluations: %d' % result['nfev'])
solution = result['x']
print('Solution: f(%s) ' % solution)