渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

MCMCの計算にStanを使ってみた（超基礎・導入編）

R 統計学 BUGS/Stan MCMC

肝心のMCMCの勉強はどこ行ったゴルァとか怒られるとアレなんですが、先にツールの使い方覚えてしまおうと思ってStanで簡単な練習をやってみました。ちなみに参考にした資料はこちら。


割とよく一緒に飲んでるid:EulerDijkstra氏のブログがとにかく役に立ちました。ありがとさんです！！！


あと、MCMCやるのはこれが初めてという人は最低限久保先生の緑本ぐらいは読んでおいて損はないと思います。ただしStanではなくWinBUGSを{R2WinBUGS}で回す系ですが。


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

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


メトロポリス法とかギブス・サンプラーについても一通りものすごーーーくわかりやすく書いてあるので、四の五の言う前にまずは読んでおきましょう。


MCMCとは何か（超絶端折ってますが）


MCMC自体の勉強はまだ途中なので、超絶端折って書いておきますが、要は

  • 統計モデリング（例えばGLMなど）で求めたいパラメータがあるけど解析的に求められない場合に、
  • ベイズ統計学の考え方に従って、
  • コンピューターによる乱数シミュレーションの手法で求める


手法です*1。正確な説明はいずれどこかで勉強が一段落したらまた書きますが、今回のところはひとまずこういう大ざっぱな理解でこのままいっちゃいましょう。


{rstan}をインストールする


インストール手順は普通にStan Projectのサイトに載ってます。英語ですが、難しいことは書いてないので読んだ通りの手順でインストールしていけば（途中コンパイルしたりするので時間がかかることがありますが）入るはずです。


ただし途中の{Rcpp}でハマる人がいる*2と思うので、僕が以前にハマってトラブルシュートした時の記事を先にどうぞ。


無事{Rcpp}が入れば、その後のインストール作業はほぼ確実にスムーズに走るはずです。途中延々と時間がかかるのはライブラリ周りのコンパイルなので、お茶でも飲みながら気長に待ちましょう*3。終わったら普通にrequireしておきます。

> require("rstan")
 要求されたパッケージ rstan をロード中です 
rstan (Version 1.3.0, packaged: 2013-04-12 21:12:02 UTC, GitRev: f57455593d14)


バージョンが1.3.0なのはご愛嬌（汗）。これでRからStanを動かす準備は整いました。


このブログで前にやってたロジスティック回帰の例で試してみる


何度かこのブログで使っているこちらのデータを使うことにしましょう。はっきり言って超簡単です。普通にインポートして"dat"という名前にでもしておきましょう。


ただし、Stanにかけたいのでcvカラムは0 or 1のnumeric型バイナリ形式に直しておく必要があります。どうやっても良いと思いますが、例えば

> dat[,8]<-as.numeric(dat[,8])-1


とすれば直せるはずです。これができたら、今度はStanに読み込ませるコードを書きましょう。文法についてはid:EulerDijkstra氏のslideshareにも詳しく載っているので、まずはそちらを。

data {
  int<lower=0> N;
  real x1[N];
  real x2[N];
  real x3[N];
  real x4[N];
  real x5[N];
  real x6[N];
  real x7[N];
  int<lower=0,upper=1> y[N];
}
parameters {
  real b0;
  real b1;
  real b2;
  real b3;
  real b4;
  real b5;
  real b6;
  real b7;
}
model {
  for (n in 1:N) 
      y[n] ~ bernoulli(inv_logit(b0 + b1*x1[n] + b2*x2[n] + b3*x3[n] + b4*x4[n] + b5*x5[n] + b6*x6[n] + b7*x7[n]));
}


ベタッと汚く書いてるのは、非エンジニアである僕の仕様です（泣）*4。

※追記その1


ベクトル演算できないのかなー、と思ってたらid:aaaazzzz036さんが教えて下さいました。

    data {
        int<lower=0> N;
        int<lower=0> M;
        matrix[N, M] X;
        int<lower=0, upper=1> y[N];
    }
    parameters {
        real beta0;
        vector[M] beta;
    }
    model {
        for (i in 1:N)
            // X[i] は row_vector, beta は vector だが, dot_product が吸収してくれる
            y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
    }


