黒木玄
2018-04-02, 2018-12-20
http://nbviewer.jupyter.org/gist/genkuroki/8a342d0b7b249e279dd8ad6ae283c5db に書いた公式のSymPyによる導出.
using Base.Meta: parse
using BenchmarkTools
using PyPlot
using SymPy
using SpecialFunctions
using Statistics
const latex = sympy[:latex]
using LaTeXStrings
latexstring_displaystyle(args...; kwargs...) =
LaTeXString(raw"$$" * prod(latex.(args; kwargs...)) * raw"$$")
latexdisp(args...; kwargs...) =
display(latexstring_displaystyle(args...; kwargs...))
const ls = latexstring_displaystyle
const ld = latexdisp
latexdisp (generic function with 1 method)
@syms a b c d
@syms x y μ μ₀ real=true
Y = [symbols("Y$i", real=true) for i in 1:5]
@syms σ λ λ₀ θ α β ξ η positive=true
(σ, λ, λ₀, θ, α, β, ξ, η)
Y
p_N = 1/√(2PI)*exp(-λ*(y-μ)^2/2)*√λ
p_μ = p_N(μ=>μ₀, y=>μ, λ=>λ*λ₀)
p_λ = exp(-λ/θ)*λ^(α-1)/(gamma(α)*θ^α)
WAICやWBICを計算するためには対数尤度函数とその二乗の事後分布に関する期待値を計算しなければいけなくなる. 共役事前分布を採用した場合には事後分布も共役事前分布と同じ形になるので、共役事前分布に関する期待値の公式を作っておけば十分である. 以下ではその計算をSymPyを使って行う.
全体を $\mu_0$ だけ平行移動すればいつでも $\mu_0=0$ となるようにできる. 以下ではSymPyでの積分計算を易しくするために $\mu_0=0$ と仮定する.
# for simplicity, put μ₀ = 0.
p_μ = p_μ(μ₀=>0)
logp_N = simplify(log(p_N))
logp_μ = simplify(log(p_μ))
logp_λ = expand(log(p_λ))
p_N = simplify(p_N)
p_μ = simplify(p_μ)
p_λ = simplify(p_λ)
logp_N = expand(logp_N)
logp_μ = expand(logp_μ)
logp_λ = expand(logp_λ)
[p_N, p_μ, p_λ, logp_N, logp_μ, logp_λ]
@time I_μ = simplify(integrate(expand(logp_N*p_μ), (μ, -oo, oo)))
0.667942 seconds (79.46 k allocations: 3.975 MiB)
# 対数尤度函数の期待値
@time Elogp_N = simplify(integrate(expand(I_μ*p_λ), (λ, 0, oo)))
1.858062 seconds (80.13 k allocations: 4.003 MiB)
A1 = -2Elogp_N
A2 = log(2PI) + α*θ*y^2 + 1/λ₀ - (polygamma(Sym(0),α) + log(θ))
simplify(A1 - A2)
expand(logp_N^2)
@time J_μ = integrate(expand(logp_N^2*p_μ), (μ, -oo, oo))
1.075623 seconds (36 allocations: 896 bytes)
# 対数尤度函数の二乗の期待値
@time Elogp_N2 = simplify(integrate(expand(J_μ*p_λ), (λ, 0, oo)))
8.517231 seconds (24 allocations: 560 bytes)
@syms ψ₀ ψ₁
B1 = 4Elogp_N2(polygamma(Sym(0),α)=>ψ₀, polygamma(Sym(1),α)=>ψ₁)
B2 = (
(log(2PI))^2
+ (α+1)*α*θ^2*y^4 + 6α*θ/λ₀*y^2 + 3/λ₀^2
+ ψ₁ + (ψ₀ + log(θ))^2
+ 2log(2PI)*(α*θ*y^2 + 1/λ₀)
- 2log(2PI)*(ψ₀ + log(θ))
- 2( y^2*(θ + α*θ*(ψ₀ + log(θ))) + (ψ₀+log(θ))/λ₀)
)
simplify(B1 - B2)
サンプルサイズ $n=5$ の場合の共役事前分布のベイズ公式の公式を求めよう. サンプルによって共役事前分布の4つのパラメーター $\mu_0, \lambda_0, \theta, \alpha$ が更新される.
この節では $\mu_0$ を復活させる.
p_μ = p_μ(μ=>μ-μ₀)
n = length(Y)
Ybar = mean(Y)
Yvar = sum((Y.-Ybar).^2)/length(Y)
# 共役事前分布
prior = simplify(p_μ*p_λ)
p5_N = simplify(prod(p_N(y=>v)^β for v in Y))
# 次の p5 を分配函数で割ったものが事後分布になる.
p5 = simplify(p5_N*p_μ*p_λ)
# 以下の式を μ, λ に関する定数倍を除いて上の共役事前分布と等しくするためには,
# 共役事前分布の4つのパラメーターをどのように変えなければいけないかを以下で計算する.
# n=5 の場合の対数尤度函数
logp5 = expand(log(p5))
# exp() の中の式における -λ の係数
c_λ = -coeff(collect(logp5, λ), λ)
@syms a b c
s = solve(c_λ - a/2*(μ-b)^2 - c, [a,b,c])[1]
ld("a = ", s[a])
ld("b = ", s[b])
ld("c = ", s[c])
# μ に関する平方完成
simplify(c_λ - s[a]/2*(μ-s[b])^2 - s[c])
# α₅ は対数尤度函数における log(λ) の係数に 1/2 を足したもの
α₅ = coeff(collect(logp5, log(λ)), log(λ)) + Sym(1)/2
λ₀₅ = s[a]
μ₀₅ = s[b]
θ₅ = 1/simplify(s[c])
α₅ - (α + n*β/2)
λ₀₅ - (λ₀+n*β)
simplify(μ₀₅ - (λ₀*μ₀+n*β*Ybar)/(λ₀+n*β))
tmp = (1/θ)*(1 + (n*β*θ/2)*(λ₀/(λ₀+n*β)*(Ybar - μ₀)^2 + Yvar))
simplify(1/θ₅ - tmp)
まとめ: 以上は $n=5$ の場合の公式. 一般の場合の正規分布モデルの事前共役分布に関するベイズ更新の公式は以下の通り:
$$ \begin{aligned} & \tilde\mu_0 = \frac{\lambda_0\mu_0 + \beta n \bar Y}{\lambda_0 + \beta n}, %\quad & & \tilde\lambda_0 = \lambda_0 + \beta n, \quad \\ & \tilde\alpha = \alpha + \frac{1}{2}\beta n, %\quad & & {\tilde\theta}^{-1} = \theta^{-1}\left[ 1 + \frac{\beta n\theta}{2}\left( \frac{\lambda_0}{\lambda_0 + \beta n}\left(\bar Y - \mu_0\right)^2 + V(Y) \right) \right]. \end{aligned} $$ここでサイズ $n$ のサンプルを $Y_1,\ldots,Y_n$ と書くと,
$$ \bar Y = \frac{1}{2}\sum_{i=1}^n Y_i, \qquad V(Y) = \frac{1}{2}\sum_{i=1}^n (Y_i - \bar Y)^2. $$ゆえに $\beta\to\infty$ で
$$ \tilde\mu_0 = \frac{\lambda_0\mu_0 + \beta n \bar Y}{\lambda_0 + \beta n} \to \bar Y, \qquad \tilde\alpha \tilde\theta = \frac{(\alpha + \frac{1}{2}\beta n)\theta} { \left[ 1 + \frac{\beta n\theta}{2}\left( \frac{\lambda_0}{\lambda_0 + \beta n}\left(\bar Y - \mu_0\right)^2 + V(Y) \right) \right] }\to\frac{1}{V(Y)}. $$このことから $\beta\to\infty$ の極限でベイズ更新は最尤法に一致していることもわかる.
注意: ガンマ分布 $p_\lambda(\lambda|\alpha,\theta)$ の平均と分散はそれぞれ $\alpha\theta$, $\alpha\theta^2$ である. 平均 $\alpha\theta$ を固定したままで分散 $\alpha\theta^2=(\alpha\theta)^2/\alpha$ を小さくするためには $\alpha$ を大きくすればよい. そのベイズ更新結果の平均と分散は $\beta\to\infty$ でそれぞれ $1/V(Y)$ と $0$ に近付く.
正規分布 $p_\mu(\mu|\mu_0,\lambda\lambda_0)$ ($\lambda\lambda_0$ は分散の逆数) の分散を小さくするためには $\lambda\lambda_0$ を大きくすればよい. そのベイズ公式結果の平均と分散はそれぞれ $\bar Y$ と $0$ に近付く.
薄く広がった共役事前分布が欲しければ, $\alpha\theta$ と $\mu_0$ を固定したままで $\alpha$ と $\lambda_0$ を小さくすればよい. 例えば, $\alpha\theta=1$, $\mu_0=0$ と固定したとき, $\alpha$, $\lambda_0$ を小さくして $\theta=1/\alpha$ と設定すれば薄く広がった事前分布が得られる.
prior_param = (μ₀=>0.0, θ=>1/α, λ₀=>0.001, α=>0.000001)
posterior_param = (μ₀=>μ₀₅, λ₀=>λ₀₅, α=>α₅, θ=>θ₅)
data_sample = (Y[1]=>0.7, Y[2]=>1.1, Y[3]=>1.2, Y[4]=>1.2, Y[5]=>1.5)
(Y1 => 0.7, Y2 => 1.1, Y3 => 1.2, Y4 => 1.2, Y5 => 1.5)
prior = (p_μ*p_λ)(prior_param...)
posterior = (p_μ*p_λ)(posterior_param...)(β=>1.0, prior_param..., data_sample...)
posterior = posterior(β=>1.0, prior_param..., data_sample...)
posterior = simplify(posterior)
# lambdify版
f_prior = lambdify(prior, [μ, λ])
f_posterior = lambdify(posterior, [μ, λ])
μs = -1:0.03:2
λs = 10 .^(-1:0.04:3)
z_prior = f_prior.(μs', λs)
z_posterior = f_posterior.(μs', λs)
sleep(0.1)
figure(figsize=(10,4))
subplot(121)
pcolormesh(μs, λs, z_prior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("prior")
subplot(122)
pcolormesh(μs, λs, z_posterior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("posterior")
tight_layout()
# eval parse 版
eval(parse("g_prior(μ,λ) = $prior"))
eval(parse("g_posterior(μ,λ) = $posterior"))
μs = -1:0.03:2
λs = 10 .^(-1:0.04:3)
z_prior = g_prior.(μs', λs)
z_posterior = g_posterior.(μs', λs)
sleep(0.1)
figure(figsize=(10,4))
subplot(121)
pcolormesh(μs, λs, z_prior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("prior")
subplot(122)
pcolormesh(μs, λs, z_posterior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("posterior")
tight_layout()
# SymPy expression
posterior
# String version of SymPy expression
"$posterior"
"0.00953627593040395*sqrt(2)*λ^2.000001*exp(-λ*(2.5005*(μ - 0.56994300569943)^2 + 0.166650670065987))/sqrt(pi)"
# 101×101 ≈ 10^4 (μ, λ)'s
μs = collect(-1:0.03:2)
λs = 10 .^(-1:0.04:3)
@show length.((μs,λs));
length.((μs, λs)) = (101, 101)
# SymPy lambdify 版
f_posterior = lambdify(posterior, [μ, λ])
@benchmark z_posterior = f_posterior.(μs', λs)
BenchmarkTools.Trial: memory estimate: 1.01 MiB allocs estimate: 51017 -------------- minimum time: 2.796 ms (0.00% GC) median time: 2.912 ms (0.00% GC) mean time: 3.189 ms (4.28% GC) maximum time: 64.073 ms (92.81% GC) -------------- samples: 1564 evals/sample: 1
# Julia eval parse 版
eval(parse("g_posterior(μ,λ) = $posterior"))
@benchmark z_posterior = g_posterior.(μs', λs)
BenchmarkTools.Trial: memory estimate: 79.89 KiB allocs estimate: 5 -------------- minimum time: 1.162 ms (0.00% GC) median time: 1.183 ms (0.00% GC) mean time: 1.241 ms (2.11% GC) maximum time: 62.258 ms (98.05% GC) -------------- samples: 4028 evals/sample: 1
N(posterior)
# Julia eval parse N 版
eval(parse("h_posterior(μ,λ) = $(N(posterior))"))
@benchmark z_posterior = h_posterior.(μs', λs)
BenchmarkTools.Trial: memory estimate: 79.89 KiB allocs estimate: 5 -------------- minimum time: 1.148 ms (0.00% GC) median time: 1.220 ms (0.00% GC) mean time: 1.377 ms (2.50% GC) maximum time: 82.969 ms (98.31% GC) -------------- samples: 3621 evals/sample: 1
ちなみにJuliaにおいて meshgrid は次のセルのようにして実現できる.
しかし, meshgridは, 同一のデータを大量に複製することによって, メモリを無駄に使用する. 可能ならば避けた方がコンピューターの効率的な利用という観点からは合理的である.
# SymPy lambdify 版 meshgrid ケース
μ_grid = repeat(μs, 1, length(λs))'
λ_grid = repeat(λs, 1, length(μs))
f_posterior = lambdify(posterior, [μ, λ])
@benchmark z_posterior = f_posterior.(μ_grid, λ_grid)
BenchmarkTools.Trial: memory estimate: 1.01 MiB allocs estimate: 51016 -------------- minimum time: 2.809 ms (0.00% GC) median time: 2.862 ms (0.00% GC) mean time: 3.066 ms (4.38% GC) maximum time: 65.052 ms (95.40% GC) -------------- samples: 1630 evals/sample: 1
figure(figsize=(5,4))
pcolormesh(μs, λs, z_posterior, cmap="CMRmap")
colorbar()
grid(ls=":")
xlabel("μ")
ylabel("λ")
yscale("log")
title("posterior")
PyObject Text(0.5,1,'posterior')
かもしれない。