MCMCと単回帰

Posted on 12/28/2017
このエントリーをはてなブックマークに追加

問題のおさらい

アヒル本の第4章のデータRStanBook/data-salary.txt を教材に、MCMCのアルゴリズムの勉強を進めています。ことの発端調査着手の経緯は過去ページをご覧になっていただくとして、ランダムウォーク・メトロポリス・ヘイスティングス法(Random Walk MH法)とハミルトニアン・モンテカルロ法(HMC法)に自分で組みながら単回帰に取り組みましたので その内容をまとめます。

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

 

単回帰問題の内容

架空データはB社の社員の年齢X(歳)と年収Y(万円)のデータで、20人分のサンプルがある。このサンプル数20という絶妙なサンプル数の単回帰を行い、予測に使おうという問題です。このうちの予測ではなく、単回帰までのところで、色々やっています。

データ: RStanBook/data-salary.txt

このデータは、$ Y_i \sim \mathcal{N}(a + b \cdot X_i, \sigma^2) $ で表現されるモデルにおいて、最尤法では $ \hat{a} = -119.697 $ , $ \hat{b} = 21.904 $ , $ \hat{\sigma} = 79.099 $ と推定されます。

視覚化のため、プロットもアヒル本の図4.3から拝借してきました(githubへリンク)。このような見た目のデータで、直線は最尤推定によるもの。灰色は信頼性区間。

最尤法に関連したおさらい

傾き $ b $ , 切片 $ a $ に関して、$ n $ をサンプル数とおいたとき、それぞれ、$ \hat{a} \sim \mathcal{N}(a, \sigma^2 \{\frac{1}{n} + \frac{\bar{X}^2}{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2} \}) $ , $ \hat{b} \sim \mathcal{N}(b, \frac{\sigma^2}{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2 }) $ , $ \hat{\sigma}^2 \sim \chi_{n-2}^2 \frac{\sigma^2}{n-2} $ ですから、 それぞれの信頼区間は以下の通り。 $$ a \in \left[\hat{a}-t_{\alpha/2,n-2} \hat{\sigma} \sqrt{\frac{1}{n} + \frac{\bar{X}^2}{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2}}, \hat{a}+ t_{\alpha/2,n-2} \hat{\sigma} \sqrt{\frac{1}{n} + \frac{\bar{X}^2}{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2}} \right] $$ $$ b \in \left[\hat{b}-t_{\alpha/2,n-2} \frac{\hat{\sigma}}{\sqrt{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2 }}, \hat{b}+ t_{\alpha/2,n-2} \frac{\hat{\sigma}}{\sqrt{\sum_{ i = 1 }^{ n } (X_n - \bar{X})^2 }} \right] $$ $$ \sigma^2 \in \left[\frac{(n-2)\hat{\sigma}}{\chi^2_{1-\alpha/2}}, \frac{(n-2)\hat{\sigma}}{\chi^2_{\alpha/2}} \right] $$

値を代入して95%信頼区間を計算します。$$ a \in \left[-262.8, 23.5 \right] $$ $$ b \in \left[18.7, 25.1 \right] $$ $$ \sigma^2 \in \left[3572, 13683 \right] $$ $$ \sigma \in \left[59.7, 117.0 \right] $$

参考

証明や理論的背景は、ガウス=マルコフの定理 あるいは赤本の第13章を参照お願いいたします。

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

 

統計値まとめ(Ground Truth)

以上をMCMC結果のGround Truthとして扱いますので、ここで一旦まとめ (信頼区間と確信区間を同じに扱っていますが、単回帰では数値的に同じになるはず)

パラメータ名 点推定 95%信頼区間下限値 95%信頼区間上限値 標準誤差
a(切片) -119.7 -262.8 23.5 68.15
b(傾き) 21.9 18.7 25.1 1.52
sigma(残差偏差) 79.1 59.7 117.0  

環境

シミュレーションで使った環境のご紹介。

  • Windows 10 - 64bit
  • Julia 0.6.1
  • mamba 0.11.1
  • R 3.4.3
  • RTools34
  • RStan Version 2.16.2

