元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

Julia tips #5: MCMC法(Metropolis-Haisting法)を実装してみた

JuliaにてMCMC法(Metropolis-Haistings法)を実装してみました。

  • コードの利用・改変・再配布は自由にして構いませんが、何か起きても責任は取りません。
  • 一部いらない変数や無駄実装が混じってます。動作に問題はありません。バイオ実験系なので許して...。
  • 一番下にコードを全部まとめたものを書いておきました。 パパっと試してみたい方はコピペしてどうぞ。

環境

実装

使うパッケージのインポート

using Distributions
using Gadfly
using StatsBase

目標分布を定義する

MCMCでよく出てくる有名なあれです。

μ = [-5.0, -5, 5.0]
σ = [1.0, 0.1, 1.0] 

function target(x::Array{Float64, 1})
    arr = zeros(length(x))
    Normal_val_1(x) = 0.45 * 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 0.45 * 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    Normal_val_3(x) = 0.1 * 1/sqrt(2π*σ[3]^2)*exp(-(x-μ[3])^2/(2σ[3]^2))
    for i in 1:length(x)
        arr[i] = Normal_val_1(x[i]) + Normal_val_2(x[i]) + Normal_val_3(x[i])
    end
    
    return arr
end

function target(x::Float64)
    arr = zeros(length(x))
    Normal_val_1(x) = 0.45 * 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 0.45 * 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    Normal_val_3(x) = 0.1 * 1/sqrt(2π*σ[3]^2)*exp(-(x-μ[3])^2/(2σ[3]^2))
    val = Normal_val_1(x) + Normal_val_2(x) + Normal_val_3(x)
    return val
end

多重ディスパッチを利用して、スカラーが来たら値を一つだけ返し、ベクトルが来たらベクトルを返すようになっています。
変な分岐処理を書かなくていいので楽ですね。

上記の関数を定義後に、このコードで分布を見ることができます。

x = Array{Float64}(linspace(-10, 10, 1000))
y = target(x)
p=plot(x=x, y=y, Geom.line)

こんな感じ。 f:id:hotoke-X:20180121050147p:plain

Sampler型を定義

Pythonのようにクラスを作って、メンバ変数としてMCMC法の状態を保持したいなーと思ったのですが、JuliaにはPythonのようなクラスがありません。*
無いならそれっぽい使い方ができればいいじゃん!
という発想です。

ここで定義した型をMCMCモジュールへ投げて、変数の保持に使います。

mutable struct Sampler
    x0::Float64
    x_present::Float64
    target::Any
    propose_next::Any
    burn_in::Int
    count::Int
    acceptance::Float64
    seed::Int
    random_state::MersenneTwister
end

各変数の説明
x0: 初期値
x_present: 現在の値を保持
target: 目標分布の密度関数
propose: 次点を生成する関数 burn_in: 焼きなまし期間
count: 現在のステップ数
acceptance: 採択回数(本当は採択率を入れたいが考え中。今回は不要。)
seed: 乱数の初期パラメータ
random_state: 乱数ジェネレータ(これを保持することで、途中からサンプリングを再開できる)

回してみる

seed = 0
rng = srand(seed)
gaussian = Normal(0,2)
uniform = Uniform(-10, 10)
inits = rand(uniform)
target = target
generator(x) = x + 0.5*rand(gaussian)
burn_in = 10000
num_sample = 1000000
cnt = 0
acceptance = 0

sampler = Sampler(inits, inits, target, generator, burn_in, cnt, acceptance, seed, rng)
burn_in = MCMC.burn_in(sampler)
samples = MCMC.run(sampler, num_sample)

結果のプロット

set_default_plot_size(9inch, 12inch/golden)

# プロット範囲
_start = 990000
_last = 1000000

p1 = plot(y=autocor(samples), Geom.line, Guide.ylabel("Autocorrelation\ncoefficient"),
    Guide.xlabel("Lag"), Guide.title("Autocorrelation"))
p2 = plot(x=_start:_last, y=samples[_start:_last], Geom.line,
    Guide.xlabel("Steps"), Guide.ylabel("Value"), Guide.title("Trajectory"))
