Как включить оптимизацию с помощью optim(), возвращающей близкую оценку?

Во-первых, мне нужно уточнить, что я прочитал следующие сообщения, но моя проблема все еще не может быть решена:

  1. R optim() L-BFGS-B нуждается конечные значения 'fn' - Вейбулл

  2. Оптимизация optim() в R ( L-BFGS-B требует конечных значений 'fn')

  3. R оптимизировать несколько параметров

  4. оптимальная функция с бесконечным значением

Ниже приведен код для моделирования и оценки максимального правдоподобия.

    #simulation
    #a0, a1, g1, b1 and d1 are my parameters
    #set true value of parameters to 
    #simulate a set of data with size 2000
    #x is the simulated data sets

    set.seed(5)
    a0 = 2.3; a1 = 0.05; g1 = 0.68; b1 = 
    0.09; d1 = 2.0; n=2000

    x = h = rep(0, n)

    h[1] = 6
    x[1] = rpois(1,h[1])

     for (i in 2:n) {

      h[i] = (a0 + a1 *
            (abs(x[i-1]-h[i-1])-g1*(x[i-1]- 
            h[i-1]))^d1 +
            b1 * (h[i-1]^d1))^(1/d1)
      x[i] = rpois(1,h[i])
    }

      #this is my log-likelihood function
       ll <- function(par) {
          h.n <- rep(0,n)
          a0 <- par[1]  
          a1 <- par[2] 
          g1 <- par[3]
          b1 <- par[4]
          d1 <- par[5]

          h.n[1] = x[1]
          for (i in 2:n) {

           h.n[i] = (a0 + a1 *
                 (abs(x[i-1]-h.n[i-1])-g1* 
                  (x[i-1]-h.n[i-1]))^d1 +
                  b1 * (h.n[i-1]^d1))^(1/d1)
            }
           -sum(dpois(x, h.n, log=TRUE))
            }

         #as my true value are a0 = 2.3; a1 
         #= 0.05; g1 = 0.68; b1 = 0.09; d1 
         #= 2.0 
         #I put the parscale to become 
         #c(1,0.01,0.1,0.01,1)
       ps <- c(1.0, 1e-02, 1e-01, 1e-02,1.0)

         #optimization to check whether 
         #estimate return near to the true 
         #value
         optim(par=c(0.1,0.01,0.1,0.01,0.1), 
          ll, method = "L-BFGS-B",
          lower=c(1e-6,-10,-10,-10, 1e- 6),
          control= list(maxit=1000,
          parscale=ps,trace=1)) 

Тогда я получу результат:

> iter   10 value 3172.782149 

> iter   20 value 3172.371186 

> iter   30 value 3171.952137 

> iter   40 value 3171.525942 

> iter   50 value 3171.174571 

> iter   60 value 3171.095186 

> Error in optim(par = c(0.1, 0.01, 0.1, 0.01, 
> 0.1), ll, method = "L-BFGS-B",  :    L-BFGS-B 
> needs finite values of 'fn'

Итак, я пытаюсь изменить нижнюю границу, и она возвращает

> > optim(par=c(0.1,0.01,0.1,0.01,0.1), ll, method = "L-BFGS-B",lower=c(1e-6,1e-6,-10,1e-6,1e-6),control=list(maxit=1000,parscale=ps,trace=1))
> 
> iter   10 value 3172.782149 
> 
> iter   20 value 3172.371186 
> 
> iter   30 value 3171.952137   
>
> iter   40 value 3171.525942   
>
> iter   50 value 3171.174571 
> 
> iter   60 value 3171.095186 
> 
> iter   70 value 3171.076036 
> 
> iter   80 value 3171.044809 
> 
> iter   90 value 3171.014010 
> 
> iter  100 value 3170.991805 
> 
> iter  110 value 3170.971857 
> 
> iter  120 value 3170.954827 
> 
> iter  130 value 3170.941397 
> 
> iter  140 value 3170.925935 
> 
> iter  150 value 3170.915694  
> 
> iter  160 value 3170.904309 
> 
> iter  170 value 3170.894642

> iter  180 value 3170.887122  

> iter  190 value 3170.880802 
> 
> iter  200 value 3170.874319 
> 
> iter  210 value 3170.870006 
> 
> iter  220 value 3170.866008 
> 
> iter  230 value 3170.865497 
> 
> final  value 3170.865422  converged  
>
> $`par` [1] 3.242429e+05
> 2.691999e-04 3.896417e-01 6.174022e-04 2.626361e+01
> 
> $value [1] 3170.865
> 
> $counts function gradient 
>      291      291 
> 
> $convergence [1] 0
> 
> $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Определенно, расчетные параметры далеки от истинного значения.

Что я могу сделать, чтобы приблизить оценки к истинному значению?


person Miyazaki    schedule 16.11.2018    source источник


Ответы (1)


Когда MLE далек от истинного значения, есть несколько возможных объяснений:

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

  2. Вы неправильно закодировали вероятность. Это сложнее диагностировать; в основном вы просто хотите прочитать его и проверить свой код.

    • I'm not familiar with your model, but this looks likely in your case: in your simulation, h[1] is always 6 and x[1] is a random value with that mean; in your likelihood, you're assuming that h[1] is equal to x[1]. That's unlikely to be true.
  3. Ваша вероятность не имеет уникального максимума, потому что параметры не идентифицируемы.

Вероятно, есть и другие.

person user2554330    schedule 16.11.2018
comment
Спасибо за объяснение, @user2554330 1. Я постараюсь увеличить размер выборки. 2. В вопросе укажу модель. Я поставил 6 в качестве начального значения h[1], чтобы он мог начать генерировать конечные данные. С вероятностью я положил x[1] равным h.n[1], так что при подгонке модели к реальным данным, которые могут содержать выбросы, первый h.n[1] будет точно таким же, как и первая точка данных на графике. Однако можно также положить h.n[1] равным нулю. 3. Я также предполагаю, что он не имеет уникального максимума, но что означает неидентифицируемые параметры? - person Miyazaki; 18.11.2018
comment
Это выходит за рамки темы для stackoverflow, но в ответ на ваш вопрос: я бы рассматривал h[1] как параметр, который нужно оценить, или известное значение. Установка его на ноль была бы плохой, потому что это говорит, что x[1] всегда будет равно нулю, но, по-видимому, это не так. Неидентифицируемый означает, что разные значения параметров дают точно такое же распределение данных, поэтому уникального MLE не существует. - person user2554330; 18.11.2018