対数尤度関数(Likelihood function)

対数尤度関数 $ \mathcal{L} (y|a, b, sigma) $ は、こちら。$$ \mathcal{L} (y|a, b, \sigma)=log(\prod_{k=1}^n \mathcal{N}(y_k | a+b \cdot x_k, \sigma^2)) \\= - \frac{n}{2} log(2 \pi) - n \cdot log(\sigma) - \frac{\sum_{k=1}^n (y_k - a - b \cdot x_k)^2}{2 \sigma^2} $$

Distributions.jlを使って書きます。

using Distributions
function likelihood(a, b, sigma)
    y_pred = a + b*x
    likelihoods = [logpdf(Normal(aa, sigma), bb) for (aa,bb) in zip(y_pred,y)]
    l_sum = sum(likelihoods)
    return l_sum
end

通常MCMCによるフィッティングでは実施しないかもしれませんが、今回は収束値が最尤法で求められるので、点推定の値を中心に、対数尤度の分布をプロットして分布を確認します。

using Plots
pyplot()
# 点推定の値(を整数に丸めた)
a = -120.0
b = 22.0
sigma = 79.0

# 分布をみる範囲
a_list = collect(-262:24)
b_list = collect(18:0.1:25)
sigma_list = collect(59:117)

ab_list = [likelihood(tmp_a,tmp_b,sigma) for tmp_a in a_list, tmp_b in b_list]
ac_list = [likelihood(tmp_a,b,    tmp_sigma) for tmp_a in a_list, tmp_sigma in sigma_list]
cb_list = [likelihood(a, tmp_b,tmp_sigma) for tmp_sigma in sigma_list, tmp_b in b_list]

levels = collect(-160:2.5:-115)

p1 = contour(b_list,a_list,ab_list,levels=levels, xlabel="b", ylabel="a")
p2 = contour(sigma_list,a_list,ac_list, levels=levels, xlabel="sigma", ylabel="a")
p3 = contour(b_list,sigma_list,cb_list, levels=levels, xlabel="b", ylabel="sigma")

plot(p1,p2,p3)

この通り、b(傾き)に対して敏感で、次点でa(切片),最後にsigmaという様子が見られます。後述する通り色々試してみた中ではRandom-Walk Metropolis Hasting法やHMC法では、 この感度が収束に影響を及ぼしているようでした。MCMCにも機械学習に於ける正規化の考え方が必要ということです。

 

ランダムウォーク=メトロポリス・ヘイスティングス法

アルゴリズムの詳細は触れなくても大丈夫ですよね。おさらい用にリンク貼っておきますね。

Github → Mamba.jl_Practice/Random-walk MH.ipynb at 13b26fe1458a3f0e41b880e25650fcadecbbee5b · Chachay/Mamba.jl_Practice

提案分布

この手法による唯一のチューニング対象(?)はランダムウォークの仕方。酔歩に性格の良い正規分布を使うことを前提にして、分散がチューニングの対象です。 直感的にコーシー・シュワルツの不等式的な発想で各パラメータの事後分布の分散の比率に提案分布が近いとサンプリング効率が良いような気がしています。

function proposal(a, b, sigma)
  return [rand(Normal(mu, sd)) for (mu, sd) in zip([a, b, sigma],[80,1,8])]
end

試しに提案分布のaの分散をいじります。

どうです。aの標準誤差65に近い80が一番良い雰囲気出していませんか?aの分散を1にしたものと60にしたもので、3000サンプル分のランダムウォークの違いを見てみます。

分散=1

分散=60

分散60のほうが心なしか良さそうですよね(点推定値への収束が悪いなぁ)

さらに提案分布を多変量正規分布にすると、より一層良くなりそうな気がしていますが、色々やっておいてなんですが、事前にパラメータの分布がわかっているのも興ざめなので、ここ迄にしておきます。

参考:How to derive variance-covariance matrix of coefficients in linear regression - Cross Validated