p3 = plot(layer(x=x, y=y, Geom.line, Theme(default_color="magenta")),
     layer(x=samples, Geom.histogram(density=true)),
     Guide.ylabel("Frequency"), Guide.xlabel("Samples"),
     Guide.title("Posterior distribution"),
     Guide.manual_color_key("Color", ["Target", "Posterir"], ["magenta", "cyan"]))

p = title(vstack(hstack(p1, p2), p3), "Metropolis-Hastings (σ=2)",  Compose.fontsize(14pt))

f:id:hotoke-X:20180121050013p:plain

できました。

全コードまとめ

using Distributions
using Gadfly
using StatsBase
using Cairo

# 目標分布の定義--------------------------
μ = [-5.0, -5, 5.0]
σ = [1.0, 0.1, 1.0] 

function target(x::Array{Float64, 1})
    arr = zeros(length(x))
    Normal_val_1(x) = 0.45 * 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 0.45 * 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    Normal_val_3(x) = 0.1 * 1/sqrt(2π*σ[3]^2)*exp(-(x-μ[3])^2/(2σ[3]^2))
    for i in 1:length(x)
        arr[i] = Normal_val_1(x[i]) + Normal_val_2(x[i]) + Normal_val_3(x[i])
    end
    
    return arr
end

function target(x::Float64)
    arr = zeros(length(x))
    Normal_val_1(x) = 0.45 * 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2))
    Normal_val_2(x) = 0.45 * 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2))
    Normal_val_3(x) = 0.1 * 1/sqrt(2π*σ[3]^2)*exp(-(x-μ[3])^2/(2σ[3]^2))
    val = Normal_val_1(x) + Normal_val_2(x) + Normal_val_3(x)
    return val
end

# Sampler型を定義(MCMCモジュール内で状態を保持する)
mutable struct Sampler
    x0::Float64
    x_present::Float64
    target::Any
    propose_next::Any
    burn_in::Int
    count::Int
    acceptance::Float64
    seed::Int
    random_state::MersenneTwister
end

# MCMC法の本体
module MCMC
    function present(sampler)
        return sampler.x_present
    end

    function next(sampler)
        proposed = sampler.propose_next(sampler.x_present)
        odds = sampler.target(proposed)/sampler.target(sampler.x_present)
    
        if rand(sampler.random_state) < odds
            sampler.x_present = proposed
            sampler.acceptance += 1
            return proposed
        else
            return sampler.x_present
        end
    end

    function burn_in(sampler)
        arr = zeros(sampler.burn_in)
        for i in 1:sampler.burn_in
            arr[i] = next(sampler)
            sampler.count += 1
        end
        sampler.acceptance = 0
        return arr
    end

    function run(sampler, n_steps)
        arr = zeros(n_steps)
        for i in 1:n_steps
            arr[i] = next(sampler)
            sampler.count += 1
        end
        return arr
    end

end

# 実行コード
seed = 0
rng = srand(seed)
gaussian = Normal(0,2)
uniform = Uniform(-10, 10)
inits = rand(uniform)
target = target
generator(x) = x + 0.5*rand(gaussian)
burn_in = 10000
num_sample = 1000000
cnt = 0
acceptance = 0

sampler = Sampler(inits, inits, target, generator, burn_in, cnt, acceptance, seed, rng)
burn_in = MCMC.burn_in(sampler)
samples = MCMC.run(sampler, num_sample)


# 結果のプロット
set_default_plot_size(9inch, 12inch/golden)

_start = 990000
_last = 1000000

p1 = plot(y=autocor(samples), Geom.line, Guide.ylabel("Autocorrelation\ncoefficient"),
    Guide.xlabel("Lag"), Guide.title("Autocorrelation"))
p2 = plot(x=_start:_last, y=samples[_start:_last], Geom.line,
    Guide.xlabel("Steps"), Guide.ylabel("Value"), Guide.title("Trajectory"))
p3 = plot(layer(x=x, y=y, Geom.line, Theme(default_color="magenta")),
     layer(x=samples, Geom.histogram(density=true)),
     Guide.ylabel("Frequency"), Guide.xlabel("Samples"),
     Guide.title("Posterior distribution"),
     Guide.manual_color_key("Color", ["Target", "Posterir"], ["magenta", "cyan"]))

p = title(vstack(hstack(p1, p2), p3), "Metropolis-Hastings (σ=2)",  Compose.fontsize(14pt))

*burn in 期間が短いときに途中から再計算させたいので。