高速なガンマ分布の最尤推定法の補足

# ディガンマ関数の逆関数 igamma() を使うため
library(distr)

### ガンマ関数のパラメータを最尤推定する関数 ###
estimate_gamma <- function(values, method = c("FOApprox", "LocalApprox")) {
  method <- match.arg(method)

  # 固定値をあらかじめ計算しておく
  mean_values <- mean(values)
  mean_log_values <- mean(log(values))
  log_mean_values <- log(mean_values)

  # 対数尤度を計算する関数
  compute_loglik <- function(values, a, b) sum(dgamma(values, a, b, log = TRUE))

  # パラメータの初期値を設定
  a <- 0.5 / (log_mean_values - mean_log_values)
  b <- mean_values / a
  # 対数尤度の初期値を計算
  log_lik_old <- compute_loglik(values, a, b)

  # 対数尤度の変化を格納する変数
  epsilon <- Inf 
  # 対数尤度の変化が小さくなるまでパラメータを更新
  while (epsilon > 1e-5) {
    if (method == "FOApprox") {
      # 式 (5) によりパラメータ a を更新
      a <- igamma(mean_log_values - log_mean_values + log(a))
    } else {
      # 式 (6) によりパラメータ a を更新
      inv_a <- 1/a + (mean_log_values - digamma(a) - log_mean_values + log(a)) / (a^2 * (1/a - trigamma(a)))
      a <- 1 / inv_a
    }
    # 式 (2) によりパラメータ b を更新
    b <- mean_values / a

    # 対数尤度を計算
    log_lik <- compute_loglik(values, a, b)
    # 対数尤度の変化を計算
    epsilon <- abs(log_lik - log_lik_old)
    # 現在の対数尤度を保管
    log_lik_old <- log_lik
  }
  c(a = a, b = b)
}

set.seed(314)
# ガンマ分布 Gamma(3, 2) からサンプリング
values <- rgamma(n = 500, shape = 3, scale = 2)

# 一次近似下限最大化法による推定
estimate_gamma(values, method = "FOApprox")

# 局所近似最大化法による推定
estimate_gamma(values, method = "LocalApprox")