メトロポリス・ヘイスティングス法の結果

チューニングしてみましたが(a, b)=(-70, 20.9)など、aが心持ち原点に寄ってくるところが解せないです。初期値を最尤法による点推定の値にしても、aが-70前後に寄ってきてしまうので頭抱えています。 実行のたびに-120を中心に-190~-70でブレるというのならわかるけど、毎回原点側に戻って来てしまうのは、 どういうことなの!?

aが原点に寄る…はっ、桁落ちでは?!

ということで、採択率の計算方法を少し変えました。$$ \alpha = exp(log(p^*) - log(p_t)) $$

を$$ \alpha = \frac{exp(log(p^*))}{exp(log(p_t))} $$ に変えました。

Before

tmp1 = posterior(proposal_draw[1],proposal_draw[2],proposal_draw[3])
tmp2 = posterior(chain[i][1], chain[i][2], chain[i][3])
alpha = min(1,exp(tmp1-tmp2))

After

tmp1 = posterior(proposal_draw[1],proposal_draw[2],proposal_draw[3])
tmp2 = posterior(chain[i][1], chain[i][2], chain[i][3])
alpha = min(1,exp(tmp1)/exp(tmp2))

何度か計算させてみても、大体、最尤法の点推定を中心にばらつきます。難点は対数尤度(負の値)が小さすぎるとゼロ割相当になってしまうことでしょうか。

数字で見て見ると $ (a, b, \sigma) = (-122.4, 22.0, 82.2) $ となり、良い!良いです!

mamba.jlでaが小さめに出る問題、桁落ちの可能性もありますね。実際mamba.jlのRandom-Walk MH法では、採択率の計算がこうなっています。きっと同じ症状が出ると思います。

他にもみていくとmamba.jlのHMCにも同じような箇所があり、少し心配になります。NUTSの実装ではHMCを使っていないように見えますので、杞憂かなぁ。

しかし、このような実装がNUTSにも…

不審ですが、これを確認するためHMCに行ってみましょう。

ハミルトニアンモンテカルロ法(HMC)

アルゴリズムについては、ここで、わざわざ書くほどでもないほど、色々な方が記事書いてますし、Stanのマニュアルも充実しています。

短く言うと $ P $ を運動量に見立てて、対数尤度関数をポテンシャル場に使いながら$$ H([a, b, \sigma], [P_a, P_b, P_{\sigma}]) = - likelihood(a, b, \sigma) + \frac{P_{a}^2+P_{b}^2+P_{\sigma}^2}{2} $$ で定義される「ハミルトニアン」を利用して無駄な酔歩を抑えようという感じの発想のアルゴリズムですね。 焼きなまし法みたいに、Iterationのたびにランダムな運動量(≒初速)を与えた玉を少し転がしてやってからサンプリングを繰り返し、 最終的に谷底の中心で跳ね回るようにしていく動作をします。

その効果がこちらのような素早い谷底への到達!!

ポテンシャルエネルギー $ U(a, b, \sigma) = -likelihood(a, b, \sigma) $ 、 運動エネルギー $ K(P_a, P_b, P_{\sigma}) = \frac{P_{a}^2+P_{b}^2+P_{\sigma}^2}{2} $ と定義すると、ハミルトニアン $ H = U + K $ になりますね。 これを手を離した振り子の要領で、系のエネルギーを保存しながら、動いてもらおうと思うと、 ハミルトン力学でいうところの正準方程式を数値解法で解いていきます。 これには、シンプレクティックな数値積分法リープ・フロッグ法を使います。

使わないと面倒なことになりそう。

