『モンテカルロ統計計算』章末問題を Julia で解く:第3章「積分法」
『モンテカルロ統計計算』(講談社、鎌谷研吾 著、駒木文保 編)の章末問題を Julia で解くシリーズ第2回で、今回は第3章「積分法」です。解析計算するだけの問題は割愛。
3.4 複合台数公式
using Plots function composite_trapezoidal_rule(N) func(x) = (cos(50x) + sin(20x)) ^ 2 f = func.(collect(range(0, 1, length=N))) (2 * sum(f) - f[1] - f[end]) / 2(N - 1) end N = 100 I = composite_trapezoidal_integral.(2:N) plot(I) hline!([0.9605764394771032]) # analytical solution savefig("3.4.png")
3.5 基本的モンテカルロ積分法でのコーシー分布の裾確率計算
using Distributions, Plots function approx_ci(dist, c, M) f_x = rand(dist, M) .>= 2.0 I = mean(f_x) σ = std(f_x) l = I - σ * c / sqrt(M) u = I + σ * c / sqrt(M) return I, l, u end function run(M) α = 0.05 c = quantile(Normal(), 1 - α / 2) dist = Cauchy() I = zeros(length(M)) l = zeros(length(M)) u = zeros(length(M)) for (index, m) in enumerate(M) I_m, l_m, u_m = approx_ci(dist, c, m) I[index] = I_m l[index] = l_m u[index] = u_m end return I, l, u end M = map(x -> ^(10, x), [1, 2, 3, 4, 5, 6]) .|> Int I, l, u = run(M) println("Estimated: ", I[end]) println("lower bound: ", l[end]) println("upper bound: ", u[end]) # Estimated: 0.147537 # lower bound: 0.14684191718484835 # upper bound: 0.14823208281515166 plot(M, I, xscale=:log10, frange=l, falpha=0.2, fcolor=:black, lab=nothing) plot!(M, I, xscale=:log10, lalpha=0, frange=u, falpha=0.2, fcolor=:black, lab=nothing) savefig("3.5.png")
信頼係数95%で小数点第2位まで求めるにはサンプル数を 1,000,000 にする必要があった。一桁少ないと精度保証できない。
3.6 コーシーの裾確率計算を利用した逆正接関数の近似
using Distributions, Plots function gen_mc(M, dist) function mc(x) I = mean(rand(dist, M) .>= x) return π * (1 / 2 - I) end mc end x = range(-10, 10, length=101) dist = Cauchy() M = map(x -> ^(10, x), [1, 2, 3, 4]) .|> Int data = zeros(length(x), length(M)) for (index, m) in enumerate(M) approx_atan = gen_mc(m, dist) data[:, index] = approx_atan.(x) end for (index, m) in enumerate(M) plot!(x, data[:, index], lab = "$m", legend = :bottomright) end plot!(x, atan.(x), lab = "truth", legend = :bottomright) savefig("3.6.png")
各点でのサンプル数が 10,000 くらいになると、真値とほとんど変わらず推論できているように見える。
3.8 コーシー分布の事後平均の自己正規化モンテカルロによる計算
using Distributions, Plots function likelihood(x::Array{Float64, 1}, θ::Float64) p(y) = pdf(Cauchy(θ), y) prod(p.(x)) end "Self-nomalized Monte Carlo" function snmc(x::Array{Float64, 1}, params::Array{Float64, 1}) l = likelihood.([x], params) numerator = sum(l .* params) denominator = sum(l) numerator / denominator end function demo() θ = 1.0 d_obs = Cauchy(θ, 1) d_prior = Cauchy() N = 100 M = map(x -> ^(10, x), [1, 2, 3, 4, 5]) .|> Int I = zeros(length(M)) for (index, m) in enumerate(M) x = rand(d_obs, N) params = rand(d_prior, m) I[index] = snmc(x, params) end println(I) plot(M, I, xscale = :log10) savefig("3.8.png") end demo()
データ点 の数(ここでは 100 で固定)を増やすと自己正規化する際の分子・分母がともに非常に小さくなり、割り算の結果が NaN となるのを確認した(各データに対する尤度が 1 未満のため、それらの積はデータ数が増えるごとに小さくなる)。数値的に不安定な方法?
3.9 重点サンプリング法
using Distributions, Plots function approx_cauchy(M::Int) x = rand(Cauchy(), M) f(y) = exp(-√abs(y)) * (sin(y)) ^ 2 * π * (1 + y ^ 2) mean(f.(x)) end function approx_normal(M::Int) x = rand(Normal(), M) f(y) = sqrt(2π) * exp(y ^ 2 / 2 -√abs(y)) * (sin(y)) ^ 2 mean(f.(x)) end function demo() M = map(x -> ^(10, x), [2, 3, 4, 5, 6]) .|> Int values_cauchy = approx_cauchy.(M) values_normal = approx_normal.(M) plot(M, values_cauchy, xscale = :log10, lab = "Cauchy") plot!(M, values_normal, xscale = :log10, lab = "Normal") savefig("3.9.png") end demo()
結果を見るとコーシー分布を使って重点サンプリングした方が安定している(正規分布の方は暴れすぎ)。いちおう被積分関数を見てみると、
となっている。