⚡ Pluto.jl ⚡

Asymmetric MCMC

The Asymmetirc Metropolis-Hastings Sampler

6.8 μs

xxxxxxxxxx
 
using Plots, Distributions, LinearAlgebra,Random,LaTeXStrings,StatsBase,DataFrames

10.2 ms

xxxxxxxxxx
 
begin
    using REPL
    REPL.REPLCompletions.latex_symbols["\\_y"] = "ᵧ"
    nothing
end

2.5 ms

Model parameter : θ = [α,β]
Data : y
Prior : p(θ)
Likelihood : p(y|θ) = Gamma(y;α,β)
Posterior : p(θ|y) = p(y|θ)p(θ) / p(y)

Gamma(y;α,β) = βᵅ yᵅ⁻¹e⁻ᵝʸ / Γ(α)
Gamma(y;α,β) = pdf(Gamma(α,β),y) : Julia

5.6 μs

Likelihood : p(y|α,β=1) = p_y_α(y,α,β)

4.1 μs

p_y_α (generic function with 1 method)

xxxxxxxxxx
 
p_y_α(y,α,β) = pdf(Gamma(α,β),y)

20.3 μs

초기값 및 상수설정

3.7 μs

xxxxxxxxxx
 
begin
    Random.seed!(0)
    n=50
    α = LinRange(0.1,5,n)
    Δα = (α[end] - α[1])/(n-1)
    #주어진 값에 가장 가까운 α 의 index구하기
    Iₐ(x) = Int(ceil(x/Δα))
    αₖ = 2
    αᵢ = repeat([ α[Iₐ(αₖ)] ],n)
    β = 1
    y = LinRange(0,10,n)
    Δy = (y[end] - y[1])/(n-1)
    #주어진 값에 가장 가까운 y 의 index구하기
    Iᵧ(x) = Int(ceil(x/Δy))
    yₖ = 1.5
    yᵢ = repeat([ y[Iᵧ(yₖ)] ],n)    
    μ = 5
end

110 ms

Likelihood 그래프 (β=1)

5.0 μs

	α	y	p
	Float64	Float64	Float64
1	0.1	0.0	Inf
2	0.1	0.204082	0.358265
3	0.1	0.408163	0.156547
4	0.1	0.612245	0.0886201
5	0.1	0.816327	0.0557771

xxxxxxxxxx
 
begin
    Dₚ = DataFrame(α=repeat(α, inner=length(α)), y=repeat(y,outer=length(y)))
    Dₚ.p = p_y_α.(Dₚ.y,Dₚ.α,β)
    first(Dₚ,5)
end

423 ms

xxxxxxxxxx
 
begin   
    plot(Dₚ.α,Dₚ.y,Dₚ.p,st=:surface,fc=:heat,
        xlabel="α",ylabel="y",zlabel="p(y|α,β=1)",camera=(40,70),
        title="Likelihood Gamma distribution, β=1 surface") 
    plot!(αᵢ,y,[p_y_α(yᵢ,αₖ,β) for yᵢ in y  ],linewidth=3,linecolor=:red,
        label="p(y|α=2,β=1)")   
end

1.8 s

100 ns

모델 파라미터 θ = [α,β]에 대한 사전확률분포(prior)를 다음과 같이 가정한다.
P(β=1) = 1
p(α) = sin(π×α)²
두 모델 파라미터의 사전확률분포는 적분값이 1이 아니기 때문에 부적절한 사전확률분포(improper priors)라고 한다.
p(α|y) = p(y|α)p(α)/p(y) 식에서 normalization constant p(y)는 계산이 어려울 수 있어나 p(y)에 대해 알지 못해도 Metropolis-Hastings 알고리즘을 사용하여 사후확률분포 p(α|y)에서 α를 샘플링 할 수 있다.
특히 정규화 상수 p(y)를 무시하고 비정규화된 사후확률분포(posterior)에서 α를 샘플링 할 수 있다.
p(α|y) ∝ p(y|α)p(α) = Gamma(y;α,β)×sin(π×α)² = pdf(Gamma(α,β),x)×sin(π×α)²

3.7 μs

Posterior : p(α,β=1|y) = p_α_y(y,α,β) = p_y_α(y,α,β) * p_α(α)

3.4 μs

p_α_y (generic function with 1 method)

xxxxxxxxxx
 
begin
    p_α(α) = sin(π*α)^2
    p_α_y(y,α,β) = p_y_α(y,α,β)*p_α(α)
end

82.1 μs

Posterior Surface

2.9 μs

	α	y	p
	Float64	Float64	Float64
