# 事後分布の例¶

2017年7月11日

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