Stan Advent Calendar 2016 の 4日目です。非マニア枠です。
本日のエントリーは Stan を普通に使える人、情報科学・データ科学な人を対象にしていません。少し R に慣れてきて、ベイズ統計モデリングにも興味を持ち始めた人向け。ずばりターゲットは次のようなレベルの人です。
いざ R を使って分析するときには使い方を忘れていて検索しながら使う。
dplyr の使い方に慣れておらず、めんどくさくなって結果エクセルを使ってしまう。
つまり私。そんな中途半端なレベルな人（私）向けに『Stan と R でベイズ統計モデリング』（以下，アヒル本とします*1）を読んでいるときにつまずくであろうところを書いておきます。アヒル本を読む・理解するには、適度に R と ggplot2 を使えるスキルが必要だからです。
さて，R を使った統計分析の本をいくつか読んで、簡単な分析をできるようになったという人が Stan を使ったベイズ統計モデリングをするときにつまずくのは次のようなところだと思います。
- 結果が出るまで待ち時間が長い。
- リスト形式でデータを渡すのを戸惑う。
- 結果が入ったオブジェクトを表示させたらズラーッと出てきてしまい、見たい結果を見るために必死に上へスクロールする。
- タイプミスする。
- 結果からほしいデータを取り出せない。
1. 待ち時間が長い
私の場合，簡単なモデルでも30秒くらいかかります。
lm()と違って，パッと結果が帰ってこないので，学習のリズムが狂いますよね。これは仕方ない。オススメは BGM をかけることです。ガンダム UC で流れていた MOBILE SUIT です。無事ちゃんと結果が出るか，ドキドキ感が増します。
2．リスト形式にとまどう
真の R の初級者はあまりリスト形式に触れることはありません。だいたいデータフレーム形式で間に合うので。lm(y ~ x, data=d)な感じで。しかし，Stan を使う場合，
> data = list(N=row(d), A=d$A, Score=d$Score/200, M=d$M)
というように分析に使うデータを別に作らなければなりません。「リスト形式ってなんだよ、めんどくさい」となるんですよ*2。
リスト形式とは、Todo リストみたいなものだと思ってください*3。Todo リストには、いろんな種類、内容を書けます。R のリスト形式も同じ。いろいろなデータをヒトマトメにできるのです。上の list だと N には1つの数値しか入っていませんが、A や Score にはたくさんの数値が入っています。そういうのをまとめて1つのものにするのがリスト。
Stan を使うには専用のデータを作ってあげる必要があって、それはリスト形式なんだと納得しましょう。慣れるしかありません。
3．結果の表示が多過ぎ！
（RStudioの場合）「サンプリングもうまくいった。よし，結果見てみよう」
> fit
どばっと結果が出ます。そしてパラメータを確認したいので上へ必死にスクロール。
こういうのを10回くらいやっていて「ひょっとしたらいい方法あるんじゃないか」と思いました。ありました。たとえば、b1・b2 というパラメータを見たい場合は次のようにします*4。
> print(fit, par="b") Inference for Stan model: model5-6b. 4 chains, each with iter=2000; warmup=1000; thin=1; post-warmup draws per chain=1000, total post-warmup draws=4000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat b[1] 3.57 0 0.10 3.39 3.51 3.57 3.63 3.76 976 1.01 b[2] 0.26 0 0.04 0.19 0.24 0.26 0.29 0.34 1820 1.00 b[3] 0.30 0 0.15 0.01 0.20 0.30 0.40 0.59 967 1.01 Samples were drawn using NUTS(diag_e) at Sat Dec 3 10:54:04 2016. 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).
事後対数尤度も見たいというなら、
> print(fit, par=c("b","lp__"))
便利ですね。Rhat を小数点とか細かく見てみたいなという場合は例えばprint(fit, par="b", digits=3)とすればOKです。詳しい人は「当たり前」と思ってしまうレベルですが、真の初級者にはこういう tips を知っておくと，途中でめんどくさくなってやめるリスクが減ります。
4．タイプミス
私の場合、1つの Stan ファイルを実行するのに4～5回くらいエラーが出ます。本当にささいなミスで。英単語のつづりが間違っていたとか、閉じカッコが足りなかったとか、引用符がなかったとか*5。これの解決法はわかりません。誰か教えてほしいです。print を pritn とやっちゃうくらいなので。
5．ほしいデータを取り出せない。
アヒル本にあるようにrstan::extract(fit)で結果を取り出せます。取り出せますが，そのデータをいじるのには慣れが必要です。「fit の中身ってどうなっているんだろうと」と思ってstr()を使ってもよくわからない。なんか表示結果に圧倒されてしまう･･･。
> ms = rstan::extract(fit) > str(ms) List of 4 $ b : num [1:4000, 1:3] 3.48 3.57 3.56 3.52 3.62 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ lambda: num [1:4000, 1:50] 3.64 3.69 3.67 3.66 3.71 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ m_pred: num [1:4000, 1:50] 34 40 54 32 43 32 32 40 36 40 ... ..- attr(*, "dimnames")=List of 2 .. ..$ iterations: NULL .. ..$ : NULL $ lp__ : num [1:4000(1d)] 6897 6896 6898 6898 6897 ... ..- attr(*, "dimnames")=List of 1 .. ..$ iterations: NULL
数字とアルファベットがたくさんあって，イラッときますね。でも落ち着いて。普通にfit = stan(～～～)としたときは合計4000個のMCMCサンプルが入っています。$b のところを見ましょう。推定したパラメータ b はb1, b2, b3の3つある（「1：3」のところ）。4000行3列にまとまっている状態です。なので，b1 のヒストグラムはhist(ms$b1[,1])で書けます。最大事後(MAP)推定値も max(ms$b1[,1]) で出ます。
ggplot を使うときはちょっとやっかいです。そもそもms$b[,1]はデータフレーム形式ではないから、そのまま ggplot には使えません。いったんデータフレーム形式に変換します。data.frame()のカッコの中が大事。データフレームを作ったときの列名をつけておきます。
> b1 = data.frame(b1 = ms$b[,1]) > str(b1) 'data.frame': 4000 obs. of 1 variable: $ b1: num 3.48 3.57 3.56 3.52 3.62 ... > head(b1) b1 1 3.484667 2 3.572416 3 3.559215 4 3.521891 5 3.615658 6 3.657821
b1という列名がついた4000行1列のデータフレームになりました。ggplot でヒストグラムを作るなら
> ggplot(data=b1) + geom_histogram(aes(x=b1)) + theme_bw()
というように、いつもどおりの方法でできます。本エントリーの内容がわかっていれば，アヒル本のコードも読み解けるはず。コツコツ取り組んでいくと、Stan への理解が深まるだけでなく、R の操作や作図能力もぐっと上がると思います。
最後に
周りに統計に詳しくてすぐに質問できる人がいない社会科学系の社会人学生の方はいっしょにがんばっていきましょう！ そして、積極的に読書会や勉強会に参加しましょう！ 意外と年齢層高いから大丈夫。
Osaka.Stan は12月23日のようです（正式な募集はまだ）。
Hommachi.Stats は12月10日です。