1	0.1	0.0	Inf
2	0.1	0.204082	0.0342113
3	0.1	0.408163	0.0149489
4	0.1	0.612245	0.00846246
5	0.1	0.816327	0.00532624

xxxxxxxxxx
 
begin
    D = DataFrame(α=repeat(α, inner=length(α)), y=repeat(y,outer=length(y)))
    D.p = p_α_y.(D.y,D.α,β)
    first(D,5)
end

58.7 ms

xxxxxxxxxx
 
begin
    plot(D.α,D.y,D.p,st=:surface,fc=:heat,xlabel="α",ylabel="y",zlabel="p(α,β=$β|y)",
        camera=(30,70),
        title="Posterior distribution") 
    # reshape α,y 순서 중요함,row,column개념으로 p 데이터가 생성 되었음
    # Iᵧ(yₖ) : yₖ에 가장 가까운 y의 index를 가져옴
    iᵧ = Iᵧ(yₖ)
    # the column-major iteration order 
    # julia는 reshape시에 column을 우선시 해서 row를 채우고 다음 column으로 넘어간다.
    # 예) p₁₁,p₁₂,p₁₃,p₂₁,p₂₂
    pᵧₖ = reshape(D.p,(length(y),length(α)))[iᵧ,:]
    plot!(α,yᵢ,pᵧₖ,linewidth=3,linecolor=:red,
            label="p(α,β=$β|y=$yₖ)")    
    plot!(α,repeat([y[end]],n),[p_α(αᵢ) for αᵢ in α  ],linewidth=3,linecolor=:blue,
            label="p(α,β=$β)")  
end

212 ms

0.0 ns

Asymmetric Proposal Distribution

$\begin{aligned} f (y; μ) = μ e^{- μ y} = \frac{1}{θ} e^{- \frac{y}{θ}}, x > 0, μ = \frac{1}{θ} \end{aligned}$

f(y;μ) = q(y,μ) = pdf(Exponential(μ),y)

3.6 μs

q (generic function with 1 method)

xxxxxxxxxx
 
q(y,μ) = pdf(Exponential(μ),y)

26.4 μs

xxxxxxxxxx
 
begin
    plot(α,α->p_α_y(yₖ,α,β),linecolor=:blue,linewidth=3,
        label="Target, p(α,β=$β|y=$yₖ)")
    plot!(α,y->q(y,μ),linewidth=3,linecolor=:red,
        label="Proposal, q(α,$μ)")
end

361 ms

Run Metropolis-Hastings Sampler

x = α 를 샘플링 한다.
β = 1
yₖ = 1.5
μ = 5

15.7 μs

xxxxxxxxxx
 
begin
    # Samples
    Nₛ = 5000
    n_burn = 500
    min_α = α[1]
    max_α = α[end]
    #Initialize sampler
    αₛ = zeros(Nₛ)
    αₛ[1] = μ
    t = 1
end

9.7 ms

xxxxxxxxxx
 
while t < Nₛ
    t = t+1
    #sample from proposal
    #제안 분포로 부터 α를 뽑는다
    # α\asteraccent => α⃰
    α⃰=rand(Exponential(μ))
    
    # Correction factor α
    c = q(αₛ[t-1],μ)/q(α⃰,μ)
    
    # Calculate the (correcte) Acceptance Ratio
    #α → α⃰ 로 전이하는 것을 받아 들이는 비율
    A = min(1,p_α_y(yₖ,α⃰,β)/p_α_y(yₖ,αₛ[t-1],β) * c)
    
    #Accept or reject?
    u = rand(Uniform(0,1))
    if u < A
        αₛ[t] = α⃰
    else
        αₛ[t] = αₛ[t-1]
    end
end

74.2 ms

xxxxxxxxxx
 
begin
    p1=plot(αₛ,1:Nₛ,linetype=:step,yflip=true,xlim=(min_α,max_α),
        label="α samples",xlabel="Samples, α",ylabel="t",
        title="Markov Chain Path",titlefontsize=10)
    plot!(α,α->n_burn,linewidth=2,linestyle=:dash,linecolor=:red,
            label="Burnin")
    
    h = fit(Histogram, αₛ[n_burn:end],α)
    p2=plot(h.edges,h.weights./sum(h.weights),st=:bar,linewidth=0.1,
        xlim=(min_α,max_α),
        title="Samples",titlefontsize=10,xlabel="Samples, α",
        label="Sampled Distribution",legendfontsize=6,
        ylabel="p(α,β=$β|y=$yₖ)")
    plot!(α,pᵧₖ./sum(pᵧₖ),linewidth=2,
        label="Target Posterior")   
​
    plot(p1,p2,layout=(2,1))
end

1.1 s