Julia でかんたんラプラス近似 - ジョンとヨーコのイマジン日記

データを生成した分布を「真の分布」と呼ぶことにする。

確率モデルを $p(x|w)$ として、確率モデルにパラメータ $w$ の最尤推定量 $\hat w$ を代入し、
「真の分布はおよそ $p(x|\hat w)$ だろう」
と推測することを最尤推定による予測分布（略して、MLE）と呼ぶことにする。

最尤推定量（またはMAP推定量）を中心にした正規分布で事後分布を近似することをラプラス近似という。

ラプラス近似の際、事後分布を近似する正規分布の分散共分散行列は負の対数尤度のヘシアンの逆行列とする。

ラプラス近似した事後分布 $f(w)$ で確率モデルを平均したもの
$E_w [p(x|w)] = \int p(x|w)f(w) \, dw$
を用いて、
「真の分布はおよそ $E_w [p(x|w)]$ だろう」
と推測することをラプラス近似による予測分布（略して、Laplace）と呼ぶことにする。

今回はMLEとLaplaceを、真の分布からのカルバック・ライブラ距離で比較する。

真の分布としてはワイブル分布とガンマ分布を設定して、確率モデルもワイブル分布とガンマ分布とした。

ワイブル分布やガンマ分布はシェイプパラメータが1より大きいときと小さいときで大きく性質が変わるので、シェイプパラメータ2のときと0.5のときの2パターンのシミュレーションを行った。

真の分布のスケールパラメータは1とした。

推定する未知パラメータはシェイプパラメータとスケールパラメータの2個である。

シミュレーションの試行回数は1000回とした。

サンプルサイズ n は 25、50、100とした。

ワイブル分布でシェイプパラメータ2のとき：

f:id:abrahamcow:20201210231329p:plain

ワイブル分布でシェイプパラメータ0.5のとき：

f:id:abrahamcow:20201210231437p:plain

ガンマ分布でシェイプパラメータ2のとき：

f:id:abrahamcow:20201210231515p:plain

ガンマ分布でシェイプパラメータ0.5のとき：

f:id:abrahamcow:20201210231535p:plain

いずれの場合もラプラス近似による予測分布のほうが、最尤推定による予測分布にくらべ、カルバック・ライブラ距離の意味で大外しすることが減っていることがわかる。

その傾向はサンプルサイズが小さいときにより顕著である。

Julia のコード

using Distributions
using Optim
using Random
using StatsFuns
using LinearAlgebra
using QuadGK
using DataFrames
using RCall
using ForwardDiff

logmeanexp(x)=logsumexp(x)-log(length(x))

function ll_wei(par,y)  
    f = Weibull(exp(par[1]),exp(par[2]))
    ll = sum(logpdf(f,y[i]) for i in eachindex(y))
    return -ll
end

function logpred_wei(x,m,η)
    return logmeanexp(logpdf(Weibull(exp(m[i]),exp(η[i])),x) for i in eachindex(m))
end

function simGE_weibull(n,m,iter)
    KL1 = zeros(iter)
    KL2 = zeros(iter)
    rgn = Random.default_rng()
    truedist = Weibull(m,1.0)
    entropy, err = quadgk(x -> -logpdf(truedist,x)*pdf(truedist,x),0,Inf)
    for i in 1:iter
y = rand(rgn,truedist,n)
opt = optimize(par -> ll_wei(par,y),zeros(2),method=NewtonTrustRegion(),autodiff=:forward)
mu = Optim.minimizer(opt)
H = Symmetric(ForwardDiff.hessian(par -> ll_wei(par,y),mu))
Sig = inv(H)
wsamp = rand(rgn,MvNormal(mu,Sig),1000)
ge1, err1 = quadgk(x -> -logpdf(Weibull(exp(mu[1]),exp(mu[2])),x)*pdf(truedist,x),0,Inf)
ge2, err2 = quadgk(x -> -logpred_wei(x,wsamp[1,:],wsamp[2,:])*pdf(truedist,x),0,Inf)
KL1[i] = ge1-entropy
KL2[i] = ge2-entropy
end
KL1, KL2
end

