16 极大似然估计

16.1 自定义极大似然函数

16.1.1 正态分布

正态分布的概率密度函数：

\[ f(x)=\frac{1}{\sqrt{2 \pi}\sigma} e ^{- \frac{(x-\mu)^{2}}{2\sigma^2}} \] 最大似然函数：

\[ \begin{aligned} L(\mu,\sigma)&=\prod_{i=1}^{n} f(x_i,\mu,\sigma) \\ &=\prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi}\sigma} e ^{- \frac{(x-\mu)^{2}}{2\sigma^2}} \\ & = (\frac{1}{\sqrt{2 \pi} \sigma})^n - \sum_{i=1}^n \frac{(x_i-\mu)^2}{2\sigma^2} \end{aligned} \]

对数最大似然函数：

\[ \begin{aligned} lnL(\mu,\sigma)&=n ln \frac{1}{\sqrt{2 \pi} \sigma} - \sum_{i=1}^{n} \frac{(x_i-\mu)^2}{2\sigma^2} \end{aligned} \]

得到对数最大似然函数后，再利用R中的optim函数求解对数最大似然函数取极值时对应的参数值 \(\mu\)和\(\sigma\)。

norm.fun<- function(theta,x) # 创建似然函数
{  
  mu<- theta[1]
  sigma<- theta[2]   
  n<- length(x)
  # 似然函数，再求其对数
  logL<- n*log( 1/( sqrt(2*pi)*sigma))  - sum((x-mu)^2)/(2*sigma^2)
  return(-logL)  # 因为R中只有最小化函数optim(),我们只需要参数的值，最大化 logL 和最小化 -logL 一致的
}  

theta0=c(0,1) #初始化两个参数作为迭代初始值
X=1:50 #抽出的样本（这里我以1至50作为数据样本）

result=optim(theta0,norm.fun,x=X)
# 用R中内置的optim()最小化函数对 似然函数优化，求出最优值的两个参数u、sigma
# 参数的值保存在result$'par' 中，有两个值，第一个是u 第二个是sigma
result

## $par
## [1] 25.48969 14.43959
## 
## $value
## [1] 204.4154
## 
## $counts
## function gradient 
##       81       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

在\(lnL\)取得最大值时（即\(-lnL\)取得最小值时），解出的\(\mu\)的参数值为25.48969，\(\sigma\)的参数估计值为14.43959，\(-lnL\)为204.4154。

16.1.2 指数分布

指数分布的概率密度函数：

\[ f(x)= \begin{cases} \lambda e^{-\lambda x}, & x>0\\ 0, & x\le0 \end{cases} \]

最大似然函数：

\[ \begin{aligned} L(\lambda)&=\prod_{i=1}^{n} f(x_i,\lambda) \\ &=\prod_{i=1}^{n}\lambda e^{-\lambda x} \\ & = \lambda^n - e^{-\lambda \sum_{i=1}^n x_i} \end{aligned} \]

对数最大似然函数：

\[ \begin{aligned} lnL(\lambda)&=nln\lambda -\lambda \sum_{i=1}^n x_i \end{aligned} \]

得到对数最大似然函数后，再利用R中的optim函数求解对数最大似然函数取极值时对应的参数值 \(\lambda\)。

exp.fun <- function(theta,x){
  n <- length(x)
  logL <- n*log(theta)-theta*sum(x)
  return(-logL)
}
x=1:50
lamda=1e-2
result=optim(par = lamda,fn = exp.fun,x=X) #运行会有一些警告，说optim（）优化的结果不一定正确，但这里是正确的，不影响

## Warning in optim(par = lamda, fn = exp.fun, x = X): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly

result 

## $par
## [1] 0.03921875
## 
## $value
## [1] 211.9339
## 
## $counts
## function gradient 
##       30       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

在\(lnL\)取得最大值时（即\(-lnL\)取得最小值时），解出的\(\lambda\)的参数值为 0.03921875，\(-lnL\)为211.9339。