ということでリープフロッグ法で対数尤度関数の偏微分を使いますので、計算します。$$ \mathcal{L} (y|a, b, \sigma)= - \frac{n}{2} log(2 \pi) - n \cdot log(\sigma) - \frac{\sum_{k=1}^n (y_k - a - b \cdot x_k)^2}{2 \sigma^2} $$ $$ \frac{\partial \mathcal{L} (y|a, b, \sigma)}{\partial a}= \frac{\sum_{k=1}^n (y_k - a - b \cdot x_k)}{\sigma^2} $$ $$ \frac{\partial \mathcal{L} (y|a, b, \sigma)}{\partial b}= \frac{\sum_{k=1}^n x_k (y_k - a - b \cdot x_k)}{\sigma^2} $$ $$ \frac{\partial \mathcal{L} (y|a, b, \sigma)}{\partial \sigma}= -\frac{n}{\sigma} + \frac{\sum_{k=1}^n (y_k - a - b \cdot x_k)^2}{\sigma^3} $$

function dlikelihood(a, b, sigma)
  y_pred = a + b*x
  
  grad_a    = sum(y - y_pred)/sigma^2
  grad_b    = dot(x, y - y_pred)/sigma^2
  grad_sigma= -1.0/sigma * length(x) + sum(abs2,y - y_pred)/sigma^3

  return [grad_a, grad_b, grad_sigma]
end

で、Leapfrogの実装は

function Leapfrog(theta, r, epsilon)
    r += 0.5 * epsilon * dlikelihood(theta[1], theta[2], theta[3])
    theta +=   epsilon * r
    r += 0.5 * epsilon * dlikelihood(theta[1], theta[2], theta[3])
    return theta, r
end

ポテンシャルエネルギーをlikelihoodの-1倍で定義したので、マイナスがキャンセルして"r-="ではなく"r+="になっています。あしからず。

Github →Mamba.jl_Practice/HMC-MassMatrix.ipynb

ステップ数 $ L = 100 $, サンプリング周期 $ \epsilon = 0.2 $ として点推定値付近を50,000 Iterationsほどウロウロさせると、良さそうな雰囲気。

最尤値(オレンジ実線)とサンプリングの比較

$ (mean(a), mean(b), mean(\sigma)) = (-118.0, 21.9, 85.4) $ で、桁落ちの影響もなく良さそうですね。mamba.jlの問題はなんだか、特定できず終いです。

ここまでの話は、この本でカバーされているようなので…ほしいけどゴクリ…。

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

 

mamba.jlとMCMCを巡る旅

Posted on 12/09/2017
このエントリーをはてなブックマークに追加

2017年Julia Advent Calendar 9日目の記事です。

Juliaとベイズ統計モデリング、mamba.jlに同時入門した先日。 先人の資料の助けもあってPkd.add("Mamba")に次いで、jupyterとijulia上でusing Mambaがスムーズに動いたとき、 それまでWindows上のPyStan(+Visual Studio)の構築などで散々イライラしていた私は、安堵感と共に、 「お、これは幸先良いで。締切まで1週間以上あるし、信頼性工学を切り口にLife Data Analysisの信頼性区間くらいは、 さっさと出して、CalenderではTest Designの話でも書くかな」なんて呑気なことを考えながら手始めに、 線形回帰でのTutorialに取り組んでいました。

そう、これが今回の始まりです…。そして、非常に恐縮ながら途中まで…という状況で、旅路の記録になっています。

旅の抜け道をご教示いただけるなら、それも助かります!!!

環境

まず、動作環境をご紹介します. Forkが使えないWindowsを使い続けているのはお約束です。

  • Windows 10 - 64bit
  • Julia 0.6.1
  • mamba 0.11.1

なぜWindowsにこだわるのか、なぜWindows subsystem for Linuxを使わないのか、そこは今日も触れないでおきましょう。

線形回帰のおさらい

線形回帰が収束しない。これだけである。

アヒル本の第4章のデータRStanBook/data-salary.txt では、$ y \sim Normal(a+bx, sigma) $ のモデルにおいて、aとbの事後平均は、それぞれ、-118, 22ほどです。 これは大体最尤法の結果と同じ。そして、sigmaは85程度。

これに対して、前回の記事では、a,b,sigma は、それぞれ、-113, 22, 33で、a,bは許せるもののsigmaはプロットを見ても収束していませんでした。

