1

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

(α−1)s−nlogΓ(α) where alpha is the shape parameter and s=∑logXi is the sufficient statistic.

Randomly draw a sample of n = 30 with a shape parameter of alpha = 4.5. Using newton_search and make_derivative, 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:

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 the make_derivative function as follows:

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

I also have created a newton_search function that incorporates the make_derivative function; However, I need to use newton_search 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:

newton_search2 <- function(f, h, guess, conv=0.001) {
    set.seed(2)  
    y0 <- guess
    N = 1000
    i <- 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
    i <- i + 1
    if (abs(y1 - y0) < conv) break
    y0 <- y1
  }
  return (p[(i-1)])
}

Hint: You must apply newton_search to the first and second derivatives (derived numerically using make_derivative) of the log-likelihood. Your answer should be near 4.5.

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

dpel
  • 1,954
  • 1
  • 21
  • 31
user114634
  • 11
  • 4
  • When I run newton_search2(gllik(), 0.0001, mean(x), conv = 0.001) I get double of what the answer should be. I should get about 4.5, but I get about 9.2. I'm not sure where my code has gone wrong. – user114634 Jul 14 '15 at 16:47
  • 1
    You show no evidence that you have read and made an effort to correct the problems that RHertel has identified. It was also puzzling that the likelihood function had no data argument. – IRTFM Jul 14 '15 at 18:20
  • The likelihood equation was given to me like that. I don't think I'm supposed to change that function. – user114634 Jul 15 '15 at 13:26
  • Also the derivative function had the requirement of a function closure that remembers f and h – user114634 Jul 15 '15 at 13:27

1 Answers1

0

I re-wrote the code and it works perfectly now (even better than what I had originally wrote). Thanks to all who helped. :-)

newton_search <- function(f, df, guess, conv=0.001) {
    set.seed(1)
    y0 <- guess
    N = 100
    i <- 1; y1 <- y0
    p <- numeric(N)
  while (i <= N) {
    y1 <- (y0 - (f(y0)/df(y0)))
    p[i] <- y1
    i <- i + 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)
user114634
  • 11
  • 4
  • Nice. But don't forget to check that what you found is a maximum. I.e. plug the solution found in `df1` and `df2` and check the sign of `df2(solution)`. – Bhas Jul 15 '15 at 20:01