これでstanコードを糞コードにしなくても済みます！有難うございましたー。


注意すべき点はあまり多くないのですが、簡潔にまとめると

  1. data部 / parameters部 / model部の順に分けて書く*5
  2. data部で変数の型と範囲を宣言
  3. parameter部で推定したい変数の型を宣言
  4. model部で推定したいモデル式を具体的に書く
  5. 代入は"<-"、確率分布などによる分散を伴う関係は"~"
  6. 大抵の場合は「目的変数が○○分布に従う」というモデルの前提*6を表現するために、「目的変数 ~ 確率分布（説明変数のモデル式）」のフォーマットで書く


というところでしょうか。今回の例だとロジスティック回帰なので、目的変数cvがベルヌーイ分布（二項分布）に従うと仮定しています。なので、bernoulli()を使って上のようなフォーマットで書いてます。例えばこれを"conflict.stan"という名前で保存しておきます。ちなみにこれが正規線形回帰モデルなら、正規分布に従うのでnormal()を使うことになります。


次に、Stanに渡すデータをR側でリスト化しておきます。

> conflict.dat<-list(x1=dat$a1,x2=dat$a2,x3=dat$a3,x4=dat$a4,x5=dat$a5,x6=dat$a6,x7=dat$a7,y=dat$cv,N=3000)


アホらしいので、普通にattachしておけば良かったかも。なお、Stanは自動で配列の長さを取ってくれたりはしないのでちゃんとNを与えておくのを忘れずに。後はstan()に計算させるだけです。

> conflict.fit<-stan(file="conflict.stan",data=conflict.dat,iter=1000,chains=4)

TRANSLATING MODEL 'conflict' FROM Stan CODE TO C++ CODE NOW.
COMPILING THE C++ CODE FOR MODEL 'conflict' NOW.
cygwin warning:
  MS-DOS style path detected: C:/PROGRA~1/R/R-30~1.2/etc/x64/Makeconf
  Preferred POSIX equivalent is: /cygdrive/c/PROGRA~1/R/R-30~1.2/etc/x64/Makeconf
  CYGWIN environment variable option "nodosfilewarning" turns off this warning.
  Consult the user's guide for more details about POSIX paths:
    http://cygwin.com/cygwin-ug-net/using.html#using-pathnames
In file included from C:/Program Files/R/R-3.0.2/library/rstan/include/rstan/rstaninc.hpp:3:0,
                 from file11ac1de92a02.cpp:6:
C:/Program Files/R/R-3.0.2/library/rstan/include/rstan/stan_fit.hpp: In constructor 'rstan::stan_fit<Model, RNG>::stan_fit(SEXP, SEXP) [with Model = model11ac79001cdd_conflict_namespace::model11ac79001cdd_conflict, RNG = boost::random::additive_combine_engine<boost::random::linear_congruential_engine<unsigned int, 40014u, 0u, 2147483563u>, boost::random::linear_congruential_engine<unsigned int, 40692u, 0u, 2147483399u> >, SEXP = SEXPREC*]':

...

SAMPLING FOR MODEL 'conflict' NOW (CHAIN 1).
Iteration: 1000 / 1000 [100%]  (Sampling)

SAMPLING FOR MODEL 'conflict' NOW (CHAIN 2).
Iteration: 1000 / 1000 [100%]  (Sampling)

SAMPLING FOR MODEL 'conflict' NOW (CHAIN 3).
Iteration: 1000 / 1000 [100%]  (Sampling)

SAMPLING FOR MODEL 'conflict' NOW (CHAIN 4).
Iteration: 1000 / 1000 [100%]  (Sampling)


これでMCMCの計算が終わりました。結果が見たいだけなら、普通に出来上がったconflict.fitを呼べば良いだけです。

> conflict.fit
Inference for Stan model: conflict.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

       mean se_mean  sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