元のIteration数が2000, Burnin(warm-up)が1000と小さめでしたので(とはいうものの、この例題でアヒル本は2000程度で上記の結果)、Samplerの性能の違いなどへの想像をしながらIterationを8000, Burninを3000へ伸ばしました。

Iterations = 3001:8000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 5000

Gelman, Rubin, and Brooks Diagnostic:
            PSRF 97.5%
        s2 1.000 1.000
        b 1.001 1.003
        a 1.002 1.003
Multivariate 1.001   NaN

収束は、していそうですね…

プロットはそこそこ雰囲気よさそうです

しかし… describe(sim).

Empirical Posterior Estimates:
    Mean         SD         Naive SE       MCSE         ESS   
s2 6986.028670 2592.8494023 18.3342139498 29.570409665 5000.00000
b   21.064966    1.2315484  0.0087083624  0.038741347 1010.54010
a  -80.859994   53.7228503  0.3798779177  1.789205605  901.56614

あれ?a、小さくないですか?!a,b,sigma、それぞれ、-80, 21, 84ですよ?

真実を求めて

Rとの比較

Iterationや初期値、modelの定義をとっかえひっかえしたりするも、事態が好転しないうちに、ひょっとしてmamba.jlのIssueに…と思いはじめました。

調べます。

ありました。linear regression example not giving reasonable results · Issue #120 · brian-j-smith/Mamba.jl

しかも、今年の7月から、まだOpenです。

After some digging it looks like the NUTS epsilon value is not getting set correctly in some cases. I will investigate further.

初期値でなんとかなる?全然納得できません。だいたい、単回帰程度のモデルで初期値に過敏に反応するようでは実務上お話になりません。NUTSの実装で $ \epsilon $ が悪い?よくわからないっす。

ここで、とっかえひっかえした経験から代表的な結果を挙げます。「おかしい」と感じる数字に*を打ってあります。mamba大丈夫ですか…?

条件 a(切片) b(傾き) sigma(残差偏差)
R glm(Ground Truth扱い) -119.7 21.9 79.1
R Stan(細かい設定なし⇔[初期値不明; Sampler NUTS; Iteration 2000]) -119 22 85
mamba.jl - 初期値~ Normal(Ground Truth, 1.0); Sampler Scheme1; Iteration 2000 -80* 20 83
mamba.jl - 初期値 a,b~Normal(0, 1), sigma~Uniform(80,85); Sampler Scheme2; Iteration 8000 -81* 21 84
mamba.jl - はじめてのmamba -113 22 33*

Ground Truth扱いの値を出すRのコード

df<-read.csv("./data-salary.txt")
data <- list(N=nrow(df),X=df$X, Y=df$Y)
ans<-lm(data$Y~data$X)
ans
summary(ans)
sigma(ans)

Rも一応初期値いじってみました(やりかたあってます?)

fit <- stan(file="./stan45.stan",
            data=data,
            seed=1234,
            verbose=TRUE,
            init=function(){ list(a=rnorm(1,0,1), b=rnorm(1,0,1), sigma=rgamma(1,1))},
            chains=4, iter=2000, warmup=1000, thin=1
)
data {
    int N;
    real X[N];
    real Y[N];
}

parameters {
    real a;
    real b;
    real<lower=0> sigma;
}

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

mamba.jlでのModelとSamplerの定義(s2ではなくてsigmaにしました)

model = Model(
    y = Stochastic(1,
        (a, b, x, sigma) -> MvNormal(a + b*x, sigma),
        false
    ),
    b = Stochastic(() -> Normal(0, 100)),
    a = Stochastic(() -> Normal(0, 100)),
    sigma = Stochastic(() -> InverseGamma(0.01, 0.01))
)
scheme1 = [NUTS([:a,:b]),Slice(:sigma,10)]
scheme2 = [NUTS([:a,:b,:sigma])]

Tutorialの確認

そして、疑いの目でMambaのTutorial — Mamba.jl 0.10.1 documentationも見てみます。

