MCMCへの道（２）

 MCMCへの道（１）からの続き

さて，MCMC法#2棄却サンプリングである。手順は次の通り。

(1) 目標分布 p(x) を事後分布ともいう。これは既知の関数。

(2) 乱数発生計算が容易な分布 q(x) で， k q(x) ≧ p(x) を満たすものを選ぶ。 kは正定数。例えば，一様分布とか正規分布。

(3) q(x) にしたがう乱数 x' を発生する。

(4) [0, k q(x)] の一様分布にしたがう乱数 y' を発生する。

(5) y' ≦ p(x') ならx' を採用する。そうでなければ棄却する。

この(3) から(5) を必要なだけ繰り返す。

using Plots

function beta(x,a,b)
  return x^(a-1)*(1-x)^(b-1)
end

x= [k for k in 0.0:0.001:1.0]
y= beta.(x,10.2,5.8)
z= [0.00013 for k in 0.0:0.001:1.0]
scatter(x,y,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :blue,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

g1=scatter!(x,z,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :orange,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

n=3000
p=rand(n)
r=0.00013*rand(n)
q=ones(n)
for i in 1:n
  f=beta(p[i],10.2,5.8)
  q[i]=ifelse(f <= r[i],0,r[i])
end

g2=scatter(p,q,xlim=(0,1.0),ylim=(0,0.00015),
  markershape = :circle,
  markersize = 1,
  markeralpha = 0.75,
  markercolor = :green,
  markerstrokewidth = 0.1,
  markerstrokealpha = 0.1,
  markerstrokecolor = :white)

plot(g1,g2,legend=:left,aspect_ratio=6000)

図：ベータ分布 p(x) と一様分布 q(x) と定数 k=0.00013 による棄却法の例