补充资料——自己实现极大似然估计(最大似然估计)MLE

这篇文章给了我一个启发，我们可以自己用已知分布的密度函数进行组合，然后构建一个新的密度函数啦，然后用极大似然估计MLE进行估计。

代码和结果演示

代码：

#取出MASS包这中的数据
data(geyser,package="MASS")
head(geyser)
attach(geyser)
par(bg='lemonchiffon')
hist(waiting,freq=F,col="lightcoral")
#freq=F要加上，否则就无法添加线了
lines(density(waiting),lwd=2,col="cadetblue4")
#根据图像，我们认为其在前后分别是两个正态分布函数的组合
#定义 log‐likelihood 函数
LL<-function(params,data){
  #参数"params"是一个向量，
  #依次包含了五个参数： p,mu1,sigma1,mu2,sigma2.
  #参数"data"，是观测数据。
  t1<-dnorm(data,params[2],params[3])
  t2<-dnorm(data,params[4],params[5])
  #f是概率密度函数
  f<-params[1]*t1+(1-params[1])*t2
  #混合密度函数
  ll<-sum(log(f))
  #log‐likelihood 函数
  return(-ll)
  #nlminb()函数是最小化一个函数的值，
  #但我们是要最大化 log‐likeilhood 函数
  #所以需要在“ ll”前加个“ ‐”号。
}
#估计函数####optim####
# debugonce(nlminb)
geyser.res<-nlminb(c(0.5,50,10,80,10),LL,data=waiting,
                   lower=c(0.0001,-Inf,0.0001,
                           -Inf,0.0001),
                   upper=c(0.9999,Inf,Inf,Inf,Inf))
#初始值为 p=0.5,mu1=50,sigma1=10,mu2=80,sigma2=10
#初始值也会被传递给LL
#LL 是被最小化的函数。
#data 是估计用的数据(传递给我们的LL)
#lower 和 upper 分别指定参数的上界和下界。
#查看拟合的参数
geyser.res$par
#拟合的效果
#解释变量
X<-seq(40,120,length=100)
#读出估计的参数
p<-geyser.res$par[1]
mu1<-geyser.res$par[2]
sig1<-geyser.res$par[3]
mu2<-geyser.res$par[4]
sig2<-geyser.res$par[5]
#将估计的参数函数代入原密度函数。
f<-p*dnorm(X,mu1,sig1)+(1-p)*dnorm(X,mu2,sig2)
#作出数据的直方图
hist(waiting,probability=T,col='lightpink3',
     ylab="Density",ylim=c(0,0.04),
     xlab="Eruption waiting times"
       )
#画出拟合的曲线
lines(X,f,col='lightskyblue3',lwd=2)
detach(geyser)

调试的说明

nlminb函数:

nlminb(c(0.5,50,10,80,10),LL,data=waiting,

                   lower=c(0.0001,-Inf,0.0001,

                           -Inf,0.0001),

                   upper=c(0.9999,Inf,Inf,Inf,Inf))

function (start, objective, gradient = NULL, hessian = NULL, 
  ..., scale = 1, control = list(), lower = -Inf, upper = Inf) 
{
  par <- setNames(as.double(start), names(start))
  n <- length(par)
  iv <- integer(78 + 3 * n)
  v <- double(130 + (n * (n + 27))/2)
  .Call(C_port_ivset, 2, iv, v)
  if (length(control)) {
    nms <- names(control)
    if (!is.list(control) || is.null(nms)) 
      stop("'control' argument must be a named list")
    pos <- pmatch(nms, names(port_cpos))
    if (any(nap <- is.na(pos))) {
      warning(sprintf(ngettext(length(nap), "unrecognized control element named %s ignored", 
        "unrecognized control elements named %s ignored"), 
        paste(sQuote(nms[nap]), collapse = ", ")), domain = NA)
      pos <- pos[!nap]
      control <- control[!nap]
    }
    ivpars <- pos <= 4
    vpars <- !ivpars
    if (any(ivpars)) 
      iv[port_cpos[pos[ivpars]]] <- as.integer(unlist(control[ivpars]))
    if (any(vpars)) 
      v[port_cpos[pos[vpars]]] <- as.double(unlist(control[vpars]))
  }
  obj <- quote(objective(.par, ...))

  rho <- new.env(parent = environment())

  assign(".par", par, envir = rho)
  grad <- hess <- low <- upp <- NULL
  if (!is.null(gradient)) {
    grad <- quote(gradient(.par, ...))
    if (!is.null(hessian)) {
      if (is.logical(hessian)) 
        stop("logical 'hessian' argument not allowed.  See documentation.")
      hess <- quote(hessian(.par, ...))
    }
  }
  if (any(lower != -Inf) || any(upper != Inf)) {
    low <- rep_len(as.double(lower), length(par))
    upp <- rep_len(as.double(upper), length(par))
  }
  else low <- upp <- numeric()
  .Call(C_port_nlminb, obj, grad, hess, rho, low, upp, d = rep_len(as.double(scale), 
    length(par)), iv, v)
  iv1 <- iv[1L]
  list(par = get(".par", envir = rho), objective = v[10L], 
    convergence = (if (iv1 %in% 3L:6L) 0L else 1L), iterations = iv[31L], 
    evaluations = c(`function` = iv[6L], gradient = iv[30L]), 
    message = if (19 <= iv1 && iv1 <= 43) {
      if (any(B <- iv1 == port_cpos)) sprintf("'control' component '%s' = %g, is out of range", 
        names(port_cpos)[B], v[iv1]) else sprintf("V[IV[1]] = V[%d] = %g is out of range (see PORT docu.)", 
        iv1, v[iv1])
    } else port_msg(iv1))
}

par最先获得了初始值

obj <- quote(objective(.par, ...))

让obj表示一个名为objective，接受形参.par和...的函数，即我们传入的对数似然函数

rho <- new.env(parent = environment())

创建新环境

assign(".par", par, envir = rho)

将par赋给rho环境中的.par变量

.Call(C_port_nlminb, obj, grad, hess, 

        rho, low, upp, d = rep_len(as.double(scale), length(par)), iv, v)

C_port_nlminb在名为stats.dll的包里，我用reflector和VS都没能看到什么东西,唉，最关键的东西是缺失的。

此时rho包含了参数的初始值，而obj接受.par和多出来的data，即参数初始值和数据。

在nlminb函数的帮助文档里提到：

...

Further arguments to be supplied to objective.

即在这里是传递给对数似然函数的更多参数，这里是data，之所以LL的第一个参数不需要传入是因为我们在源码看到其是.par，也可以在帮助文档看到:

objective

Function to be minimized. Must return a scalar value. The first argument to objective is the vector of parameters to be optimized, whose initial values are supplied through start(start 参数). Further arguments (fixed during the course of the optimization) to objective may be specified as well (see ...(参数)).

---------------------------------------

参考：

http://cos.name/2009/07/maximum-likelihood-estimation-in-r/

原文档：

null

附件列表