条件 a(切片) b(傾き) sigma(残差偏差)
R glm(Ground Truth扱い) 0.6 0.8 0.73
mamba.jl Tutorial(Iteration 5000の後) 0.6 0.8 1.08(⇔s2=1.18)

s2はカイ二乗分布に従うので自由度5くらいの例だと言いにくいですが、こちらの残差分散も怪しそうです。私が使い方を理解していない…ということ以上のことが起きているのが、いよいよ濃厚に感じられてきました。 Mamba.jlのMCMCに何か起きていそうです。ライブラリの中に潜るしかなさそうな様子です。

次ステップの逡巡

ここで、私にはいくつかの選択肢がありました。カレンダー締め切りが迫る中、妻の協力を得つつ、検証時間を確保した私は考えます。

ベイズ統計初学者、MCMCに至っては入門すらしていない私は、ろくな初期化もせず、2000やそこいらのIterationでGLMに近い答えを出してしまうRStanがすごすぎるのか、mamba.jlの熟成が足りておらず、何らかの間違いを含んでいるのか見極めないと何も言えません。 RStanには、SamplerにNUTSを採用していることが書かれており、 バージョンもいくつかありそうなことがわかります。 NUTSのパラメータをRStan, Mambaで比べてみると、若干違うことも気になります。

algorithm One of sampling algorithms that are implemented in Stan. Current options are "NUTS" (No-U-Turn sampler, Hoffman and Gelman 2011, Betancourt 2017), "HMC" (static HMC), or "Fixed_param". The default and preferred algorithm is "NUTS".

サーベイを軽くしながら、思い浮かんだことリスト。

  • とりあえず結果を…型
    • 事前分布をもっと広げて、うすーくしたり、事後分布とほぼ1:1にしてしまったりする
    • mamba.jlの他のSamplerをチューンして、とりあえずGLMに近い結果が出るまで、ガチャガチャを回す
    • mamba.jlで$ sigma = \frac{10000}{1+exp(-dummySigma)} $ を用意して下限値を設定する. またはそれに準じたmodel定義を.
  • 動作状況確認型
    • StanのSampler実装を読み込んで、mamba.jlのNUTS.jlと比べる.
    • いや、そもそもSamplerじゃないじゃないかもしれないじゃないか。ひょっとして事後確率分布の計算とかで違ってる?
    • NUTS(No-U-Turn Sampler)を自分で組んでみて比較する
      • mamba.jlのユーザー定義Samplerに自分で実装
      • MCMCも含めてjuliaで書いてみる?
  • 基礎テスト実行型
    • 単回帰よりさらに単純な条件での評価

MCMCへダイブ

次第に苦しくなってきましたが、結局、NUTSを自分で組んで確認する方向を取ってみることにしました。CourseraでNg先生の機械学習で線形回帰の機械学習やNNを組んだように、 単回帰に特化すればMCMCとて意外と書けるかもしれません。

どこから始めるか

NUTS(No-U-Turn Sampler)は、HMC(Hamiltonian Monte Calro)を改良したもので、HMCは、Metropolis-Hasting Algorithmの提案関数を効率的にしたものということです。 幸いなことに、MH法はWikipediaもあるくらいですし、 ここから始めましょう。

ランダムウォークMH法

MH法の提案分布のうち、正規分布のような対称な分布を使ものを、Random-walk MH algorithm(ランダムウォークMH法)と呼ぶとのこと。日本語のWikipediaではメトロポリス・アルゴリズムと呼んでいるようです。 実装はMITの資料を見つけたのでこちらにならいます。

まず、概要は非常にスッキリ。こちらのランダムウォークの部分、そして、alphaの計算がNUTSでは複雑になっていくという心構えで見ていけば、よさそうな感触がつかめます。コメントを付けた通り、ランダムウォークMH方では、 提案分布が対称であるため、alphaの計算がシンプルにできます。

#Random-walk MH法
# モデルはy ~ N(a+b*x, sigma)
# a, b, sigmaの初期値; Chain 1つ分
Init = Real[-200, 1, 1] 
Iteration = 90000
BurnIn    = 10000
chain = Array[Init]

