StanとJuliaでベイズ統計モデリング

確率的プログラミング言語を学ぶ機運を感じたので勉強している。
以下の本を元に勉強しているが、写経するのも面白くないのでJuliaでやっていく。
環境はUbuntu16.04、Julia0.6.2。

StanとRでベイズ統計モデリング (Wonderful R)

インストール

Stan.jlはCmdStanを呼び出すため、CmdStanとStan.jlの両方のインストールが必要である。
Stan.jlのインストール時にいい感じにやってくれんかなと思ったけどそんなことはないし、やられたらそれはそれでうざい。

CmdStanのインストール

aptとかにもないので公式ページから落としてくる必要がある。
Stan - CmdStan

私は~/bin以下に解凍した。
次に.zshrcに以下を追記。

export CMDSTAN_HOME=/home/youraccount/bin/cmdstan-2.17.1

いつも通りsource .zshrcして完了。

Stan.jlのインストール

CmdStanが設定されていればPkg.addするだけである。

Pkg.add("Stan.jl")

Juliaのlm関数で推定

上記書籍でも一旦、Rのlmを使用して回帰させているので、こちらも一旦それに習う。
GLM.jlを使用し同じことを行う。

using CSV
using GLM
using DataFrames

data = CSV.read("data-salary.txt")
profit = lm(@formula(Y ~ X), data)
new_x = DataFrame(X=23:60)
                                                                                      
@show predict(profit, new_x) 

ちなみに予測区間を出すにはpredict(profit, new_x, :confint)とやれば良さそうにREADMEには書いてあるが、
DataFrameからMatrixの変換はStatsModelsで実装されており、StatsModels側の実装がSymbolを受け取るようなI/F
ではないため、できない。READMEにはかいてあるのに。
issueをみると1年くらい前に治そうとしている形跡はある。
私も直してやろうかと思ったらそもそもテストが落ちているので、まずそこを直そうとなり、これがそもそもこの仕様やめね？という話に飛び火してつらみを出してる。
巷でJulia最高と言っている人間はこういう辛い面を全く言わないので多分違うJuliaを使っている気がする。

Stanで実装、Juliaからの実行方法

上記書籍のmodule4-5.stanを写経する。
下からダウンロードしても良い。
RStanBook/model4-5.stan at master · MatsuuraKentaro/RStanBook · GitHub

後はStan.jlからモデルを呼び出す。

using CSV
using Stan
d = CSV.read("data-salary.txt")
data = [Dict("N"=>size(d, 1), "X"=>d[:X], "Y"=>d[:Y])]
f = open("model4-5.stan", "r")
model = readstring(f)
close(f)
stanmodel = Stanmodel(model = model, random=Stan.Random(1234))
rc, sim = stan(stanmodel, data)

実行結果であるsimはMambaというMCMCのライブラリで扱われるデータ型になる。

実行結果を見る。

using Mamba
describe(sim)

最後にグラフに表示する。
MambaはGadflyでグラフ表示することを意識して作っているが、
Plots(GRBackend)が好き(というかそれ以外であまり上手く行ったことがない)なため、これで描画する。
今回もGadflyで描画しようとしてクラッシュしたので逃げたとも言う。

Traceを表示する。
といっても表示要素が多いためどれを表示したいかを選ぶ。

sim.names

表示させたいものが決まったらグラフに出す。

using Plots
gr()

nrows, nvars, nchains = size(sim.value)
i = 9 # bを表示
Plots.plot(1:nrows, sim.value[:, i, :], title=sim.name[i])

bが21付近をうろうろしているグラフが得られ、1年ごとに21万も昇給するわけねぇだろという気持ちになったので終わります。