『モンテカルロ統計計算』章末問題を 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()

結果を見るとコーシー分布を使って重点サンプリングした方が安定している（正規分布の方は暴れすぎ）。いちおう被積分関数を見てみると、

となっている。