Random.seed!(1)
@time KL1, KL2 = simGE_weibull(25,2,1000)
df1 = DataFrame(n=25,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_weibull(50,2,1000)
df2 = DataFrame(n=50,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_weibull(100,2,1000)
df3 = DataFrame(n=100,MLE=KL1,Laplace=KL2)

dfc_wei2 = stack([df1;df2;df3],[:MLE,:Laplace])
rename!(dfc_wei2,Dict(:variable => "method", :value => "KL"))

R"
library(ggplot2)
ggplot($dfc_wei2,aes(x=factor(n),y=KL,colour=method))+
geom_violin()+
xlab('n')+
theme_bw()
ggsave('wei2.png')
"

@time KL1, KL2 = simGE_weibull(25,0.5,1000)
df1 = DataFrame(n=25,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_weibull(50,0.5,1000)
df2 = DataFrame(n=50,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_weibull(100,0.5,1000)
df3 = DataFrame(n=100,MLE=KL1,Laplace=KL2)

dfc_wei05 = stack([df1;df2;df3],[:MLE,:Laplace])
rename!(dfc_wei05,Dict(:variable => "method", :value => "KL"))

R"
library(ggplot2)
ggplot($dfc_wei05,aes(x=factor(n),y=KL,colour=method))+
geom_violin()+
xlab('n')+
theme_bw()
ggsave('wei05.png')
"

function ll_gam(par,y)
    f = Gamma(exp(par[1]),exp(par[2]))
    ll = sum(logpdf(f,y[i]) for i in eachindex(y))
    return -ll
end

function logpred_gam(x,a,b)
    return logmeanexp(logpdf(Gamma(exp(a[i]),exp(b[i])),x) for i in eachindex(a))
end

function simGE_gamma(n,a,iter)
    KL1 = zeros(iter)
    KL2 = zeros(iter)
    rgn = Random.default_rng()
    truedist = truedist = Gamma(a,1.0)
    entropy, err = quadgk(x -> -logpdf(truedist,x)*pdf(truedist,x),0,Inf)
    for i in 1:iter
y = rand(rgn,truedist,n)
opt = optimize(par -> ll_gam(par,y),zeros(2),method=NewtonTrustRegion(),autodiff=:forward)
mu = Optim.minimizer(opt)
H = Symmetric(ForwardDiff.hessian(par -> ll_gam(par,y),mu))
Sig = inv(H)
wsamp = rand(rgn,MvNormal(mu,Sig),1000)
ge1, err1 = quadgk(x -> -logpdf(Gamma(exp(mu[1]),exp(mu[2])),x)*pdf(truedist,x),0,Inf)
ge2, err2 = quadgk(x -> -logpred_gam(x,wsamp[1,:],wsamp[2,:])*pdf(truedist,x),0,Inf)
KL1[i] = ge1-entropy
KL2[i] = ge2-entropy
end
KL1, KL2
end

Random.seed!(1)
@time KL1, KL2 = simGE_gamma(25,2,1000)
df1 = DataFrame(n=25,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_gamma(50,2,1000)
df2 = DataFrame(n=50,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_gamma(100,2,1000)
#549.924845 seconds (48.92 M allocations: 29.257 GiB, 0.63% gc time)
df3 = DataFrame(n=100,MLE=KL1,Laplace=KL2)

dfc_gam2 = stack([df1;df2;df3],[:MLE,:Laplace])
rename!(dfc_gam2,Dict(:variable => "method", :value => "KL"))


R"
library(ggplot2)
ggplot($dfc_gam2,aes(x=factor(n),y=KL,colour=method))+
geom_violin()+
xlab('n')+
theme_bw()
ggsave('gam2.png')
"

@time KL1, KL2 = simGE_gamma(25,0.5,1000)
df1 = DataFrame(n=25,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_gamma(50,0.5,1000)
df2 = DataFrame(n=50,MLE=KL1,Laplace=KL2)
@time KL1, KL2 = simGE_gamma(100,0.5,1000)
df3 = DataFrame(n=100,MLE=KL1,Laplace=KL2)

dfc_gam05 = stack([df1;df2;df3],[:MLE,:Laplace])
rename!(dfc_gam05,Dict(:variable => "method", :value => "KL"))

R"
library(ggplot2)
ggplot($dfc_gam05,aes(x=factor(n),y=KL,colour=method))+
geom_violin()+
xlab('n')+
theme_bw()
ggsave('gam05.png')
"

いくつかメモを。

最尤推定（尤度関数の最大化）を行っているのはこの一行です。

opt = optimize(par -> ll_wei(par,y),zeros(2),method=NewtonTrustRegion(),autodiff=:forward)

autodiff=:forward を指定することにより、自動微分で正確な微分を行います。

対数尤度のヘシアンを求めるのはこの一行です。

H = Symmetric(ForwardDiff.hessian(par -> ll_wei(par,y),mu))

Symmetricという関数に入れることにより、対称行列であることが明示され、逆行列を計算したときの数値エラーが出にくくなります。

dfc_wei2 = stack([df1;df2;df3],[:MLE,:Laplace])

stack は R でいう tidyr のgather みたいな関数です。

バイオリンプロットの描画はRCallでRをよんでいます。