あいかわらずWindows環境で苦行を続けております。私は元気です。

さて、ベイズ統計を実施するにあたって、環境を色々模索しておりましたが、 Windows + Python + PyStanは茨の道すぎて困難を極め、PyMC3も色々環境を壊してくれてありがとうなことから、 黒木先生おすすめのJuliaで着手。 「変わらんな 今度はJuliaに目を付けたか」という感じですね。すみません。

使うのはRでもStanでもないですが、ベイズ統計の考え方自体は教科書として松浦先生のご本を参考にします。

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

• 作者: 松浦健太郎, 石田基広
• 出版社/メーカー: 共立出版
• 発売日: 2016/10/25

では環境づくりからスタート

環境構築(インストール)

• windows10

Python系環境について

Python系の環境はAnaconda 3.6.0::Anaconda 4.3.0(64bit)使用しています。 Anaconda 5.0.1はmatplotlib.pyplotのimportでクラッシュして、なんだったんだ、あれは。

JULIAの環境づくり

黒木 玄先生のWindowsでJuliaをJupyterから使う をよく参考にさせていただいて、JupyterでJuliaを使えるようにします。 Juliaはインストール済みとして進めます。

1. 環境変数の設定

金言に従ってPKGDIRをしっかり設定。

JULIA_HOME=D:\Users\UserName\AppData\Local\Julia-0.6.1\binJULIA_PKGDIR=C:\Users\UserName\.julia\v0.6

JULIA_PKGDIRは、下記メソッドで調べられます.

Pkg.dir()

※ 実行はJuliaのターミナルで(スタートメニュー>Julia)

2. IJuliaの導入～jupyter notebookへの設定

工夫なくIJulia入れるだけ。Windowsのパスを設定するときは生文字リテラルraw"文字列"が便利そう。

Pkg.add("IJulia")ENV["PYTHON"]=raw"C:\Users\UserName\Anaconda3\python.exe"
Pkg.add("PyPlot")

あとは、Jupyter notebookを起こしてjuliaカーネルのノートが使えるか確認するだけ。

jupyter notebook

しばらくすると、jupyterが起動してブラウザでアクセスができます。

3. 一応確認

Notebook上で一応確認。
Plots/GR: グラフ package のおすすめ · julia について

using Plots
gr()
plot(randn(100,3))

【補足】便利かと思って足しておいたパッケージのご紹介

Pkg.add("Plots")
Pkg.add("GR")
Pkg.add("Pandas")

参考:RユーザーのためのJulia100問100答 - りんごがでている

線形回帰

アヒル本ことStanとRでベイズ統計モデリング (Wonderful R)の4章で取り上げられている例をもとに、線形回帰。 分析対象のデータはGitHubに公開されているものを使います→RStanBook/data-salary.txt

モデル

見本のStanでのモデリング(model4-5.stan)

data {
    int N;
    real X[N];
    real Y[N];
}

paramters {
    real a;
    real b;
    real sigma;
}

model{
    for (n in 1:N) {
        Y[n] ~ normal(a + b*X[n], sigma);
    }
}

これをmamba.jlで書くと

model = Model(
    y = Stochastic(1,
        (a, b, x, s2) -> MvNormal(a + b*x, sqrt(s2)),
        false
    ),
    b = Stochastic(() -> Normal(0, 100.0)),
    a = Stochastic(() -> Normal(0, 100.0)),
    s2 = Stochastic(() -> InverseGamma(0.01,0.01))
)

sigmaについては $ s^2 = sigma*sigma $ として代用。Stanとは異なり、無情報事前分布の指定が必要で、アヒル本(p.31)で触れられています。 本書によるとStanの場合は十分広い範囲の一様分布とのことですが、 gitの記事Prior Choice Recommendations · stan-dev/stan Wiki には、事前分布の選び方に議論があり、Stanのリポジトリも軽く探してみましたが良くわかりません。 宿題ということにして、次に進みます。

a, bには正規分布。s2には逆ガンマ関数を利用。定数の選び方もちょっと苦しいので詰められるときついっす。逆ガンマ関数については、正規分布とセットで使うことが多いとかとか→ベイズ統計

上記モデルでは、 $ y = a + bx + \epsilon $ の表現で基本モデルを定義しましたが、線形回帰 - Wikipedia でも言及されている通り、ベクトル・行列記法を用いて、 $ Y = X \beta + \epsilon $ とも表せます。 この場合、mambaのtutorialの通り、 こちらのモデルになります。

