1

スケール パラメーター 1 のガンマ分布の対数尤度は、次のように記述できます。

(α−1)s−nlogΓ(α) ここで、alpha は形状パラメーターでs=∑logXiあり、十分な統計量です。

n = 30 のサンプルを、形状パラメーター alpha = 4.5 でランダムに描画します。とを使用newton_searchしてmake_derivative、アルファの最尤推定を見つけます。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対数尤度関数の 2 次導関数を使用する必要があり、それを行うために次のコードを修正する方法がわかりません。

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の 1 次および 2 次導関数 ( を使用して数値的に導出) に適用する必要があります。make_derivativeあなたの答えは 4.5 に近いはずです。

を実行するnewton_search2(gllik(), 0.0001, mean(x), conv = 0.001)と、答えが2倍になります。

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 に答える