16 极大似然估计

16.1 自定义极大似然函数

16.1.1 正态分布

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

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

$$\begin{array}{cc}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{array}$$

对数最大似然函数：

$$\begin{array}{cc}lnL(\mu ,\sigma )& =nln\frac{1}{\sqrt{2\pi }\sigma }-\sum _{i=1}^n\frac{(x_i-\mu {)}^2}{2{\sigma }^2}\end{array}$$

得到对数最大似然函数后，再利用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\left(x\right)=\{\begin{array}{cc}\lambda e^{-\lambda x},& x>0\\ 0,& x\le 0\end{array}$$

最大似然函数：

$$\begin{array}{cc}L\left(\lambda \right)& =\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{array}$$

对数最大似然函数：

$$\begin{array}{cc}lnL\left(\lambda \right)& =nln\lambda -\lambda \sum _{i=1}^n{x}_i\end{array}$$

得到对数最大似然函数后，再利用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。