Как найти вторую производную в R и при этом использовать метод Ньютона с численным выводом

Логарифмическая вероятность гамма-распределения с масштабным параметром 1 может быть записана как:

(α−1)s−nlogΓ(α), где альфа — параметр формы, а s=∑logXi — достаточная статистика.

Случайным образом нарисуйте образец n = 30 с параметром формы альфа = 4,5. Используя newton_search и make_derivative, найдите оценку максимального правдоподобия альфа. Используйте моментную оценку альфы, т. е. среднее значение x в качестве начального предположения. Функция логарифмического правдоподобия в R:

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

Я создал функцию make_derivative следующим образом:

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

Я также создал функцию newton_search, которая включает в себя функцию make_derivative; Однако мне нужно использовать newton_search во второй производной функции логарифмического правдоподобия, и я не уверен, как исправить следующий код, чтобы он это сделал:

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)])
}

Подсказка: вы должны применить newton_search к первой и второй производным (полученным численно с использованием make_derivative) логарифмического правдоподобия. Ваш ответ должен быть около 4,5.

когда я запускаю newton_search2(gllik(), 0.0001, mean(x), conv = 0.001), я получаю ответ вдвое больше, чем должен быть.


r derivative gamma-distribution newtons-method
person user114634    schedule 14.07.2015    source источник
comment
Когда я запускаю newton_search2(gllik(), 0,0001, mean(x), conv = 0,001), я получаю в два раза больше, чем должен быть ответ. Я должен получить около 4,5, но я получаю около 9,2. Я не уверен, где мой код пошел не так.   -  person user114634    schedule 14.07.2015
comment
Вы не предоставляете никаких доказательств того, что вы прочитали и предприняли попытку исправить проблемы, выявленные RHertel. Также вызывало недоумение то, что у функции правдоподобия не было аргумента данных.   -  person IRTFM    schedule 14.07.2015
comment
Уравнение правдоподобия мне дали так. Я не думаю, что я должен изменить эту функцию.   -  person user114634    schedule 15.07.2015
comment
Также у производной функции было требование замыкания функции, которое запоминает f и h.   -  person user114634    schedule 15.07.2015

Ответы (1)


arrow_upward
0
arrow_downward

Я переписал код, и теперь он отлично работает (даже лучше, чем то, что я написал изначально). Спасибо всем, кто помог. :-)

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)
person user114634    schedule 15.07.2015
comment
Ницца. Но не забудьте проверить, что то, что вы нашли, является максимальным. т.е. подключите решение, найденное в df1 и df2, и проверьте знак df2(solution). - person Bhas; 15.07.2015