for i in 1:Iteration
    # 前回のa,b,sigmaの推定値からランダムウォークで次の推定値候補を引き出す
    # chain[i][1] : a
    # chain[i][2] : b
    # chain[i][3] : sigma
    # 戻り値は配列値で[aの候補, bの候補, sigmaの候補]
    proposal_draw   = proposal(chain[i][1], chain[i][2], chain[i][3])
    
    # 提案分布が対象であることから、MITの資料 数式(1)によってalphaの計算が非常に簡単になる
    # posteriorは目標分布の確率密度関数(の対数値)
    tmp1 = posterior(proposal_draw[1],proposal_draw[2],proposal_draw[3])
    tmp2 = posterior(chain[i][1], chain[i][2], chain[i][3])
    alpha = min(1,exp(tmp1-tmp2))
    
    # 採択するかどうか決める!
    u = rand(Uniform(0,1))
    if u < alpha
        push!(chain,proposal_draw)
    else
        push!(chain,chain[i])
    end
end

# BurnInは除外
chain = chain[BurnIn+1:end]
# 結果の取り出し
a = [chain[i][1] for i in 1:length(chain)]
b = [chain[i][2] for i in 1:length(chain)]
sigma = [chain[i][3] for i in 1:length(chain)]

println("MCMC complete")

初期値[a,b,sigma]=[-100, 1, 1]; Iteration = 200000; BurnIn = 90000など、いくつか特徴を見てみましたが、切片と残差偏差が、トレードオフのような関係にあってなかなかGround Truthに近づきません。 たとえば先述の条件では(a, b, sigma)=(-70, 20.9, 85.2)でした。 mamba.jlの方に近い!!

こちらのコードはGithubにてご紹介しています → Mamba.jl_Practice/Random-walk MH.ipynb at 4ce5b1fd3069cbf0ac6bc313489e62dad2ac0625 · Chachay/Mamba.jl_Practice

お詫びと感想

今回はMCMCおよびベイズ統計の勉強が足りず、みなさまにお見苦しい記事をお見せしてしまいました。

ただ、Random-walk MH法は是非動かしていただきたいのですが、Pythonでは考えられないほどのIteration速度で、ただ、ただ、驚くばかりです。これぞJuliaを使う素晴らしさでしょう。 しかしながら、PythonやRと違って、開発環境やパッケージの洗練、コミュニティはの厚みには、 もう少し時間がかかりそうなこともも感じました。 エラーの発生時も箇所の特定に少しコツがいりますし、mamba.jl 1つとっても使用例がやはり少ないです。 これから、みなさまとご一緒できればと思いますので、暖かく見守っていただけると幸いです。

続きについて

HMCとNUTSは続きでやっていきます。今回は時間切れで、オチなし…ということでした。完全に見積もりミスです。すみません。 次回は旅行記録ではなく、報告形式で!

続き

一旦のケリをつけたのでリンク。

Juliaでベイズ統計: Mambaのトライアル(アヒル本 第4章を参考に…)

Posted on 12/02/2017
このエントリーをはてなブックマークに追加

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

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

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

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

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

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

  • 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

PythonのGUI Toolkit比較(Tkinter, PyQt5, wxPython)

Posted on 8/26/2017
このエントリーをはてなブックマークに追加
PythonにはGUIツールキットが、いっぱいあって、定番を決められないうちに、いくつも使うようになってしまいました。Pythonのドキュメントにも、Tkinter, wxPython, PyQt, Kivy, さらにpygletだけにとどまらず、いくつものGUIツールキットが紹介されており、あたかも全てが公式のようです。サポートの面からも、いずれもコミュニティがあまり生き生きしていないようで、「このツールキットなら困ったときに、すぐ教えてもらえる!」みたいな選び方も難しいように感じます。

と愚痴りつつ、これらのうち、Tkinter, wxPython, PyQt, Kivy, Pygletは少しいじりますので、うち最初の3つで、簡単なGUIを作って比較します。