Julia tips #5: MCMC法(Metropolis-Haisting法)を実装してみた
JuliaにてMCMC法(Metropolis-Haistings法)を実装してみました。
- コードの利用・改変・再配布は自由にして構いませんが、何か起きても責任は取りません。
- 一部いらない変数や無駄実装が混じってます。動作に問題はありません。バイオ実験系なので許して...。
- 一番下にコードを全部まとめたものを書いておきました。 パパっと試してみたい方はコピペしてどうぞ。
環境
- Windows 10 HOME
- Julia 0.6.2
実装
使うパッケージのインポート
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)
こんな感じ。
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))
できました。
全コードまとめ
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 期間が短いときに途中から再計算させたいので。