b0     -1.4     0.0 0.3   -1.9   -1.6   -1.4   -1.2   -0.9  1053    1
b1      1.1     0.0 0.2    0.7    1.0    1.1    1.2    1.4  1355    1
b2     -0.5     0.0 0.2   -0.9   -0.7   -0.6   -0.4   -0.2  1456    1
b3      0.1     0.0 0.2   -0.2    0.0    0.1    0.2    0.4  1425    1
b4     -3.0     0.0 0.2   -3.5   -3.2   -3.0   -2.9   -2.6  1040    1
b5      1.5     0.0 0.2    1.2    1.4    1.5    1.7    1.9  1384    1
b6      5.4     0.0 0.2    5.0    5.2    5.4    5.5    5.8   826    1
b7      0.1     0.0 0.2   -0.3    0.0    0.1    0.2    0.4  1425    1
lp__ -526.3     0.1 2.0 -530.8 -527.4 -526.0 -524.8 -523.4   727    1

Samples were drawn using NUTS2 at Wed Nov 06 12:38:51 2013.
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).


ちなみにglm()関数でロジスティック回帰した時の結果と見比べてみると、

> dat.glm<-glm(cv~.,dat,family=binomial)
> summary(dat.glm)

Call:
glm(formula = cv ~ ., family = binomial, data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.6404  -0.2242  -0.0358   0.2162   3.1418  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.37793    0.25979  -5.304 1.13e-07 ***
a1           1.05846    0.17344   6.103 1.04e-09 ***
a2          -0.54914    0.16752  -3.278  0.00105 ** 
a3           0.12035    0.16803   0.716  0.47386    
a4          -3.00110    0.21653 -13.860  < 2e-16 ***
a5           1.53098    0.17349   8.824  < 2e-16 ***
a6           5.33547    0.19191  27.802  < 2e-16 ***
a7           0.07811    0.16725   0.467  0.64048    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4158.9  on 2999  degrees of freedom
Residual deviance: 1044.4  on 2992  degrees of freedom
AIC: 1060.4

Number of Fisher Scoring iterations: 7


ほぼ同じ偏回帰係数が得られているのが分かるかと思います。ただしStanの結果から有意性*7を確認するには、2.5%点か97.5%点を自分の目で見て判断する必要があります。


{coda}パッケージで見やすくプロットする


ただ、WinBUGS + {R2WinBUGS}と違ってこのままではMCMCサンプリングされたデータの分布を見ることはできません。そこで{coda}パッケージを使ってこんなことをやってみましょう。

> require("coda")
 要求されたパッケージ coda をロード中です 

 次のパッケージを付け加えます: ‘coda’ 

 以下のオブジェクトはマスクされています (from ‘package:rstan’) : 

     traceplot 

> conflict.fit.coda<-mcmc.list(lapply(1:ncol(conflict.fit),function(x) mcmc(as.array(conflict.fit)[,x,])))


これで、{coda}パッケージのplot()関数でMCMCサンプルの分布を見ることができます。やってみるとこうなります。

> plot(conflict.fit.coda)


f:id:TJO:20131106124829p:plain

f:id:TJO:20131106124839p:plain

f:id:TJO:20131106125651p:plain


という感じで、分布を確認することができます。今回は割と綺麗な分布の形をしていますが、収束の悪いデータとかだと普通に多峰性になったりします。

※追記その2


id:aaaazzzz036さんの記事で、{coda}を使わずに可視化する方法が紹介されています。


{reshape}, {ggplot2}の組み合わせでいけるみたいです。こちらも多謝！

最後に


何でいの一番に普通の正規線形回帰モデルじゃなくてロジスティック回帰を選んだのか？みたいなことは突っ込まれても困るので念のため（笑）。


ちなみに、Stanのマニュアルにはその他の回帰や隠れマルコフモデル、果てはARMAなどの時系列モデルに至るまでStanで推定する方法が載っているので、おいおい試してみて結果をupしていこうかなと思います。気が向けば、ですが。。。

*1:端折り過ぎてほとんどデタラメ書いてる気がするのは多分その通りですごめんなさい

*2:Python(x,y)入れてると大体ハマる

*3:遅いマシンだと結構かかる。自宅のトロいノートPCでは15分ぐらいかかった記憶が。。。

*4:いつまでも言い訳するなというツッコミはご勘弁

*5:順番を守らないとエラーになる

*6:久保先生のサイトでも言及されてますが、○○分布に従うのは基本的には目的変数そのものですよー

*7:有意に0から離れているかどうか