INPUTしたらOUTPUT!

忘れっぽいんでメモっとく

第9回 「データ解析のための統計モデリング入門」 読書会に参加してきた

統計

ニコ生見直して復習していたので遅くなったけど以下メモ

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)


第9章 「マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデル」

前半


  • P.202でloadされているd.RData、現在はNot Foundになってる。

    • 北大のサーバーメンテの影響かと思ってたけど違うみたい
    • 著者のHPには下記のとおり全ファイルまとめたものもアップされているのでそちらから入手することは可能

      以下の全ファイル (+ 作図関連の余計なファイル) をまとめてダウンロードできるようにしました: kubobook_2012.zip (26 MB, 2013-12-27)

  • ベイズモデルの事後分布は尤度x事前分布に比例する(P.196)

    • {p \left( \beta_1 ,\beta_2 \mid Y \right) \propto p \left( Y \mid \beta_1 ,\beta_2 \right) p \left( \beta_1 \right) p \left( \beta_2 \right)}
      • 事後分布 : {p \left( \beta_1 ,\beta_2 \mid Y \right)}
      • 尤度 : {p \left( Y \mid \beta_1 ,\beta_2 \right)}
      • 切片{\beta_1}の事前分布 : {p \left( \beta_1 \right)}
      • 傾き{\beta_2}の事前分布 : {p \left( \beta_2 \right)}
    • 切片{\beta_1}, 傾き{\beta_2}を適切に指定すれば事後分布が得られる

  • 切片{\beta_1}、傾き{\beta_2}の確率分布は分からない

    • {-\infty}から{\infty}までの範囲で好きな値をとる分布（無情報事前分布）を仮定する

  • Stanの書き方が分からなくて困ってました

    • me too!
    • glmer2stanパッケージ使うとlme4形式でコーディングできる
      • デフォルトだとchain=1なので注意する
      • @stanmodelでstanのモデルコードを表示してくれる


