在 R 中编程 Newton Raphson 以进行最大似然估计答案

【问题标题】：Programming Newton Raphson in R for Maximum Likelihood estimation在 R 中编程 Newton Raphson 以进行最大似然估计
【发布时间】：2017-07-29 17:40:44
【问题描述】：

我需要在 R 中编写 Newton-Raphson 方法来估计泊松分布的参数。我刚刚开始使用编程和 R。当我使用模拟数据运行程序时，R 返回一些错误。

Error in if (abs(x1 - x0) < stoptol) break : 
missing value where TRUE/FALSE needed
In addition: Warning messages:
1: In log(mu) : NaNs produced
2: In log(mu) : NaNs produced
3: In log(mu) : NaNs produced

这是我目前所拥有的：

#
#  NEWTON-RAPHSON METHOD 
#

#generate the data 
lambda=3.2
y=rpois(500,lambda)

#declare the log likehood function
poisson.lik<-function(mu,ydata=y){
   n<-length(ydata)
   logl<-sum(ydata)*log(mu)-n*mu-sum(lfactorial(ydata))
   return(-logl)
}

## Newton-Raphson 

NR<-function(initval, f, stoptol=1e-05, imax=25){
   i=0
   h=1e-05
   x0=initval-0.1
   x1=initval

   while(i<=imax){
     df.dx=(f(x0+h)-f(x0))/h
     x1=(x0-(f(x0)/df.dx))
     i=i+1
     if(abs(x1-x0)<stoptol) break
     x0=x1
    }
    list(nstep=i, initial=initval, final=x1, fctval=f(x1))
}

NR(initval=3,poisson.lik)

据我了解，一个问题来自参数 mu 在 NR 函数的迭代和对数计算中所取的值。也许我应该强制 mu 只采用一系列值... 另一个错误是关于条件“如果”（停止标准），但我真的不知道问题是什么。

【问题讨论】：

您看到的错误与您在poisson.lik 函数中执行的log 操作有关。在你的“while”循环中，当x0 为负数时，f(x0) 将返回NaN，因为负数的log 是不可能的，给你NaN（不是数字）。对NaN 的操作将给您另一个NaN，因此您的if 语句将失败。（例如尝试if(NaN < 0) 1+1）

标签： r optimization poisson newtons-method

【解决方案1】：

您的 NaN 问题来自您的 poisson.lik 函数。在 mu 为负数的情况下，您需要有 log(abs(mu))。

#
#  NEWTON-RAPHSON METHOD 
#

#generate the data 
lambda <- 3.2
ydata <- rpois(500, lambda)

#declare the log likehood function
poisson.lik <- function(mu = ""){
    n <- length(ydata)
    loglik <- -n * mu - sum(log(factorial(ydata))) + log(abs(mu)) * sum(ydata) 
    return(-loglik)
}

#creating the Newton Raphson loop
NR <- function ( mu = "", initval = "", f = "", stoptol = 1e-05, imax = "") {
    i = 0
    h = .1
    mu0 = initval - 0.1
    mu1 = initval
    df.dx <- double(1) #predeclare your variable types     
    while ( abs(f(mu) - f(mu1)) > stoptol && i <= imax ) {
        mu0 <- mu1
        df.dx <- (f(mu0 + h) - f(mu0)) / h
        mu1 <- mu0 + (mu - mu1) / abs(df.dx)
        i <- i + 1        
    }
    return(list("nstep" = i, "Final" = mu1, "fctval" = f(mu1)))
} 

mu <- 3.2
initval <- 1
f <- poisson.lik
stoptol <- 0.0001
imax <- 10^4
NR(mu, initval, f, stoptol, imax)

我看到的主要问题是你需要知道你想让函数去哪里。当您了解数据的某些方面时，通常会使用它。例如，如果您知道您想要的对数似然、偏度或某种类型的统计量，那么您就可以求解产生相同统计量的参数。上面的代码解决了这个问题。希望这可以帮助。

【讨论】：