title: How to enable optimization by optim() return close estimate?
body: <p>First, I need to clarify that I have read the following posts but my problem still cant be solved:</p>
  
  <ol>
  <li><p><a href="https://stackoverflow.com/questions/39747569/r-optim-l-bfgs-b-needs-finite-values-of-fn-weibull">R optim() L-BFGS-B needs finite values of &#39;fn&#39; - Weibull</a></p></li>
  <li><p><a href="https://stackoverflow.com/questions/41541528/optimization-of-optim-in-r-l-bfgs-b-needs-finite-values-of-fn">Optimization of optim() in R ( L-BFGS-B needs finite values of &#39;fn&#39;)</a></p></li>
  <li><p><a href="https://stackoverflow.com/questions/46462483/r-optimize-multiple-parameters">R optimize multiple parameters</a></p></li>
  <li><p><a href="https://stackoverflow.com/questions/45171900/optim-function-with-infinite-value">optim function with infinite value</a></p></li>
  </ol>
  
  <p>Below are the code to do simulation and proceed maximum likelihood estimation.</p>
  
  <pre><code>    #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 &lt;- function(par) {
            h.n &lt;- rep(0,n)
            a0 &lt;- par[1]  
            a1 &lt;- par[2] 
            g1 &lt;- par[3]
            b1 &lt;- par[4]
            d1 &lt;- 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 &lt;- 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)) 
  </code></pre>
  
  <p>Then I will get the result of:</p>
  
  <pre><code>&gt; iter   10 value 3172.782149 
  
  &gt; iter   20 value 3172.371186 
  
  &gt; iter   30 value 3171.952137 
  
  &gt; iter   40 value 3171.525942 
  
  &gt; iter   50 value 3171.174571 
  
  &gt; iter   60 value 3171.095186 
  
  &gt; Error in optim(par = c(0.1, 0.01, 0.1, 0.01, 
  &gt; 0.1), ll, method = "L-BFGS-B",  :    L-BFGS-B 
  &gt; needs finite values of 'fn'
  </code></pre>
  
  <p>So I try to change the lower bound, and it returns</p>
  
  <pre><code>&gt; &gt; 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))
  &gt; 
  &gt; iter   10 value 3172.782149 
  &gt; 
  &gt; iter   20 value 3172.371186 
  &gt; 
  &gt; iter   30 value 3171.952137   
  &gt;
  &gt; iter   40 value 3171.525942   
  &gt;
  &gt; iter   50 value 3171.174571 
  &gt; 
  &gt; iter   60 value 3171.095186 
  &gt; 
  &gt; iter   70 value 3171.076036 
  &gt; 
  &gt; iter   80 value 3171.044809 
  &gt; 
  &gt; iter   90 value 3171.014010 
  &gt; 
  &gt; iter  100 value 3170.991805 
  &gt; 
  &gt; iter  110 value 3170.971857 
  &gt; 
  &gt; iter  120 value 3170.954827 
  &gt; 
  &gt; iter  130 value 3170.941397 
  &gt; 
  &gt; iter  140 value 3170.925935 
  &gt; 
  &gt; iter  150 value 3170.915694  
  &gt; 
  &gt; iter  160 value 3170.904309 
  &gt; 
  &gt; iter  170 value 3170.894642
  
  &gt; iter  180 value 3170.887122  
  
  &gt; iter  190 value 3170.880802 
  &gt; 
  &gt; iter  200 value 3170.874319 
  &gt; 
  &gt; iter  210 value 3170.870006 
  &gt; 
  &gt; iter  220 value 3170.866008 
  &gt; 
  &gt; iter  230 value 3170.865497 
  &gt; 
  &gt; final  value 3170.865422  converged  
  &gt;
  &gt; $`par` [1] 3.242429e+05
  &gt; 2.691999e-04 3.896417e-01 6.174022e-04 2.626361e+01
  &gt; 
  &gt; $value [1] 3170.865
  &gt; 
  &gt; $counts function gradient 
  &gt;      291      291 
  &gt; 
  &gt; $convergence [1] 0
  &gt; 
  &gt; $message [1] "CONVERGENCE: REL_REDUCTION_OF_F &lt;= FACTR*EPSMCH"
  </code></pre>
  
  <p>Definitely, the estimated parameters are far from the true value.</p>
  
  <p>What could I do to get close estimates to true value?</p>