『RとStanではじめるベイズ統計モデリングによるデータ分析入門』
統計モデリング系のお勉強は、昔『みどり本』*1を読んだ後『アヒル本』*2を半分くらい読んでほったらかしにしていました。実務で必要に迫られたこともなかったので、いつかやろうと先延ばしにしていたのですが、今回良い機会なので、この本から再入門することにしました。
勉強した本の紹介方法は、現状手探りでやっていこうかなと思っていて、今回は各部毎の内容を一部抜粋でまとめつつ、最後に全体の感想を書くようなスタイルでやってみようかなと思います。*3
尚、私の読解力、知識不足により誤読している箇所があるかもしれないので、もし見つけられた方は優しく諭して頂けますと幸いです。
内容紹介
1部 【理論編】ベイズ統計モデリングの基本
1部は理論編ということで、確率や統計の基礎事項、モデリングとはなんぞや、事後分布に従う乱数を作成するためのMCMC、などをコンパクトに紹介した部です。
”5章からはちゃんと読んでね”とある通り、ここから統計モデル構築のための準備が始まります。大事なのがベイズ更新の考え方ですかね。
ベイズ更新
が事前確率、つまりデータが得られる前に想定していた確率ですね。
が事後確率、つまりデータが得られた後に想定する確率となります。これらを結びつけるのが、
の部分です。分子が尤度、分母が周辺尤度となります。つまり考えうる様々な仮定の元でデータが得られる確率のうち、今回想定している仮定
の元でデータが得られる確率を考えることで、事前分布を事後分布に更新しているわけです。
ベイズ推論とMCMC
パラメータの事後分布が推定できれば、万事解決のように見えるんですが、まだ問題があります。それが、事後分布が複雑過ぎて扱いにくい問題です。例えばパラメータの期待値を求めるのに
この計算をしないといけないですが、先程の事後分布が複雑過ぎ問題があるため容易に計算ができなかったりします。そこで、事後分布に従う乱数をいくつか生成して、その平均を取ることで代替しようということを考えます。そのために使うのがMCMCというわけですね。
乱数発生アルゴリズムについては、本書ではMH法とHMC法が紹介されています。MH法のざっくりした考え方としては、以下のような感じです。
手順の2.~3.を適当な回数繰り返して乱数を得ます。カーネルというのは、事後分布の分子であるのことです。分母の周辺尤度は無くて良いのかと思われるかもしれませんが、これは上記の手続きで”比”をとっているため必要がないのですね。この点は結構大事で、そもそも周辺尤度の計算自体が難儀なため、使わないで済むととても有り難いわけです。
2部 【基礎編】RとStanによるデータ分析
2部は前半は基本的なR、可視化周り、stanの使い方の解説が中心です。後半では実際にstanを使って簡単な事例に取り組みます。そうそう、こんな感じだったなあとか言いながら写経しました。対数密度加算文形式の表記もちゃんと紹介されています。
target += normal_lpdf(データ|パラメータ)
3部 【実践編】一般化線形モデル
いよいよ実践編です。まずは一般化線形モデルについて概要がまとめられており、その後実際にStanによる分析例をいくつか追っていくような形式です。こういうのは慣れもあるので、適度な分量を写経できる感じになってて有り難いです。
最後の交互作用の章にある売り上げと製品数、店員数の関係性の分析は、交互作用を考えることで関係性が浮き彫りになる面白い例です。また、結果が因果関係ではないことを指摘されているのもとても大事なポイントだと思います。*4
線形予測子・リンク関数
以下の一般化線形モデルで理解するのがわかりやすいです。
が線形予測子で、
がリンク関数ですね。本書の例ではリンク関数は恒等関数になってますが、ここを色々弄れるわけですね。
brms
brmsパッケージを使うとstanファイルを書かなくてもお手軽にベイズ統計モデリングが出来ます。ただ、状態空間モデルなど現状未対応なモデルもあるため、stanファイルで記述するベーシックな書き方もちゃんとマスターしといた方が良いようです。
以下に実行イメージを掲載しておきます。*5
library(rstan) library(brms) # このパッケージ rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) # データ読み込み data <- read_csv("hoge.csv") # モデル作成 result_brms <- brm( formula = y ~ X, family = gaussian(link = "identity"), # 確率分布 data = data, seed = 1 ) # MCMCサンプルの取得 # as.mcmc(result_brms, combine_chaines = TRUE) # 事後分布の確認 plot(result_brms)
4部 【応用編】一般化線形混合モデル
ここからは階層ベイズモデルを扱います。3部で扱ったポアソン回帰で 過分散が起きていた場合の対処方法として、一般化線形混合モデルを用いて対処する方法が紹介されています。
固定効果・ランダム効果
3部で推定してきたパラメータ（切片、”天気が釣獲尾数にもたらす効果”、”気温が釣獲尾数にもたらす効果”といった部分）が固定効果です。 一方線形予測子に加えられた調査ごとに変化するようなランダムな影響をランダム効果と言います。
固定効果、ランダム効果両方を考慮するモデルを混合モデルと言います。 つまり一般化線形混合モデルというのは、一般化線形モデルにランダム効果を加味したモデルということですね。
5部 【応用編】状態空間モデル
最後は状態空間モデルです。ローカルレベルからトレンド、周期性、自己回帰モデル、動的一般化線形モデルと様々なモデルを試せるようになっています。
状態モデル・観測モデル
状態空間モデルは”状態モデル”と”観測モデル”から成ります。まず状態モデルは以下のように表せます。
つまり、t時点の状態がt-1時点の状態に影響されるということですね。この状態に誤差が加わったものが観測値と考えるわけです。つまり観測モデルは以下のように表せます。
ローカルレベルモデル
ローカルレベルモデルは以下のようなものです。
を過程誤差、
を観測誤差と言います。
は正規ホワイトノイズに従う過程誤差が積み重なっているような系列です。つまりランダムウォークですね。因みにホワイトノイズというのは、期待値が0、分散が一定で、同時刻以外の自己相関が0となっているものです。さらに”正規”なのでこれはi.i.d系列になっているのでした。
基本構造時系列モデル
観測値が「トレンド成分 + 周期成分 + ホワイトノイズ」で表現されるモデルを基本構造時系列モデルと言います。例えば、1次のトレンドに周期kの周期成分が入ったモデルは以下のようになります。
がトレンド、
が周期成分ですね。周期成分について補足をすると、周期kのk個の周期成分の合計が満たす性質、
に注意して、
となります。
全体の感想
今まで統計モデリングの本というと、『みどり本』と『アヒル本』が2強という感じでした。私は『アヒル本』はまだ全部読んではいないんですが、入門という観点だと本書がベストな気がしています。
というのも、とにかく例が豊富で、コード解説も初心者にもわかりやすく工夫されているので、本書を実際に手を動かしながら読み進めることで、stanの基本的な使い方や統計モデリングの勘所が自然と身に付く構成になっていると思います。ただ、当然入門書なので、細かいところ*6は追加で適宜勉強していく必要があります。
とはいえ本書だけでもまあまあお腹一杯感があるんですが（笑）これで満足せず、引き続きベイズ統計モデリングを学んでいければなと思います。
