Julia tips #9: CVODE (Sundials.jl) の数値積分が終わらない時の回避方法
追記
どうやら、せっかく設定したhminは内部でフラグチェックが行われないようです。
solve関数のmaxitersオプションで対応するしかないかもしれません。 また後程検証してみます。
以下は検証が十分ではありません。
ソルバ自体の動作に支障はありませんが、期待通りに動いているかは検証する必要があります。
そのことをご理解いただいたうえでご覧ください。
2018/03/08当時の環境
- Julia 0.6.2
- Sundials 1.2.0
事の発端
ODEモデルのパラメータ推定をしていたら、非線形なモデルだとどうしても途中でソルバがフリーズしてしまう...。
原因はソルバが無限に数値積分のステップを続けてしまうことでした。
数値計算に疎いのでどうしてそんなことになってしまうのかはちゃんとわかっていないのですが、stiffなODEを解こうとすると、数値積分時に細かいステップを刻まないといけなくなるのが原因です。
「DifferentialEquationsから簡単にSundialsが使える！」と喜んだのもつかの間。
DifferentialEquations.jlでもこの件について触れられていました。
Frequently Asked Questions · DifferentialEquations.jl
公式によると、現状アルゴリズムを変えたり、ソルバを変えたりするしかないようです。
しかし、stiffなODEが解けないような状況においてCVODE以外ではお話にならないのです（経験上）。
そもそもODEが解けないようなパラメータは捨ててしまえばよいので、途中で計算を止めてなんとか捨てたいわけです。
Sundials? CVODE?
Sundialsは常微分方程式のソルバ、CVODEを作っているところです（他に幾つも作ってる）。多用しないのでわかりませんが、良く名前を聞くのでシェアは大きいと思います、MATLABにもあります。
stiffな問題をCVODEに解かせたら他のソルバより圧倒的に高速で数値積分してくれます。
stiffでなくてもめっちゃ速いです。過去記事中のリンク先でJuliaのソルバで速度比較をした方がいます。
ちなみにCVODEとよく比較されるソルバにLSODAっていうのがあります。自分の場合はCVODEの方が5倍くらい速かったです。問題にもよると思うので、CVODEがダメだったら試す価値はあるかもしれません。JuliaならPkg.add("LSODA")ってやればDiffEq.jlのインターフェースから使えるようになります。 以下のブログ記事ではLSODAも良いぞ！って言ってますね。
CVODEが途中で止まる
最強のCVODEといえど限界はあります。積分のステップが小さすぎると計算が終わらなくなります。途中でエラー吐けよって思うんですが、デフォルト設定だと無限にステップを刻めることになっています。こいつをちょちょいといじれば良いのです。
以下のhminというオプションをいじる（CVODE公式ドキュメント, p37より）↓
The default value is 0.0となっていますね。
これがいかんわけです。
Sundials.jlを書き換える（自己責任で）
DiffEq.jlのインターフェースには残念なことに、hminを設定するオプションがありません。しかし、Sundials.jl自体ではちゃんとラップしてくれているのです。じゃあ直接Sundials.jlを使えばいいんですが、exampleを見てるとわかるんですがjuliaの可読性が損なわれます。何より面倒です。
なので、DiffEq.jlのインターフェースから使えるようにSundials.jlを少しだけいじりました。
自分でやる場合は自己責任でお願いします。 ...といっても、ミスっても書き換えたところを戻せば元通りです（2敗）
DiffEq.jlのインターフェースとなるsolve.jlをSundialsのパッケージディレクトリから探します。
hminを引数としてsolveに追加（後ろにkwargs...があるので多分これやらなくても動くけど）。
以下のようにhminをCVODEの設定として反映させる（ハイライト部分）。
これで完了。
SundialsはCのライブラリなんですが、Sundilas.jlの作者さんがきれいにラップしてくれているおかげでインターフェース部分を少し書き換えるだけで良いです。作者さんありがとうございます。
後は動作確認したらよいです。
パッケージを書き換えたので、初めにプリコンパイルが走ります。
（以下ではCVODEが固まるほどのstiffなODEがどんなものかわからないので、ソルバが動くかどうかだけ確かめています。）
using DifferentialEquations using Sundials using Plots function parameterized_lorenz(du,u,p,t) du[1] = p[1]*(u[2]-u[1]) du[2] = u[1]*(p[2]-u[3]) - u[2] du[3] = u[1]*u[2] - p[3]*u[3] end u0 = [1.0,0.0,0.0] tspan = (0.0,100.0) p = [10.0,28.0,8/3] alg = CVODE_BDF() prob = ODEProblem(parameterized_lorenz, u0, tspan, p) sol = solve(prob, hmin=1.0e-8) plot(sol,vars=(1,2,3))
以下のようにグラフがプロットされればとりあえず正常に動作しています。
しかし、本当に設定が反映されたのだろうか...
（どなたか、CVODEをフリーズさせる程のめちゃめちゃ硬いODEを教えてください）
Julia tips #8: パッケージの追加に失敗したとき
自宅PCにDifferentialEquations.jlを追加しようとしたら以下のようなエラーが...
julia> Pkg.add("DifferentialEquations") ERROR: resolve is unable to satisfy package requirements. The problem was detected when trying to find a feasible version for package DiffEqPhysics. However, this only means that package DiffEqPhysics is involved in an unsatisfiable or difficult dependency relation, and the rootof the problem may be elsewhere. (you may try increasing the value of the JULIA_PKGRESOLVE_ACCURACY environment variable)
パッケージのアップデートしたら解決するだろうと
Pkg.update()
もしてみたが解説せず。
同じことで困っている人を探したら以下のような投稿を発見。
いっぱいやり取りしてくれてますが、要は
julia>ENV["JULIA_PKGRESOLVE_ACCURACY"] = 10 julia>Pkg.resolve()
とするだけ。
何が起きてるのかはよくわかりませんが、バージョンの依存関係をいい感じに直してくれてるみたい？
後はPkg.addで好きなパッケージを追加していきましょう。
Julia tips #7: Juliaでプロットのlegendを手動で変更したい
Juliaのプロットにて、legendの変更で躓いたので色んなバックエンドで簡単に試してみた。
Plotlyでのプロットが消えてしまいましたが、まぁいいでしょう。
legendをプロットの外へ移動したい場合はどうすればいいんだろうか...
解決しました。
これではだめですか？ https://t.co/KTMiog5fvY
— ktrst (@ktrst) 2018年2月12日
@ktrstさん、ありがとうございます。
Jupyter notebookのフォントを変える（Windows 10）
Juliaで下付き文字を使うときにあまりにも見にくかったので、フォントを変えて多少なり誤魔化したいというのが事の発端。
幾らでも情報は転がっていますが、ぱっと調べた感じだとWindowsでいじっているのを見かけなかったので、念のためのメモ書きです。
以下のように記述したファイルをcustom.cssという名前で.jupyter/custom以下に保存します。
.jupyterファイルは大抵の人はユーザフォルダ直下にあると思います。
.CodeMirror pre { font-family: Consolas; font-size: 12pt; }
WindowsなのでConsolasフォントを指定しています。
好きなフォントにしたらいいと思います。
また、フォントサイズを少し大きくして多少誤魔化しています。
Windows 10でJulia専用のPython仮想環境を作る
Juliaでは、Pythonを呼び出して使うことができます。しかし、何も考えずに使おうとするとMinicondaのPython 2.7がインストールされてしまいます（2018/02/01現在）。
そこで、Anacondaで作ったJulia専用の仮想環境を作って、その中でJulia用のPythonパッケージ管理をしようという話です。
以下のような人を対象としています。
Anacondaをダウンロードしてインストール
既にAnacondaをインストールしている人は読み飛ばしてください。
以下からWindows用のAnacondaインストーラをダウンロードしてインストールします。 https://www.anaconda.com/download/
インストールの仕方は公式のこちらで確認できます。 https://docs.anaconda.com/anaconda/install/windows
インストールの際、AnacondaのPythonをPathへ追加するかどうかのチェック画面が現れます。
WindowsにはデフォルトでPythonが入っていないので追加しちゃってよいでしょう。
両方にチェックを入れときましょう。
（嫌な人はAnacondaを起動するときは、スタートメニューからAnaconda プロンプトを起動してください。）
Julia専用のAnaconda仮想環境の構築
コマンドプロンプト（またはAnacondaプロンプト）にて、以下のように入力して仮想環境を作ります。
ついでにjupyterも一緒にインストールしておきましょう。
C:\Users\username>conda create -n env_name python=3 jupyter
env_nameに自分の好きな仮想環境名を入れてください。
次に、activateして仮想環境の中に入ります。
また、Juliaで使いそうなパッケージをインストールしておきます。
C:\Users\username>activate env_name C:\Users\username>conda install numpy pandas scipy matplotlib
念のため、Jupyter notebookのカーネルに作った仮想環境を追加しておく。
（多分やらなくても良い）
※必ず仮想環境に入った状態で行ってください。
C:\Users\username>ipython kernel install --user --name=env_name --display-name=env_name
仮想環境のPythonのパスを確認
C:\Users\username>where python C:\Users\username\Anaconda3\envs\env_name\python.exe C:\Users\username\Anaconda3\python.exe
仮想環境とAnacondaインストール時に入るPythonのパスが表示されます。
envsの下にあるpythonが仮想環境のpythonなのでパスをコピーしておきます。
仮想環境から抜けて、jupyterのパスを取得します。
C:\Users\username>deactivate C:\Users\username>where jupyter C:\Users\username\Anaconda3\Scripts\jupyter.exe
表示されたパスを環境変数JUPYTERへ追加します。
コントロールパネル＞システムとセキュリティ＞システム＞システムの詳細設定＞環境変数(N)
で環境変数の設定画面に行けます。
上段のユーザー環境変数に新規でJUPYTERという名前で環境変数を作成し、jupyterのパスを入力します。
Juliaが勝手に環境変数を読んでjupyterを見つけてくれます。
jupyter.exeまで入力しないと無視されるので注意
以上でPython側の準備は終わりです。
Juliaをインストール
公式からWindows用のJuliaをダウンロードしてインストール。
インストール場所はどこでも良いです。
https://julialang.org/downloads/
このままだとPCがjuliaを見つけられず
C:\Users\username>julia
と入力しても何も起きません。パスを教えてやる必要があります。
インストールしたフォルダ/binをユーザー環境変数Pathへ追加。
（julia.exeがあるのでわかりやすいと思います）
juliaのパッケージをインストールするフォルダの場所をデフォルトから変えたい場合は、環境変数JULIA_PKGDIRを作って好きなパスを入力して下さい。
追加の仕方は環境変数JUPYTERの設定の時と同じです。
次に、コマンドプロンプトを立ち上げなおしてjuliaを起動します。
まず初めにJuliaの環境変数を書き換えて、先ほど作った仮想環境のPythonの場所を教えておきます。
julia>ENV["PYTHON"] = "C:/Users/userneme/Anaconda3/envs/env_name/python.exe"
juliaの設定を初期化します。
PyPlotを入れるとPyCallが一緒についてきて、その際に仮想環境のpythonを見つけて使ってくれるようになります。
julia>Pkg.init() julia>Pkg.add("IJulia") julia>Pkg.add("PyPlot") julia>using PyCall
以上で仮想環境のpythonを使うことができました。
まとめ
- 環境が汚れなくてハッピー！
- Julia自体でも仮想環境を作れたら楽。
Julia tips #6: Multiple-Try Metropolisを実装してみた
Multiple-Try Metropolis (MTM)って？
メトロポリス・ヘイスティングスの多数決バージョンです。
今いる場所から次の候補を複数生成し、その値周辺が今いる場所より尤もらしければ候補の中の1つへ移動します。
細かい話を書く元気がないので、理論的な背景は気が向いたらまたまとめます。
というわけで実装して行きましょう。
MTMアルゴリズム
任意の提案関数を定義し
を満たすとする。
また、を
で対称な任意の関数であるとする。
マルコフ連鎖の現在の状態をとおくと、MTMアルゴリズムは次のようになる。
から
を独立に生成する
重み
に従って
から
をサンプリングして計算する。
から
をサンプリングし、
とする。
を以下の確率で受容する。
今回はベイズ計算統計学（統計解析スタンダード）p93 例3.7を実装してみた。提案関数と対称関数は以下のようになっている。
$$T(x, y) = x_t + \mathcal{N}\left(0, 5\sqrt{2}\right)$$ $$\lambda(\boldsymbol{y}, \boldsymbol{x})=\left(\frac{T(x, y)+T(y, x)}{2}\right)^{-1}$$
目的関数は、
$$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$$
である。
使うパッケージのインポート
using Distributions using StatsBase import Plots; plt = Plots; plt.pyplot() using LaTeXStrings # グラフ周りの設定 fntsm = plt.font("serif") fntlg = plt.font("serif", 20) plt.default(titlefont=fntlg, guidefont=fntsm, tickfont=fntsm, legendfont=fntsm) plt.default(size=(800, 600))
目標分布の定義
$$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$$
const μ = [-5.0, 5.0] const σ = [1.0, 1.0] function target(x::Array{Float64, 1}) arr = zeros(length(x)) Normal_val_1(x) = 1/sqrt(2π*σ[1]^2)*exp(-(x-μ[1])^2/(2σ[1]^2)) Normal_val_2(x) = 1/sqrt(2π*σ[2]^2)*exp(-(x-μ[2])^2/(2σ[2]^2)) for i in 1:length(x) arr[i] = 1/3 * Normal_val_1(x[i]) + 2/3 * Normal_val_2(x[i]) end return arr end function target(x::Float64) 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)) val = 1/3 * Normal_val_1(x) + 2/3 * Normal_val_2(x) return val end x = Array{Float64}(linspace(-10, 10, 1000)) y = target(x) l = plt.@layout [a b; c{0.2h}] eq = L"$\frac{1}{3} \mathcal{N}(-5,1) + \frac{2}{3} \mathcal{N}(5, 1)$" p = plt.plot(x, y, title="Target distribution", lab=eq)
実装
まず、MTMに必要な関数やら値を保持する型を作っとく。
mutable struct Sampler
x_present::Float64 # 現在の値
target::Function # 目標分布関数
propose_next::Function # 提案関数
q::Function # (x_t+1, x_t)の順で渡す # 推移核（上で書いたT(x, y)）
λ::Function # 対称関数
burn_in::Int # 焼きなましの回数
count::Int # 何回目の更新か保持
random_state::MersenneTwister # 乱数ジェネレータ
end
MTM本体
module MTM using StatsBase function present(sampler) # 現在の値を返す関数 return sampler.x_present end function next(sampler, n_candidate) # MTM w_x2y = zeros(n_candidate) w_y2x = zeros(n_candidate) proposed = zeros(n_candidate) proposed_r = zeros(n_candidate) accept = 0 # Σⱼ w(yⱼ,xₜ) for i in 1:n_candidate # q(yⱼ|xₜ)でyⱼをサンプリング proposed[i] = sampler.propose_next(sampler.x_present) # q(yⱼ|xₜ) と q(xₜ|yⱼ) を計算 q_x2y = sampler.q(proposed[i], sampler.x_present) q_y2x = sampler.q(sampler.x_present, proposed[i]) # w(yⱼ, x) = π(yⱼ)q(yⱼ|xₜ)λ(yⱼ, x) w_x2y[i] = sampler.target(proposed[i]) * q_x2y * sampler.λ(q_x2y, q_y2x) end # w(y,xₜ)したがってyをサンプリング y = sample(proposed, Weights(w_x2y)) # Σⱼ w(xₜ, yⱼ) for i in 1:(n_candidate-1) # q(xⱼ|y)でxⱼをサンプリング proposed_r[i] = sampler.propose_next(y) # q(xⱼ|y) と q(y|xⱼ) を計算 q_y2x = sampler.q(proposed_r[i], y) q_x2y = sampler.q(y, proposed_r[i]) # w(xⱼ, y) = π(xⱼ)q(xⱼ|y)λ(xⱼ, y) w_y2x[i] = sampler.target(proposed_r[i]) * q_y2x * sampler.λ(q_y2x, q_x2y) end # set xₘ=x proposed_r[end] = sampler.x_present # q(xₘ|y) と q(y|xₘ) を計算 q_y2x = sampler.q(proposed_r[end], y) q_x2y = sampler.q(y, proposed_r[end]) w_y2x[end] = sampler.target(proposed_r[end]) * q_y2x * sampler.λ(q_y2x, q_x2y) # αの計算 odds = sum(w_x2y)/sum(w_y2x) # 受容するか棄却するか決定 if rand(sampler.random_state) < odds sampler.x_present = y accept = 1 return y, accept else return sampler.x_present, accept end end function burn_in(sampler, n_candidate) # 焼きなまし期間 arr = zeros(sampler.burn_in) accept = zeros(sampler.burn_in) for i in 1:sampler.burn_in arr[i], accept[i] = next(sampler, n_candidate) sampler.count += 1 end return arr end function run(sampler, n_steps, n_candidate) #好きなだけ回す arr = zeros(n_steps) accept = zeros(n_steps) for i in 1:n_steps arr[i], accept[i] = next(sampler, n_candidate) sampler.count += 1 end return arr, accept end end
上で作った MTMを動かすところ。
const σ_proposal = 5*sqrt(2) const seed = 0 rng = srand(seed) # 乱数シードの固定とジェネレータの取得 gaussian(xₜ) = Normal(xₜ, σ_proposal^2) # 推移核の分布関数 generate(xₜ) = rand(gaussian(xₜ)) # 推移核の分布関数からサンプリングする関数 uniform = Uniform(-10, 10) # 一様分布関数 inits = rand(uniform) # 初期値の生成 target = target # 意味ないけど、見た目のため q(y, xₜ) = pdf(gaussian(xₜ), y) λ(x, y) = ((x+y)/2)^(-1) # MTM-inv scheme（下の参考文献参照） const burn_in = 1000000 const num_sample = 100000 cnt = 0 sampler_MTM = Sampler(inits, target, generate, q, λ, burn_in, cnt, rng)
MTMの実行
ベイズ計算統計学（統計解析スタンダード）p93 例3.7の例と同じように、候補数を1, 5, 10と変えてみた。
候補数1は単純なメトロポリス・ヘイスティングス法になります。
焼きなまし
const n_candidate=10 # MTMでつかう候補の数（ここでは10） samples_burn_in = MTM.burn_in(sampler_MTM, n_candidate)
本番。ステップごとのサンプルと受容フラグが返ってくる。
samples, accept = MTM.run(sampler_MTM, num_sample, n_candidate)
結果の確認
n_slice = 3000 lags = collect(1:200) auto_correlation = autocor(samples[end-n_slice:end], lags) layout = plt.@layout [a b; c] p1 = plt.bar(auto_correlation, title="Auto correlation", lab="", xlab="Lag", ylab="Auto\ncorrelation (A.U.)", xlim=(-10,200), ylim=(-0.1,1)) p2 = plt.plot(samples[end-n_slice:end], title="Trajectory", lab="") p3 = plt.histogram(samples , title="Target distribution", lab="", norm=true, bins=128) p3 = plt.plot!(x,y, color="magenta", lab="") p=plt.plot(p1, p2, p3, layout=layout) print(sum(accept[end-n_slice:end]) / length(accept[end-n_slice:end])) # 採択率を算出
できました。
候補を増やすと、目標分布の2つの山を頻繁に行き来しているのがわかります。
採択確率も候補を増やしたほうが良くなっています。
ただ、候補の数分計算量が増えるので注意ですね。
まとめ
実装は簡単だったが、はてなブログに書くのが非常に疲れt
間違いがあったらご指摘ください。
ちなみに以前書いたメトロポリス・ヘイスティングスは少し間違ってます（笑）
参考文献
Pandolfi, S., Bartolucci, F. & Friel, N. A generalized multiple-try version of the Reversible Jump algorithm. Comput. Stat. Data Anal. 72, 298–314 (2014).
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 期間が短いときに途中から再計算させたいので。
Julia tips #4: Cairoのインストールに失敗する
環境
- Windows 10 HOME
- Julia 0.6.2
JuliaでグラフをプロットするライブラリとしてはGadlyが有名？ですが、プロットをpng形式で保存するためにはCairoというパッケージが必要です。
そこでPkg.add("Cairo")したらインストールに失敗したというお話。
解決法も書いておきます。
エラーメッセージ
Pkg.add("Cairo")したら何やらダウンロードエラーが...
WARNING: Unknown download failure, error code: 2148270086 WARNING: Retry 1/5 downloading: https://cache.julialang.org/http://download.opensuse.org/repositories/windows:/mingw:/win64/openSUSE_Leap_42.2/noarch/mingw64-zlib1-1.2.8-8.15.noarch.rpm WARNING: Unknown download failure, error code: 2148270086 WARNING: Retry 2/5 downloading: https://cache.julialang.org/http://download.opensuse.org/repositories/windows:/mingw:/win64/openSUSE_Leap_42.2/noarch/mingw64-zlib1-1.2.8-8.15.noarch.rpm WARNING: Unknown download failure, error code: 2148270086 WARNING: Retry 3/5 downloading: https://cache.julialang.org/http://download.opensuse.org/repositories/windows:/mingw:/win64/openSUSE_Leap_42.2/noarch/mingw64-zlib1-1.2.8-8.15.noarch.rpm WARNING: Unknown download failure, error code: 2148270086 WARNING: Retry 4/5 downloading: https://cache.julialang.org/http://download.opensuse.org/repositories/windows:/mingw:/win64/openSUSE_Leap_42.2/noarch/mingw64-zlib1-1.2.8-8.15.noarch.rpm WARNING: Unknown download failure, error code: 2148270086 WARNING: Retry 5/5 downloading: https://cache.julialang.org/http://download.opensuse.org/repositories/windows:/mingw:/win64/openSUSE_Leap_42.2/noarch/mingw64-zlib1-1.2.8-8.15.noarch.rpm INFO: try running WinRPM.update() and retrying the install
大事なのは最後の1行のメッセージでした。
この通りにやってみるとインストールできました。
julia> using WinRPM
julia> WinRPM.update()
最後にさっきインストールに失敗したCairoをビルドし直します。
Pkg.build(Cairo)
これでうまくインストールできました。
Julia tips #3: ODE solverの速度比較 with Python
JuliaではSundials.jlをインストールするだけで、SundialsのODEソルバ、CVODEが使えます！！（マニアックすぎ？）
未知システムのパラメータ推定をする際、数値積分の区間をうまいこと刻まないと値が飛びすぎて計算できなかったり、刻む必要のないところで刻んでしまって計算不能に陥ることがあります。
そこでCVODEなどのAdaptive stepなODEソルバが必要なわけです。
何個か種類がありますが、有名どころでは
- LSODA
- CVODE
が2強といったところで、Pythonでよく使われるodeintもadaptive stepです（内部的にはLSODEとかいうやつらしいけど...良く知りません）。
じゃあ、JuliaでCVODEを呼び出して使うのと、Pythonでodeint使うのはどっちがどんくらい速いのよ？と思ったので比べてみました。
なお、Juliaのソルバについては比較してくれている方がいました。
CVODEの圧勝です。
Juliaのコードもこちらから借りてきて、Pythonと比較したいと思います。
実装
Julia（上記リンク先から部分的に拝借）
using Sundials # 非線形パラメータ mu = 0.2 # 関数定義 function vdp_sun(t, u, du) # dx/dt du[1] = u[2] # dy/dt du[2] = mu * (1.0 - u[1]^2.0) * u[2] - u[1] du end # パラメータ設定 u0 = [0.2; 0.2] t = [0.0:0.1:10.0;]
Python (3.6)
import numpy as np from scipy.integrate import odeint # 非線形パラメータ mu = 0.2 # 関数定義 def vdp_sun(u, t): # dx/dt dxdt = u[1] # dy/dt dydt = mu * (1.0 - u[0]**2.0) * u[1] - u[0] return [dxdt, dydt] # パラメータ設定 u0 = [0.2, 0.2] t = np.arange(0,10.1,0.1) # メインの計算 res = odeint(vdp_sun, u0, t)
速度比較
Juliaは1ループごとに出力、Pythonのtimeitではループにかかった平均時間を出力していますが、大雑把な比較はこれで良いでしょう。
見やすさのため10回だけ回しています。
Julia
res = Sundials.cvode(vdp_sun, u0, t); for i in 1:10 @time res = Sundials.cvode(vdp_sun, u0, t); end
Python (3.6)
%timeit -n 10 odeint(vdp_sun, u0, t)
結果
Julia
0.000206 seconds (1.81 k allocations: 43.734 KiB) 0.000132 seconds (1.81 k allocations: 43.734 KiB) 0.000130 seconds (1.81 k allocations: 43.734 KiB) 0.000117 seconds (1.81 k allocations: 43.734 KiB) 0.000118 seconds (1.81 k allocations: 43.734 KiB) 0.000115 seconds (1.81 k allocations: 43.734 KiB) 0.000115 seconds (1.81 k allocations: 43.734 KiB) 0.000115 seconds (1.81 k allocations: 43.734 KiB) 0.000116 seconds (1.81 k allocations: 43.734 KiB) 0.000124 seconds (1.81 k allocations: 43.734 KiB)
Python (3.6)
716 µs ± 5.63 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
まとめると以下のようになります。
- Julia: 128.8 µs ± 26.43 µs
- Python: 716 µs ± 5.63 µs
Juliaの方が5倍ほど速い！
シミュレーションを何千回、何万回と繰り返すものにとっては5倍の高速化は結構大きいと思います。
予想ではもっと変わるかなと思っていたのですがこんなものでした。
Pythonを高速化
Pythonを高速化したら、Sundialsに迫るのでは？？
と考えたので、numbaで高速化したPythonで戦ってみた。
odeintに渡す微分方程式をjit+型指定すると速くなる。
from scipy.integrate import odeint import numpy as np from numba import jit # 型指定jit, tが逐次的に渡されることに注意 @jit("float64[:](float64[:], float64)", nopython=True) def vdp_sun_jit(u, t): # dx/dt dxdt = u[1] # dy/dt dydt = mu * (1.0 - u[0]**2.0) * u[1] - u[0] return np.array([dxdt, dydt]).astype(np.float64) # パラメータ設定 u0 = np.array([0.2, 0.2]).astype(np.float64) # 渡す型に合わせる t = np.arange(0,10.1,0.1).astype(np.float64) # 渡す型に合わせる # メインの計算 res = odeint(vdp_sun, u0, t) %timeit -n 10 odeint(vdp_sun_jit, u0, t)
結果
387 µs ± 11.6 µs
うーん、速くはなりましたが2倍行かないくらいでしたね。
まとめ
やっぱSundials凄い。
数値計算ではJuliaを使おう！
Julia tips #2: ファイル検索
スタックオーバフローより。
「組み込み関数はないけどワンライナーで書けるよ」
とのこと。
関数化したほうが便利ですね。以下コードです。
searchdir(path,key) = filter(x->contains(x,key), readdir(path))
readdir (path)
渡したpathの中にあるディレクトリ、ファイルを文字列として出力してくれる
contains (x, key)
文字列 x に 文字列 key が含まれているかtrue or falseで返してくれる
filter (function, collection)
collectionを受け取る関数functionがtrueを返したら、そのcollectionのコピーを返す
例）
a = 1:10 julia> filter(isodd, a) 1 3 5 7 9