後半


  • Windows 8ではWinBugs動かなかった
    • 代わりにStanにした
  • Stanの場合{\beta_1}、{\beta_2}が定数なのか実数なのか指定しないと怒られる
  • codaで可視化
    • 1段目 ：{\beta_1}
    • 2段目 ：{\beta_2}
    • 3段目 ：対数尤度
  • サンプルの結果を抽出するにはextract()
  • 信頼区間の考え方
    • {\beta_1} の事後確率が95%で1.81〜2.13になる
      • 正しい
    • {\beta_1} が1.81〜2.13にある確率は95%
      • 誤り
    • {\beta_2}の信頼区間が-0.05〜0.22となっている
      • 0を挟んでいるので十分に0から離れているとは言えないのではないか？
  • パラメータが複数ある場合は基本的に1つずつサンプリングする
    • サンプリングの順序としてはどれから行っても問題ない
      • {\beta_2}を定数としてまず{\beta_1}をサンプリングする
  • メトロポリス法とギブスサンプリングの違い
    • メトロポリス法 ・・・ 新しい値の候補をあげ、それに変化するかどうかを決める
      • 尤度が上がれば値を更新、上がらなくてもランダムで更新
    • ギブスサンプリング ・・・ 新しい値の確率分布を作り、その確率分布のランダムサンプルを新しい値とする
  • ギブスサンプリングのメリット
    • 更新前後の値の相関が比較的小さい
      • 相関が高いとサンプルをとっても捨てられてしまう
    • 共役事前分布が使えるので計算が高速
      • 尤度 : 平均{\lambda}のポアソン分布
      • 事前分布 : パラメータ{\alpha}、{\beta}のガンマ分布
        • とすると
      • 事後分布 : パラメータ{\alpha ' = \alpha+N}、{\beta ' = \beta + \sum y_i}のガンマ分布
        • となりガンマ分布の乱数の発生方法は既知なので計算しなくてもサンプルがとれる(＝高速）
  • 共役である = 事前分布と事後分布が同じ種類の確率分布となるような場合
  • 共役事前分布って使った方が良い？


LT

マルコフ連鎖モンテカルロ法入門-１



可視化で理解するマルコフ連鎖モンテカルロ法


  • メトロポリス法だと尤度が上がらない場合、値を更新しないのが視覚的によく分かった！
  • ギブスサンプリングで複数パラメータを推定する場合、一方を定数としてランダムサンプルするのが視覚的によく分かった！
  • 今回は受付として参加させて頂いた
    • 運営側の苦労が少し分かりました。。。
    • ピザの分量は最適化されたよう
  • MacではWinBugs動かないのでRStan使ってMCMCの計算をやってみた(1) - データ解析勉強日記を参考にStanで動かしてみた。
    • GLMで推定した場合と結果はほぼ同じ
    • Rhat = 1なので収束しているとみなせる
      • Rhat > 1.1ではサンプル列間のばらつきが大きいので定常分布・事後分布の推定はできないと判断する
> #まずはGLM
> load("d.RData")
> y <- d$y
> x <- d$x - mean(d$x)
> model <- glm(y~x, family=poisson)
> summary(model)

Call:
glm(formula = y ~ x, family = poisson)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.52168  -0.53195   0.06417   0.40797   1.57939  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.98276    0.08318  23.836   <2e-16 ***
x            0.08334    0.06838   1.219    0.223    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 15.660  on 19  degrees of freedom
Residual deviance: 14.171  on 18  degrees of freedom
AIC: 94.036

Number of Fisher Scoring iterations: 4

> # Stanで実行
(中略）
> d.fit
Inference for Stan model: kubo9.
3 chains, each with iter=1600; warmup=100; thin=3; 
post-warmup draws per chain=500, total post-warmup draws=1500.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta1   1.97    0.00 0.08   1.81   1.92   1.97   2.03   2.13  1500    1
beta2   0.08    0.00 0.07  -0.05   0.04   0.08   0.13   0.22  1500    1
lp__  143.97    0.03 0.98 141.23 143.60 144.27 144.65 144.94  1143    1

Samples were drawn using NUTS(diag_e) at Fri Sep 09 08:50:27 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

> # CRANにはないのでGithubからインストールして読み込み
> library(devtools)
> install_github("rmcelreath/glmer2stan")
> library(glmer2stan)

> # Stan実行
> res <- glmer2stan(y~x, family="poisson", iter=1600, warmup=100, thin=3, chain=3)
> (res)
Inference for Stan model: y ~ x [poisson].
3 chains, each with iter=1600; warmup=100; thin=3; 
post-warmup draws per chain=500, total post-warmup draws=1500.

            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
Intercept   1.98    0.00 0.08   1.81   1.92   1.98   2.03   2.13  1439    1
beta_x      0.08    0.00 0.07  -0.04   0.04   0.08   0.13   0.21  1416    1
dev        92.01    0.06 2.06  90.09  90.62  91.32  92.69  97.71  1273    1
lp__      143.99    0.03 1.03 141.14 143.65 144.33 144.68 144.95  1273    1

Samples were drawn using NUTS(diag_e) at Sat Sep 13 12:50:34 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

# stanのモデルコードを表示
> res@stanmodel
S4 class stanmodel 'y ~ x [poisson]' coded as follows:
data{
    int N;
    int y[N];
    real x[N];
}

parameters{
    real Intercept;
    real beta_x;
}

model{
    real vary[N];
    real glm[N];
    // Priors
    Intercept ~ normal( 0 , 100 );
    beta_x ~ normal( 0 , 100 );
    // Fixed effects
    for ( i in 1:N ) {
        glm[i] <- Intercept
                + beta_x * x[i];
        glm[i] <- exp( glm[i] );
    }
    y ~ poisson( glm );
}

generated quantities{
    real dev;
    real vary[N];
    real glm[N];
    dev <- 0;
    for ( i in 1:N ) {
        glm[i] <- Intercept
                + beta_x * x[i];
        dev <- dev + (-2) * poisson_log( y[i] , exp(glm[i]) );
    }
}