Rでラプラス近似

仕組み

\(f(\theta)\) をモード周りで2次までテイラー展開すると、モードでは1次微分がゼロであるため

\[f(\theta) \approx f(\hat\theta) + \frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2\]

となり、正規分布の対数密度と同じ形になります。

これによって、正規化された確率密度関数 \(p(\theta)\) は次のような正規分布で近似できます。

\[p(\theta) \approx N\left(\hat\theta,\ -\frac{1}{f''(\hat\theta)}\right)\]

つまり「モードの位置」と「モードでの曲率(2階微分)」の2つだけから、分布全体を正規分布で近似できます。

特徴

• 必要なのは モード探索(最適化) と 2階微分(Hessian) だけですので、MCMCなどのサンプリング手法より計算コストが低くなります。
• サンプルサイズが大きいほど事後分布は正規分布に近づくため(漸近正規性)、近似精度が上がります。
• 分布が非対称(歪んでいる)の場合や多峰性の場合は精度が落ちます。
• 境界のあるパラメータ(0〜1など)では、logitなど無制約な空間に変換してから近似すると改善することがあります。

導出

以下はモード周りでのテイラー展開の導出です。

1. 一般のテイラー展開

任意の点 \(\theta\) を、モード \(\hat\theta\) の近くの点として、\(f\) を \(\hat\theta\) の周りで展開します。テイラーの定理より

\[f(\theta) = f(\hat\theta) + f'(\hat\theta)(\theta-\hat\theta) + \frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2 + \frac{1}{6}f'''(\hat\theta)(\theta-\hat\theta)^3 + \cdots\]

これは単に \(f\) を \(\hat\theta\) での微分係数を使って多項式で書き直しただけです。

2. 2次の項までで打ち切る

\(\theta\) が \(\hat\theta\) に十分近い範囲では \((\theta-\hat\theta)^3\) 以降の項は急速に小さくなるとみなし、無視します。

