Программирование Ньютона Рафсона в R для оценки максимального правдоподобия

Мне нужно запрограммировать метод Ньютона-Рафсона в R для оценки параметра распределения Пуассона. Я только начинаю программировать и работать с 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 принимать только диапазон значений... Другая ошибка связана с условием "если" (критерий остановки), но я действительно не знаю, в чем проблема.


r optimization poisson newtons-method
person Community    schedule 08.03.2017    source источник
comment
Ошибка, которую вы видите, связана с операцией log, которую вы выполняете в функции poisson.lik. В вашем цикле while, когда x0 отрицательное, f(x0) вернет NaN, потому что log отрицательного числа невозможно, давая вам NaN (не число). Операция над NaN даст вам еще один NaN, так что ваш оператор if потерпит неудачу. (например, попробуйте if(NaN < 0) 1+1)   -  person ialm    schedule 09.03.2017

Ответы (1)


arrow_upward
0
arrow_downward

Ваша проблема с NaN исходила из вашей функции poisson.lik. Вам нужно иметь 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)

Основная проблема, которую я увидел, заключалась в том, что вам нужно знать, куда вы хотите, чтобы функция перешла. Обычно он используется, когда вы знаете какой-то аспект данных. Например, если бы вы знали, что вам нужно, логарифмическая вероятность была бы асимметрией или каким-либо типом статистики, чтобы вы могли найти параметр, который дает ту же статистику. Приведенный выше код решает проблему. Надеюсь это поможет.

person Rex    schedule 08.02.2019