title: How to find the second derivative in R and while using newton's method with numerical derivation
body: <p>The log-likelihood of the gamma distribution with scale parameter 1 can be written as:</p>
  
  <p><code>(α−1)s−nlogΓ(α)</code>
  where alpha is the shape parameter and <code>s=∑logXi</code> is the sufficient statistic.</p>
  
  <p>Randomly draw a sample of n = 30 with a shape parameter of alpha = 4.5. Using <code>newton_search</code> and <code>make_derivative</code>, find the maximum likelihood estimate of alpha. Use the moment estimator of alpha, i.e., mean of x as the initial guess. The log-likelihood function in R is:</p>
  
  <pre><code>x &lt;- rgamma(n=30, shape=4.5)
  gllik &lt;- function() {
    s &lt;- sum(log(x))
    n &lt;- length(x)
    function(a) {
      (a - 1) * s - n * lgamma(a)
      }
  }
  </code></pre>
  
  <p>I have created the make_derivative function as follows:</p>
  
  <pre><code>make_derivative &lt;- function(f, h) {
    (f(x + h) - f(x - h)) / (2*h)
  }
  </code></pre>
  
  <p>I also have created a <code>newton_search</code> function that incorporates the <code>make_derivative</code> function; However, I need to use <code>newton_search</code> on the second derivative of the log-likelihood function and I'm not sure how to fix the following code in order for it to do that:</p>
  
  <pre><code>newton_search2 &lt;- function(f, h, guess, conv=0.001) {
      set.seed(2)  
      y0 &lt;- guess
      N = 1000
      i &lt;- 1; y1 &lt;- y0
      p &lt;- numeric(N)
    while (i &lt;= N) {
      make_derivative &lt;- function(f, h) {
    (f(y0 + h) - f(y0 - h)) / (2*h)
      }
      y1 &lt;- (y0 - (f(y0)/make_derivative(f, h)))
      p[i] &lt;- y1
      i &lt;- i + 1
      if (abs(y1 - y0) &lt; conv) break
      y0 &lt;- y1
    }
    return (p[(i-1)])
  }
  </code></pre>
  
  <p>Hint: You must apply <code>newton_search</code> to the first and second derivatives (derived numerically using <code>make_derivative</code>) of the log-likelihood. Your answer should be near 4.5.</p>
  
  <p>when I run <code>newton_search2(gllik(), 0.0001, mean(x), conv = 0.001)</code>, I get double what the answer should be.</p>