\[f(\theta) \approx f(\hat\theta) + f'(\hat\theta)(\theta-\hat\theta) + \frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2\]

3. モードであることを使う

\(\hat\theta\) は \(f\) の最大値を与える点(モード)ですので、定義より

\[f'(\hat\theta) = 0\]

が成り立ちます（内部の最大点では微分がゼロになるため）。これを代入すると、1次の項が消えて

\[f(\theta) \approx f(\hat\theta) + \frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2\]

だけが残ります。

4. 指数を取って正規分布の形に対応づける

\(p(\theta) \propto \exp(f(\theta))\) ですので、両辺の指数を取ると

\[p(\theta) \propto \exp\left(f(\hat\theta) + \frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2\right) = \exp(f(\hat\theta)) \cdot \exp\left(\frac{1}{2}f''(\hat\theta)(\theta-\hat\theta)^2\right)\]

\(\exp(f(\hat\theta))\) は \(\theta\) に依存しない定数ですので、\(\theta\) に依存する部分だけを見ると

\[p(\theta) \propto \exp\left(-\frac{1}{2}\cdot\frac{(\theta-\hat\theta)^2}{-1/f''(\hat\theta)}\right)\]

これはまさに正規分布 \(N(\hat\theta,\ \sigma^2)\) の密度のカーネル \(\exp\left(-\frac{(\theta-\mu)^2}{2\sigma^2}\right)\) と同じ形です。両者を見比べると

\[\sigma^2 = -\frac{1}{f''(\hat\theta)}\]

が対応することがわかります。\(\hat\theta\) はモード（\(f\) の最大点）なので \(f''(\hat\theta) < 0\) であり、\(\sigma^2 > 0\) となって分散として矛盾しない点も確認できます。

以上より、

\[p(\theta) \approx N\left(\hat\theta,\ -\frac{1}{f''(\hat\theta)}\right)\]

が導かれます。「1次の項がモードの定義によりゼロになる」という点が、この導出でいちばん重要なステップです。

具体例の設定: 二項分布の尤度 + ベータ事前分布

• データ
  • \(n\)回の試行中 \(y\)回成功
• モデル
  • \(y \mid \theta \sim \text{Binomial}(n, \theta)\)
• 事前分布
  • \(\theta \sim \text{Beta}(a, b)\)

このとき事後分布は解析的に \(\text{Beta}(y+a,\ n-y+b)\) になることが知られています。

つまり「真の答え」がわかっているため、ラプラス近似がどれくらい正確であるか検証することができます。

あえて 1 に近い側へ非対称(歪んだ)な事後分布になるデータを選び、正規近似の弱点が見えるようにします。

n <- 20
y <- 18
a <- 2
b <- 2

cat(sprintf("データ: n=%d, y=%d,  事前分布: Beta(%d,%d)\n", n, y, a, b))
cat(sprintf("真の事後分布: Beta(%d, %d)\n\n", y + a, n - y + b))

# 対数事後密度 (正規化定数を除いたもの) f(theta)
#   p(theta|y) ∝ theta^(y+a-1) * (1-theta)^(n-y+b-1)
#   f(theta) = (y+a-1)*log(theta) + (n-y+b-1)*log(1-theta)
f_theta <- function(theta) {
  (y + a - 1) * log(theta) + (n - y + b - 1) * log(1 - theta)
}

データ: n=20, y=18,  事前分布: Beta(2,2)
真の事後分布: Beta(20, 4)

数値的な2階微分(曲率)を求める関数

中心差分による2階微分の近似 : \[f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}\]

second_derivative <- function(f, x, h = 1e-4) {
  (f(x + h) - 2 * f(x) + f(x - h)) / h^2
}

パターンA: theta の空間で直接ラプラス近似する

\(\theta \in (0,1)\)ですので、モード探索は1次元の有界最適化 optimize() を使います。

theta_star <- optimize(f_theta, interval = c(1e-8, 1 - 1e-8), maximum = TRUE)$maximum
curvature_theta <- second_derivative(f_theta, theta_star)
var_theta <- -1 / curvature_theta

cat("--- パターンA: theta空間で直接ラプラス近似 ---\n")
cat(sprintf("モード theta_hat        = %.5f\n", theta_star))
cat(sprintf(
  "解析解 (y+a-1)/(n+a+b-2) = %.5f  (検算用)\n",
  (y + a - 1) / (n + a + b - 2)
))
cat(sprintf("近似分散 -1/f''(theta_hat) = %.6f\n\n", var_theta))

laplace_theta_direct <- function(theta) {
  dnorm(theta, mean = theta_star, sd = sqrt(var_theta))
}

--- パターンA: theta空間で直接ラプラス近似 ---
モード theta_hat        = 0.86364
解析解 (y+a-1)/(n+a+b-2) = 0.86364  (検算用)
近似分散 -1/f''(theta_hat) = 0.005353

パターンB: ロジット変換した空間でラプラス近似する

\(\theta\) は \((0,1)\) に制約されているため、正規分布で直接近似すると0や1を跨いだ「あり得ない」領域にまで裾が広がってしまいます。

そこで

\[\psi = \operatorname{logit}(\theta) = \log\left(\frac{\theta}{1-\theta}\right)\]

という無制約の変数に変換してからラプラス近似をかけ、その後 \(\theta\) に戻す方が精度が良くなることがあります。

変数変換をする場合、確率密度にはヤコビアン(変換の微分)を掛ける必要があります:

\[p(\psi) = p(\theta(\psi)) \cdot \left| \frac{d\theta}{d\psi} \right|\]

\(\theta = \frac{1}{1 + \exp(-\psi)}\) （ロジスティック関数）のとき、その微分は \(\frac{d\theta}{d\psi} = \theta(1-\theta)\) となります。

これにより、\(\psi\) 空間における新たな対数(未正規化)密度 \(g(\psi)\) は、元の対数未正規化密度 \(f(\theta)\) を用いて以下のように表せます。 \[g(\psi) = f(\theta(\psi)) + \log(\theta(\psi)) + \log(1 - \theta(\psi))\]

g_psi <- function(psi) {
  th <- plogis(psi)
  f_theta(th) + log(th) + log(1 - th)
}

psi_star <- optimize(g_psi, interval = c(-15, 15), maximum = TRUE)$maximum
curvature_psi <- second_derivative(g_psi, psi_star)
var_psi <- -1 / curvature_psi

cat("--- パターンB: logit(theta)空間でラプラス近似してtheta空間へ戻す ---\n")
cat(sprintf("モード psi_hat = logit(theta_hat) = %.5f\n", psi_star))
cat(sprintf("対応する theta = plogis(psi_hat)  = %.5f\n", plogis(psi_star)))
cat(sprintf("近似分散 -1/g''(psi_hat)          = %.6f\n\n", var_psi))

# psi空間の正規分布を、変数変換(ヤコビアン |dpsi/dtheta| = 1/(theta(1-theta)))
# を通じて theta 空間の密度に戻す
laplace_theta_via_psi <- function(theta) {
  psi <- qlogis(theta)
  jac <- 1 / (theta * (1 - theta)) # dpsi/dtheta
  dnorm(psi, mean = psi_star, sd = sqrt(var_psi)) * jac
}

--- パターンB: logit(theta)空間でラプラス近似してtheta空間へ戻す ---
モード psi_hat = logit(theta_hat) = 1.60944
対応する theta = plogis(psi_hat)  = 0.83333
近似分散 -1/g''(psi_hat)          = 0.300000

真の事後分布と2種類のラプラス近似を重ねて可視化

library(ggplot2)

theta_grid <- seq(0.001, 0.999, length.out = 500)

df_compare <- data.frame(
  theta = rep(theta_grid, 3),
  density = c(
    dbeta(theta_grid, y + a, n - y + b),
    laplace_theta_direct(theta_grid),
    laplace_theta_via_psi(theta_grid)
  ),
  type = rep(
    c(
      "真の事後分布 Beta(y+a, n-y+b)",
      "ラプラス近似 (theta空間で直接)",
      "ラプラス近似 (logit変換してから)"
    ),
    each = length(theta_grid)
  )
)

p1 <- ggplot(df_compare, aes(x = theta, y = density, color = type, linetype = type)) +
  geom_line(linewidth = 1) +
  labs(
    title = sprintf("ラプラス近似の比較 (n=%d, y=%d, Beta(%d,%d)事前分布)", n, y, a, b),
    x = expression(theta), y = "density", color = NULL, linetype = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p1)

真の事後分布が非対称(歪んでいる)の場合、logit変換してから近似する方が精度が良い事を確認できます。

ラプラス法による「積分(正規化定数)」の近似

真の値(ベータ関数)

\[\int_0^1 \theta^{y+a-1} (1-\theta)^{n-y+b-1} d\theta = B(y+a,\ n-y+b)\]

正しく変数変換を行えば真の積分値は不変ですが、ラプラス近似によって得られる近似値は、どの空間で近似を施したかによって変わります。

一般に、事後分布がより正規分布に近い形になる変数（この場合は logit スケール）の空間で近似した方が、真の積分値に近くなります。

true_integral <- beta(y + a, n - y + b)

# theta空間で直接ラプラス近似した積分
laplace_integral_theta <- exp(f_theta(theta_star)) * sqrt(2 * pi * var_theta)

# psi空間でラプラス近似した積分
laplace_integral_psi <- exp(g_psi(psi_star)) * sqrt(2 * pi * var_psi)

cat("--- ラプラス法による正規化定数(積分)の近似 ---\n")
cat(sprintf("真の値                 B(y+a, n-y+b)      = %.8f\n", true_integral))
cat(sprintf(
  "ラプラス近似 (theta空間) exp(f)*sqrt(2*pi*var) = %.8f  (誤差 %.3f%%)\n",
  laplace_integral_theta,
  100 * (laplace_integral_theta - true_integral) / true_integral
))
cat(sprintf(
  "ラプラス近似 (psi空間)   exp(g)*sqrt(2*pi*var) = %.8f  (誤差 %.3f%%)\n\n",
  laplace_integral_psi,
  100 * (laplace_integral_psi - true_integral) / true_integral
))

--- ラプラス法による正規化定数(積分)の近似 ---
真の値                 B(y+a, n-y+b)      = 0.00002823
ラプラス近似 (theta空間) exp(f)*sqrt(2*pi*var) = 0.00002869  (誤差 1.626%)
ラプラス近似 (psi空間)   exp(g)*sqrt(2*pi*var) = 0.00002763  (誤差 -2.126%)

ラプラス近似は変数変換に対して不変ではなく、パラメータ化依存性がある事を確認できます。

サンプルサイズを増やすとラプラス近似がどう改善するか

事後分布の漸近正規性(ベルンシュタイン・フォンミーゼスの定理)により、データ数 n が大きくなるほど事後分布は正規分布に近づき、ラプラス近似の精度も上がっていくことを確認します。

p_true <- 0.9 # 成功確率の目安 (y/nがこの値に近くなるように設定)
n_values <- c(20, 50, 200, 1000)

build_df_for_n <- function(n_i) {
  y_i <- round(n_i * p_true)
  f_i <- function(theta) (y_i + a - 1) * log(theta) + (n_i - y_i + b - 1) * log(1 - theta)
  theta_hat_i <- optimize(f_i, interval = c(1e-8, 1 - 1e-8), maximum = TRUE)$maximum
  var_i <- -1 / second_derivative(f_i, theta_hat_i)

  data.frame(
    theta = rep(theta_grid, 2),
    density = c(
      dbeta(theta_grid, y_i + a, n_i - y_i + b),
      dnorm(theta_grid, theta_hat_i, sqrt(var_i))
    ),
    type = rep(c("真の事後分布", "ラプラス近似"), each = length(theta_grid)),
    n = sprintf("n = %d", n_i)
  )
}

df_n <- do.call(rbind, lapply(n_values, build_df_for_n))
df_n$n <- factor(df_n$n, levels = sprintf("n = %d", n_values))

p2 <- ggplot(df_n, aes(x = theta, y = density, color = type, linetype = type)) +
  geom_line(linewidth = 1) +
  facet_wrap(~n, scales = "free") +
  labs(
    title = "サンプルサイズ n が増えるほどラプラス近似は真の分布に近づく",
    subtitle = sprintf("成功確率の目安 p ≈ %.1f, 事前分布 Beta(%d,%d) は固定", p_true, a, b),
    x = expression(theta), y = "density", color = NULL, linetype = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p2)

LaplaceApproximation {LaplacesDemon} との比較

library(LaplacesDemon)

Data <- list(
  y = y, n = n, a = a, b = b,
  mon.names = "theta", # 監視する量の名前(元スケールのthetaを見る)
  parm.names = "psi" # 最適化するパラメータの名前(無制約のpsi)
)

Model <- function(parm, Data) {
  psi <- parm[1]
  theta <- plogis(psi)
  LP <- (Data$y + Data$a - 1) * log(theta) +
    (Data$n - Data$y + Data$b - 1) * log(1 - theta) +
    log(theta) + log(1 - theta) # 事後密度 + ヤコビアン
  Dev <- -2 * LP
  Modelout <- list(
    LP = LP, Dev = Dev, Monitor = theta,
    yhat = theta, parm = parm
  )
  return(Modelout)
}

Initial.Values <- c(psi = 0) # 最適化の初期値 (theta=0.5に相当)

Fit <- LaplaceApproximation(
  Model,
  parm = Initial.Values, Data = Data,
  Iterations = 200, Method = "BFGS", CovEst = "Hessian", sir = FALSE
)

psi_star_pkg <- Fit$Summary1["psi", "Mode"]
var_psi_pkg <- Fit$Covar[1, 1]

cat("--- LaplacesDemon::LaplaceApproximation() による結果 ---\n")
cat(sprintf("収束したか: %s\n", Fit$Converged))
cat(sprintf(
  "モード psi_hat = %.6f   (theta換算 = %.6f)\n",
  psi_star_pkg, plogis(psi_star_pkg)
))
cat(sprintf(
  "分散  var_psi  = %.6f   (SD = %.6f)\n\n",
  var_psi_pkg, sqrt(var_psi_pkg)
))

Sample Size:  1 
Laplace Approximation begins...
Iteration:  20  of  200 ,   LP:  -10.8 
Estimating the Covariance Matrix
Creating Summary from Point-Estimates
Estimating Log of the Marginal Likelihood
Laplace Approximation is finished.

--- LaplacesDemon::LaplaceApproximation() による結果 ---
収束したか: TRUE
モード psi_hat = 1.609411   (theta換算 = 0.833330)
分散  var_psi  = 0.299996   (SD = 0.547719)

comparison <- data.frame(
  項目 = c("psi_hat (モード)", "theta_hat (モード, 換算)", "var_psi (分散)"),
  自作実装 = c(psi_star, plogis(psi_star), var_psi),
  パッケージ = c(psi_star_pkg, plogis(psi_star_pkg), var_psi_pkg)
)
comparison$差の絶対値 <- abs(comparison$自作実装 - comparison$パッケージ)

cat("--- 比較表 ---\n")
print(comparison, digits = 8)

--- 比較表 ---
                      項目   自作実装 パッケージ    差の絶対値
1         psi_hat (モード) 1.60943575 1.60941087 2.4878711e-05
2 theta_hat (モード, 換算) 0.83333303 0.83332958 3.4554101e-06
3           var_psi (分散) 0.29999954 0.29999582 3.7249065e-06

\(\theta\) 空間での密度を重ねて最終確認

theta_grid <- seq(0.001, 0.999, length.out = 500)

density_manual <- function(theta) {
  psi <- qlogis(theta)
  jac <- 1 / (theta * (1 - theta))
  dnorm(psi, psi_star, sqrt(var_psi)) * jac
}

density_pkg <- function(theta) {
  psi <- qlogis(theta)
  jac <- 1 / (theta * (1 - theta))
  dnorm(psi, psi_star_pkg, sqrt(var_psi_pkg)) * jac
}

df_final <- data.frame(
  theta = rep(theta_grid, 3),
  density = c(
    dbeta(theta_grid, y + a, n - y + b),
    density_manual(theta_grid),
    density_pkg(theta_grid)
  ),
  type = rep(c("真の事後分布", "自作のラプラス近似", "LaplacesDemonパッケージ"),
    each = length(theta_grid)
  )
)

p_final <- ggplot(df_final, aes(x = theta, y = density, color = type, linetype = type)) +
  geom_line(linewidth = 1) +
  labs(
    title = "自作実装とLaplacesDemonパッケージのラプラス近似の比較",
    subtitle = sprintf("n=%d, y=%d, Beta(%d,%d)事前分布 (logit変換空間で近似)", n, y, a, b),
    x = expression(theta), y = "density", color = NULL, linetype = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(p_final)

以上です。