1

尺度参数为 1 的 gamma 分布的对数似然可以写为：

(α−1)s−nlogΓ(α) 其中 alpha 是形状参数，s=∑logXi是充分的统计量。

随机抽取 n = 30 的样本，形状参数为 alpha = 4.5。使用newton_search和make_derivative，找到 alpha 的最大似然估计。使用 alpha 的矩估计，即 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)，我得到的答案应该是双倍的。

4

1 回答 1

0

我重新编写了代码，它现在可以完美运行（甚至比我最初编写的更好）。感谢所有帮助过的人。:-)

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)
于 2015-07-15T15:54:56.667 回答