How to find the second derivative in R and while using newton's method with numerical derivation -

the log-likelihood of gamma distribution scale parameter 1 can written as:

(α−1)s−nlogΓ(α) alpha shape parameter , s=∑logxi sufficient statistic.

randomly draw sample of n = 30 shape parameter of alpha = 4.5. using newton_search , make_derivative, find maximum likelihood estimate of alpha. use moment estimator of alpha, i.e., mean of x initial guess. log-likelihood function in r is:

x <- rgamma(n=30, shape=4.5) gllik <- function() {   s <- sum(log(x))   n <- length(x)   function(a) {     (a - 1) * s - n * lgamma(a)     } } 

i have created make_derivative function follows:

make_derivative <- function(f, h) {   (f(x + h) - f(x - h)) / (2*h) } 

i have created newton_search function incorporates make_derivative function; however, need use newton_search on second derivative of log-likelihood function , i'm not sure how fix following code in order that:

newton_search2 <- function(f, h, guess, conv=0.001) {     set.seed(2)       y0 <- guess     n = 1000     <- 1; y1 <- y0     p <- numeric(n)   while (i <= n) {     make_derivative <- function(f, h) {   (f(y0 + h) - f(y0 - h)) / (2*h)     }     y1 <- (y0 - (f(y0)/make_derivative(f, h)))     p[i] <- y1     <- + 1     if (abs(y1 - y0) < conv) break     y0 <- y1   }   return (p[(i-1)]) } 

hint: must apply newton_search first , second derivatives (derived numerically using make_derivative) of log-likelihood. answer should near 4.5.

when run newton_search2(gllik(), 0.0001, mean(x), conv = 0.001), double answer should be.

i re-wrote code , works (even better had wrote). helped. :-)

newton_search <- function(f, df, guess, conv=0.001) {     set.seed(1)     y0 <- guess     n = 100     <- 1; y1 <- y0     p <- numeric(n)   while (i <= n) {     y1 <- (y0 - (f(y0)/df(y0)))     p[i] <- y1     <- + 1     if (abs(y1 - y0) < conv) break     y0 <- y1   }   return (p[(i-1)]) }  make_derivative <- function(f, h) {   function(x){(f(x + h) - f(x - h)) / (2*h)   } }  df1 <- make_derivative(gllik(), 0.0001) df2 <- make_derivative(df1, 0.0001) newton_search(df1, df2, mean(x), conv = 0.001)