model = Model(
    y = Stochastic(1,
        (mu, s2) ->  MvNormal(mu, sqrt(s2)),
        false
    ),
    mu = Logical(1,
        (xmat, beta) -> xmat * beta,
        false
    ),
    beta = Stochastic(1,
        () -> MvNormal(0, 100)
    ),
    s2 = Stochastic(
        () -> InverseGamma(0.01, 0.01)
    )
)

閑話休題で、$ y = a + bx + \epsilon $ の形式のモデルで続けます。

Samplingスキームの設定

a,bについては、RStanと同様にNUTSを使います.

s2については、Tutorial — Mamba.jl 0.10.1 documentationの Sampling Scheme2と同様にすべてNUTSを試してみたけど、うまくいかずScheme1を踏襲。

# Data Conversion, PandasデータフレームとしてData読み込み済みなので、mamba向けに変換
Dat = Dict{Symbol, Any}(
    :x => Array(Data[:X]),
    :y => Array(Data[:Y])
)
# Initial Values
inits = [
    Dict{Symbol, Any}(
        :y => Dat[:y],
        :a => rand(Normal(0, 1)),
        :b => rand(Normal(0, 1)),
        :s2 => rand(Gamma(1,1))
    ) for i in 1:4
]
## MCMC Simulations
sim = mcmc(model, Dat, inits, 2000, burnin=1000, thin=1, chains=4)
Mamba.describe(sim)

MCMC実行

mcmcのパラメータ設定はRStanのデフォルトをできるだけmimicするようにしたつもり

sim = mcmc(model, Dat, inits, 2000, burnin=1000, thin=1, chains=4)

burnin

RStanではwarmupと呼ばれている

第4引数

上記では2000. Iteration回数。

実行結果の転記

s2の結果が収束せず、アヒル本より小さめ(sigma 85前後⇔s2 7225程度に対して、1103)になり、どうしたものかという状況。

	Mean	SD	Naive SE	MCSE	ESS	2.50%	25.00%	50.00%	75.00%	97.50%
a	-113.3	29.0	0.5	1.5	353.5	-171.2	-131.4	-112.9	-95.9	-52.1
b	21.8	0.6	0.0	0.0	390.8	20.4	21.4	21.8	22.2	23.1
s2	1103.6	112.3	1.8	17.8	40.0	830.9	1047.4	1102.1	1182.5	1279.1

プロット

s2…悪いですね… うーん？？a,bの信頼区間も狭いような？？

追記:mamba.jlでも上手く行った(2018-01-03)

続く数日のポスト（MCMCと単回帰 - 夜間飛行）で勉強しながら、しばらく検証に時間をかけていましたが、mamba.jlでも良い結果が出るようになってきました。Stan モデリング言語: ユーザーガイド・リファレンスマニュアル9章2節の事前分布に関する項目が要因にあたっていたように思います。

BUGSの例題(Lunn et al., 2012)は全般に, スケールについてのデフォルトの事前分布を下のようにしていますが, これは使わないようにしましょう. $$ \sigma^2∼InvGamma(0.001,0.001) $$ このような事前分布は, 妥当な事後分布の値の外側に, あまりに大きく集中した確率の山があります. 正規分布につける対称的で広い事前分布とは異なり, これにより事後分布をゆがめる深刻な影響が発生する可能性があります. 例題と議論はGelman (2006)を参照してください.

手直し

model = Model(
    y = Stochastic(1,
        (a, b, x, sigma) -> MvNormal(a + b*x, sigma),
        false
    ),
    b = Stochastic(() -> Normal(0, 1000)),
    a = Stochastic(() -> Normal(0, 1000)),
    sigma = Stochastic(() -> Uniform(0, 1e6))
)

事前分布変更後の結果

実質1000サンプリング程度だと目が粗いですが、なかなか良い結果に落ち着いてきました

Empirical Posterior Estimates:
          Mean        SD       Naive SE     MCSE       ESS
sigma   84.498765 15.5085775 0.490424280 0.68927615 506.24112
    b   21.955192  1.6355915 0.051721944 0.17563254  86.72408
    a -120.771614 72.9060900 2.305492998 7.94019602  84.30730

サンプリング増し

上手く行ったので50,000サンプル取ります

きれいですね

手直しNUTS版はこちらに置きました→Mamba.jl_Practice/NUTS-mamba-Retry.ipynb