事後分布の例

黒木玄

2017年7月11日

渡辺澄夫著『ベイズ統計の理論と方法』の1.4節の内容を再現したい。

In [1]:
# 乱数のseed値を指定

srand(9869)
Out[1]:
MersenneTwister(UInt32[0x0000268d], Base.dSFMT.DSFMT_state(Int32[-1985563288, 1073428835, -1636231178, 1073332724, -1816778846, 1073204032, 158087370, 1073236640, 762020198, 1072698358  …  -426643226, 1073209542, -1875029544, 1072985977, 1548904078, 104009346, -1256843749, -457095338, 382, 0]), [1.65121, 1.4691, 1.49154, 1.0495, 1.66047, 1.99551, 1.64181, 1.70037, 1.03812, 1.47842  …  1.8145, 1.55174, 1.25966, 1.98041, 1.07741, 1.38703, 1.90453, 1.72587, 1.2752, 1.19044], 382)
In [2]:
using Distributions
import PyPlot
plt = PyPlot
hist = plt.plt[:hist]
hist2D = plt.plt[:hist2d]

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end
Out[2]:
meshgrid (generic function with 1 method)

混合正規分布モデル

$$ Y \sim (1-a)\,\mathrm{Normal}(0,1) + a\,\mathrm{Normal}(b,1) $$
In [3]:
mixmodel(a,b) = MixtureModel(Normal, [(0.0,1.0), (b,1.0)], [1-a, a])
Out[3]:
mixmodel (generic function with 1 method)
In [4]:
log_likelihood(a,b,sample) = loglikelihood(mixmodel(a,b),sample)

# Sample is generated by the model for the parameters a0 and b0
#
function plot_sample_and_posterior(a0,b0,N)
    println("a0 = $a0, b0 = $b0, N = $N")
    
    sample = rand(mixmodel(a0,b0),N)
    
    plt.figure(figsize=(5,3))
    hist(sample, bins=Int(floor(sqrt(N))))
    plt.title("Sample histogram")
    
    a = [linspace(0,1,100);]
    b = [linspace(-1,2*b0+1,100);]
    aa, bb = meshgrid(a,b)
    c = log_likelihood(a0,b0,sample)
    z = ((a,b) -> exp(log_likelihood(a,b,sample) - c)).(aa, bb)

    plt.figure()
    plt.surf(a,b,z, cmap=plt.ColorMap("rainbow"))
    plt.title("Posterior distribution")
    plt.xlabel("a")
    plt.ylabel("b")
    
    plt.figure()
    plt.axes(facecolor="black")
    plt.pcolormesh(a,b,z, cmap="rainbow")
    plt.title("Posterior distribution")
    plt.xlabel("a")
    plt.ylabel("b")
    plt.colorbar()
end
Out[4]:
plot_sample_and_posterior (generic function with 1 method)
In [5]:
@time plot_sample_and_posterior(0.5, 3.0, 100);
a0 = 0.5, b0 = 3.0, N = 100
  6.553598 seconds (39.17 M allocations: 727.143 MiB, 0.99% gc time)
In [6]:
@time plot_sample_and_posterior(0.5, 1.0, 100);
a0 = 0.5, b0 = 1.0, N = 100
  1.905923 seconds (37.90 M allocations: 658.169 MiB, 2.88% gc time)
In [7]:
@time plot_sample_and_posterior(0.5, 0.5, 100);
a0 = 0.5, b0 = 0.5, N = 100
  1.857465 seconds (37.90 M allocations: 658.169 MiB, 2.25% gc time)
In [8]:
@time plot_sample_and_posterior(0.5, 0.2, 100);
a0 = 0.5, b0 = 0.2, N = 100
  1.843625 seconds (37.90 M allocations: 658.169 MiB, 2.26% gc time)
In [9]:
@time plot_sample_and_posterior(0.5, 0.1, 100);
a0 = 0.5, b0 = 0.1, N = 100
  1.912925 seconds (37.90 M allocations: 658.169 MiB, 2.30% gc time)
In [10]:
@time plot_sample_and_posterior(0.5, 0.1, 1000);
a0 = 0.5, b0 = 0.1, N = 1000