Матрица Гессе в optim () в R

У меня возникли проблемы с использованием optim () в R для определения вероятности, включающей интеграл, и получения матрицы Гесси для двух параметров. Алгоритм сходился, но я получаю сообщение об ошибке, когда использую опцию hessian = TRUE в optim (). Ошибка:

Ошибка интегрирования (подынтегральное выражение1, нижний = s1 [i] - 1, верхний = s1 [i]): неконечное значение функции.

Также было предупреждение о НП

Вот мой код:

s1=c(1384,1,1219,1597,2106,145,87,1535,290,1752,265,588,1188,160,745,237,479,39,99,56,1503,158,916,651,1064,166,635,19,553,51,79,155,85,1196,142,108,325  
 ,135,28,422,1032,1018,128,787,1704,307,854,6,896,902)

LLL=function (par) {

  integrand1 <- function(x){ (x-s1[i]+1)*dgamma(x, shape=par[1], rate=par[2]) }
  integrand2 <- function(x){ (-x+s1[i]+1)*dgamma(x, shape=par[1],rate=par[2]) }



  likelihood = vector() 

  for(i in 1:length(s1)) {likelihood[i] = 
    log( integrate(integrand1,lower=s1[i]-1,upper=s1[i])$value+ integrate(integrand2,lower=s1[i],upper=s1[i]+1)$value )  
  }

  like= -sum(likelihood)
  return(like)

}




optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0))
optim(par=c(0.1,0.1),LLL,method="L-BFGS-B", lower=c(0,0), hessian=TRUE)

Спасибо за вашу помощь!


r optimization integral gamma-distribution nonlinear-optimization
person Y. Ma    schedule 09.01.2017    source источник
comment
Я полагаю, вы оцениваете dgamma(0, shape=.1, scale=.1), что равно Inf, при интеграции для i == 2 от s1[2]-1 (0) до s1[2] (1). Что произойдет, если вы исключите 1 из s1?   -  person Thales    schedule 09.01.2017
comment
@Thales Спасибо.   -  person Y. Ma    schedule 09.01.2017
comment
Вчера вам сказали использовать lower=c(.001,.001). Почему вы не используете это в качестве аргумента lower? Измените lower на c(0.01,0.01), и вы увидите гессен.   -  person Bhas    schedule 09.01.2017
comment
@Bhas Да, я использовал lower = c (0,01, 0,01) в соответствии с вашим предложением (скопировал неправильный код в этом сообщении). Есть идеи, почему параметр скорости в гамме всегда сходится к нижней границе и верна ли оценка? Спасибо.   -  person Y. Ma    schedule 09.01.2017

Ответы (1)


arrow_upward
0
arrow_downward

optim минимизирует функцию. Вы можете построить функцию правдоподобия с аргументом rate. Чтобы получить сюжет, нужно немного повозиться. Делай это так:

z2 <- function(rate) {
    par <- numeric(2)
    par[1] <- .68
    par[2] <- rate
    y <- LLL(par)
    y
}

z1 <- Vectorize(z2,vectorize.args="rate")

curve(z1, from=.001,to=1)

Вы увидите, что функция минимальна для самого низкого значения rate. То же самое, если вы измените from на .1. Я не могу судить, верна ли оценка.

person Bhas    schedule 09.01.2017
comment
Еще раз спасибо за вашу огромную помощь. - person Y. Ma; 09.01.2017