不完全ベータ関数の近似はまとめるほど理解できなかったので、また時間があるときにやろう。
今日はラプラス近似が何をしてるかを少しグラフィカルに検討した。
\[
\int_{\mathbf b\in R^q}\exp\{h(\mathbf b)\}d\mathbf b \simeq \exp\{h(\hat{\mathbf b})\} \times \int_{\mathbf b\in R^q}\exp\left\{ \frac{1}{2}(\mathbf b - \hat{\mathbf b})^t h^{\prime\prime}(\hat{\mathbf b})(\mathbf b - \hat{\mathbf b})\right\} d\mathbf b
\]
を元に考えてみる。
シンプルに $\exp\{h(\mathbf b)\}$ を変量効果の分布が一様分布な全尤度関数 $p(\mathbf Y=\mathbf y \mid \mathbf b)$ と考えた場合、$h(\mathbf b)=\log \{ p(\mathbf Y=\mathbf y \mid \mathbf b) \}$ になる。
上式の積分を取っ払ってみると
\[
p(\mathbf Y = \mathbf y \mid \mathbf b) \simeq
p(\mathbf Y = \mathbf y \mid \hat{\mathbf b})\times \exp\left\{ \frac{1}{2}(\mathbf b - \hat{\mathbf b})^t \frac{d^2}{d\mathbf b d\mathbf b^t}\log\{p(\mathbf Y = \mathbf y \mid \mathbf b)\}\mid_{\mathbf b=\hat{\mathbf b}}(\mathbf b - \hat{\mathbf b})\right\}
\] と書ける。
式の形状からして、尤度関数を正規分布の形で近似している事が分かる。
従って、ラプラス近似をする場合は、尤度関数が正規分布に近いほど近似が良くなるという事になる (尤度関数の定義域が違う場合も当然誤差が増える)。
$p(\mathbf Y=\mathbf y \mid \hat{\mathbf b})$ は、$p(\mathbf Y=\mathbf y \mid \mathbf b)$ の最大値であり、$\hat{\mathbf b}$ は最大値を与える引数だったので、最大値の部分はぴったり一致する事も分かる。
分散に当たる項も Fisher 情報行列にマイナスをつけて、$\hat{\mathbf b}$ をぶち込んだような形になっている (期待値を取ってるわけではないので同じではない)。
二項分布の尤度関数を近似してみる。
実際の積分値 は
\[
I(n, y) = 1/(n+1)
\] ラプラス近似は
\[
\tilde{I}(n, y)=\frac{\sqrt{2\pi}{n \choose y}(y/n)^y(1-y/n)^{n-y}}{\sqrt{\frac{n^3}{y(n-y)}}}
\] になる。
$n$ が大きく $0.5$ に近い所を中心に分布している場合は、正規分布によく似た形をしている。
こういう場合は近似がよいはず。
\[
I(500, 240)=1.996\times 10^{-3}, \quad \tilde{I}(500, 240)=1.999\times 10^{-3}
\]
$n$ が小さく端っこに分布している場合で、形がかなり異なる。
こういう場合は近似がよくないはず。
\[
I(30, 2)=3.226\times 10^{-2}, \quad \tilde{I}(30, 2)=3.198\times 10^{-2}
\]
########################################## # 近似がいい場合 n <- 500 y <- 240 Bin <- function(x) { exp(lchoose(n, y) + y * log(x) + (n-y) * log(1-x)) } LAP <- function(x) { d2lf <- function(n, y) { -exp((3 * log(n) - (log(n - y) + log(y)))) } lap <- Bin(y/n) * exp((x - y/n) * d2lf(n, y) * (x - y/n) / 2) } postscript( file=paste("c:/lap_1.eps"), width = 9, height = 18, pointsize = 20, horizontal = FALSE, onefile = FALSE, paper = "special", family="Palatino" ) par(mar=c(4.6, 4.6, 2.1, 1.6), pch=1) par(mfrow=c(2, 1)) curve( Bin, 0.35, 0.6, n = 1000, lty = 1, xlab = expression(paste(theta)), ylab = "binomial likelihood", main = "n = 500, y = 240" ) curve( LAP, 0.35, 0.6, n = 1000, lty = 2, xlab = expression(paste(theta)), ylab = "laplace appriximation likelihood" ) dev.off() ########################################## # 近似が良くない場合 n <- 30 y <- 2 Bin <- function(x) { exp(lchoose(n, y) + y * log(x) + (n-y) * log(1-x)) } LAP <- function(x) { d2lf <- function(n, y) { -exp((3 * log(n) - (log(n - y) + log(y)))) } lap <- Bin(y/n) * exp((x - y/n) * d2lf(n, y) * (x - y/n) / 2) } postscript( file=paste("c:/lap_2.eps"), width = 9, height = 9, pointsize = 20, horizontal = FALSE, onefile = FALSE, paper = "special", family="Palatino" ) par(mar=c(4.6, 4.6, 2.1, 1.6), pch=1) par(mfrow=c(1, 1)) curve( Bin, 0, 1, n = 1000, lty = 1, xlab = expression(paste(theta)), ylab = "binomial & laplace appriximation likelihood", main = "n = 30, y = 2" ) curve( LAP, 0, 1, n = 1000, lty = 2, add = TRUE, xlab = expression(paste(theta)), ylab = "laplace appriximation likelihood" ) legend( 0.77, 0.28, c("Bin", "LAP"), lty = c(1, 2), merge = TRUE ) dev.off()
R で描いたら結構綺麗だった。
