Vfrdtyky

First, I need to clarify that I have read the following posts but my problem still cant be solved:

1. R optim() L-BFGS-B needs finite values of 'fn' - Weibull




2. Optimization of optim() in R ( L-BFGS-B needs finite values of 'fn')




3. R optimize multiple parameters




4. optim function with infinite value

Below are the code to do simulation and proceed maximum likelihood estimation.

    #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)) 

Then I will get the result of:

> 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'

So I try to change the lower bound, and it returns

> > 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"

Definitely, the estimated parameters are far from the true value.

What could I do to get close estimates to true value?

When an MLE is far from the true value, there are several possible explanations:

1. You don't have enough data to get an accurate estimate. Try using a much larger sample size and see if thing come out closer.



2. 
   

   You have coded the likelihood incorrectly. This is harder to diagnose; basically you just want to read it over and check your coding.

   
   
   
   
   
   • 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. Your likelihood doesn't have a unique maximum, because the parameters are not identifiable.

There are probably others, too.