title: Fitting Gamma distribution to data in R using optim, ML
body: <p>Im kinda new to R. I have a dataset, which also includes data of family income and I have to fit a Gamma distribution to this data, using the Maximum Likelihood Estimates. It is specifically told that we need to use the package optim, and not fitdistr. So this is my code: </p>
  
  <pre><code>t1 &lt;- sum(log(newdata$faminc)) 
  t2 &lt;- sum(newdata$faminc)
  obs &lt;- nrow(newdata)
  lh.gamma &lt;- function(par) {
    -((par[1]-1)*t1 - par[2]*t2 - obs*par[1]*log(par[2]) - obs*lgamma(par[1]))
  }
  
  #initial guess for a = mean^2(x)/var(x) and b = mean(x) / var(x) 
  a1 &lt;- (mean(newdata$faminc))^2/var(newdata$faminc)
  b1 &lt;- mean(newdata$faminc)/var(newdata$faminc)
  
  init &lt;- c(a1,b1)
  q &lt;- optim(init, lh.gamma, method = "BFGS")
  q
  </code></pre>
  
  <p>Also tried filling in just values in the init vector, and including this piece of code; </p>
  
  <pre><code>  dlh.gamma &lt;- function(par){
    cbind(obs*digamma(par[1])+obs*log(par[2])-t2,
       obs*par[1]/par[2]-1/par[2]^2*t1)
  }
  </code></pre>
  
  <p>and then the optim would look like :</p>
  
  <pre><code> q &lt;- optim(init, lh.gamma, dhl.gamma, method="BFGS")
  </code></pre>
  
  <p>None of it 'works'. First, when I tried the code at school computers, it gave me very huge numbers for the shape and rate parameters, which was not possible. Now, trying at home, I get this:</p>
  
  <pre><code>&gt; q &lt;- optim(init, lh.gamma, method = "BFGS")
  Error in optim(init, lh.gamma, method = "BFGS") : 
    non-finite finite-difference value [2]
  In addition: There were 50 or more warnings (use warnings() to see the first 50)
  &gt; q
  function (save = "default", status = 0, runLast = TRUE) 
  .Internal(quit(save, status, runLast))
  &lt;bytecode: 0x000000000eaac960&gt;
  &lt;environment: namespace:base&gt;
  </code></pre>
  
  <p>q is not even 'created'. Except for when I include the dlh.gamma part above, but then I just get huge numbers again and no convergence. </p>
  
  <p>Anybody who knows what goes wrong/what to do? </p>
  
  <p>Edit: </p>
  
  <pre><code>&gt; dput(sample(newdata$faminc, 500))
  c(42.5, 87.5, 22.5, 17.5, 12.5, 30, 30, 17.5, 42.5, 62.5, 62.5, 
  30, 30, 150, 22.5, 30, 42.5, 30, 17.5, 8.75, 42.5, 42.5, 42.5, 
  62.5, 42.5, 30, 17.5, 87.5, 62.5, 150, 42.5, 150, 42.5, 42.5, 
  42.5, 6.25, 62.5, 87.5, 6.25, 87.5, 30, 150, 22.5, 62.5, 42.5,    
  150, 17.5, 42.5, 42.5, 42.5, 62.5, 22.5, 42.5, 42.5, 30, 62.5, 
  30, 62.5, 87.5, 87.5, 42.5, 22.5, 62.5, 22.5, 8.75, 30, 30, 17.5, 
  87.5, 8.75, 62.5, 30, 17.5, 22.5, 62.5, 42.5, 30, 17.5, 62.5, 
  8.75, 62.5, 42.5, 150, 30, 62.5, 87.5, 17.5, 62.5, 30, 62.5, 
  87.5, 42.5, 62.5, 30, 62.5, 42.5, 87.5, 150, 12.5, 42.5, 62.5, 
  42.5, 62.5, 62.5, 150, 30, 87.5, 12.5, 17.5, 42.5, 62.5, 30, 
  6.25, 62.5, 42.5, 12.5, 62.5, 8.75, 17.5, 42.5, 62.5, 87.5, 8.75, 
  62.5, 30, 62.5, 87.5, 42.5, 62.5, 62.5, 12.5, 150, 42.5, 62.5,  
  12.5, 62.5, 42.5, 62.5, 62.5, 87.5, 42.5, 62.5, 30, 42.5, 150, 
  42.5, 30, 62.5, 62.5, 87.5, 42.5, 30, 62.5, 62.5, 42.5, 42.5, 
  30, 62.5, 42.5, 42.5, 62.5, 62.5, 150, 42.5, 30, 42.5, 62.5, 
  17.5, 62.5, 17.5, 150, 8.75, 62.5, 30, 62.5, 42.5, 42.5, 22.5, 
  150, 62.5, 42.5, 62.5, 62.5, 22.5, 30, 62.5, 30, 150, 42.5, 42.5, 
  42.5, 62.5, 30, 12.5, 30, 150, 12.5, 8.75, 22.5, 30, 22.5, 30, 
  42.5, 42.5, 42.5, 30, 12.5, 62.5, 42.5, 30, 22.5, 42.5, 87.5, 
  22.5, 12.5, 42.5, 62.5, 62.5, 62.5, 30, 42.5, 30, 62.5, 30, 62.5, 
  12.5, 22.5, 42.5, 22.5, 87.5, 30, 22.5, 17.5, 42.5, 62.5, 17.5, 
  250, 150, 42.5, 30, 42.5, 30, 62.5, 17.5, 87.5, 22.5, 150, 62.5, 
  42.5, 6.25, 87.5, 62.5, 42.5, 30, 42.5, 62.5, 42.5, 87.5, 62.5, 
  150, 42.5, 30, 6.25, 22.5, 30, 42.5, 42.5, 62.5, 250, 8.75, 150, 
  42.5, 30, 42.5, 30, 42.5, 42.5, 30, 30, 150, 22.5, 62.5, 30, 
  8.75, 150, 62.5, 87.5, 150, 42.5, 30, 42.5, 42.5, 42.5, 30, 8.75, 
  42.5, 42.5, 30, 22.5, 62.5, 17.5, 62.5, 62.5, 42.5, 8.75, 42.5, 
  12.5, 12.5, 150, 42.5, 42.5, 17.5, 42.5, 62.5, 62.5, 42.5, 42.5, 
  30, 42.5, 62.5, 30, 62.5, 42.5, 42.5, 42.5, 22.5, 62.5, 62.5, 
  62.5, 22.5, 150, 62.5, 42.5, 62.5, 42.5, 30, 30, 62.5, 22.5, 
  62.5, 87.5, 62.5, 42.5, 42.5, 22.5, 62.5, 62.5, 30, 42.5, 42.5, 
  8.75, 87.5, 42.5, 42.5, 87.5, 30, 62.5, 17.5, 62.5, 42.5, 17.5, 
  22.5, 62.5, 8.75, 62.5, 22.5, 22.5, 22.5, 42.5, 17.5, 22.5, 62.5, 
  42.5, 42.5, 42.5, 42.5, 42.5, 30, 30, 8.75, 30, 42.5, 62.5, 22.5, 
  6.25, 30, 42.5, 62.5, 17.5, 62.5, 42.5, 8.75, 22.5, 30, 17.5, 
  22.5, 62.5, 42.5, 150, 87.5, 22.5, 12.5, 62.5, 62.5, 62.5, 30, 
  42.5, 22.5, 62.5, 87.5, 30, 42.5, 62.5, 22.5, 87.5, 30, 30, 22.5, 
  87.5, 87.5, 250, 30, 62.5, 250, 62.5, 42.5, 42.5, 62.5, 62.5, 
  42.5, 6.25, 62.5, 62.5, 62.5, 42.5, 42.5, 150, 62.5, 62.5, 30, 
  150, 22.5, 87.5, 30, 150, 17.5, 8.75, 62.5, 42.5, 62.5, 150, 
  42.5, 22.5, 42.5, 42.5, 17.5, 62.5, 17.5, 62.5, 42.5, 150, 250, 
  22.5, 42.5, 30, 62.5, 62.5, 42.5, 42.5, 30, 150, 150, 42.5, 17.5, 
  17.5, 42.5, 8.75, 62.5, 42.5, 42.5, 22.5, 150, 62.5, 30, 250, 
  62.5, 87.5, 62.5, 8.75, 62.5, 30, 30, 8.75, 17.5, 17.5, 150, 
  22.5, 62.5, 62.5, 42.5)
  </code></pre>
  
  <p>The faminc variable is in 1000s</p>
  
  <p>Edit2: </p>
  
  <p>Okay, the code is good, but now I try to fit the distribution over the histogram using the following: </p>
  
  <pre><code>x &lt;- rgamma(500,shape=q$par[1],scale=q$par[2])
  hist(newdata$faminc, prob = TRUE)
  curve(dgamma(x, shape=q$par[1], scale=q$par[2]), add=TRUE, col='blue') 
  </code></pre>
  
  <p>It just produces a flat blue line at the x-axis.. </p>