ども。久しぶりの投稿です。
今更ながらベイズでのパラメータ推定をしっかり学んでみようと、
アヒル本を読み進めて、一旦いいところまで来たので自分用思考整理メモ。
まぁ、数多あるStan解説blog/メモblogの一つだと思っていただければ・・・
なお、コードとかの話よりは概念に近いです。
余裕があったら別エントリでコード載せるかもしれません。
Stanとは
Stan とは MCMCを用いたパラメータ推定の為のソフトウェアー。
HPはこちらStan、アヒル本にあるようにRで使う人はRStanを、PythonユーザーはPyStanをどうぞ。
Stan 自体はC++で書かれているコンパイル型なのですが、色々な言語に対応しているようで、MatlabやJuliaにも対応しています。
ですが、SinyStanやrstantoolなど、拡張ツールがRベースなので、基本的にはRと親和性が高いのだと思います。
何ができるの？
MCMCのchainを複数生成し、chainからサンプルされた分布を用いて、回帰式などのパラメータの推定ができます。ベイズによるパラメータ推定(と言っていいのか？)です。
ベイズによるパラメータ推定とは
パラメータの推定法(最小二乗法、最尤推定、ベイズ推定)はそれぞれ噛み砕くとこういった思想だと解釈しています。
最小二乗法「回帰式から得た予測とモデルの差を最小にするものが良いパラメータである」
最尤推定の「得られたサンプルの分布(ヒストグラムなどで確認できる)に一番当てはまりのいいものがよいパラメータである」
ベイズ推定「平均などのパラメータは確率分布しており、その分布(事後分布)は事前分布と尤度の積に比例しているよ」
R.A.フィッシャー以来のメインストリームである「真の平均/母集団の平均は一定である」と言う思想に対して、
ベイズ推定は「真の平均も一意ではないよ。確率的に定まっているよ」と言う考え方を採用して作られている為、割と思想の根本で対立していますが、まぁ実務面では現実に起きている現象を適切にモデル化できればなんだっていいわけで、その辺は学会の専門家に任せたいと思います。
MCMC 法とchainが収束するは？
Markov Chain Monte Carlo 法が正式名称です。MCを使ったMC法の事です。
丁寧に話すと、
MC(マルコフ連鎖)1時点前の確率変数にのみ依存し現在の値が決まる確率過程です
MC(モンテカルロ)乱数を使ったシミュレーション/数値推定方法の全体です
ん？まだ難しい？そんな人は
を見てください。ともあれ、MCMCとはマルコフチェインを使ったモンテカルロ法の事です。
ぶっちゃけた話、直感的には難しいことしてなくて、無茶苦茶たくさん乱数ベースのシミュレーションして、そこからパラメータの分布(事後分布)を推定しようぜー。って話です。
パラメータ推定には複数のchainを生成し、収束した部分のサンプルを抽出して利用します。
なので、以下の注意点があります。
・収束していない段階のサンプル(warm up)を使ってはならない
・chainが収束しない事もあるので、その場合は推定してはならない
上は一度図を見れば理解できるともいます。 アヒル本や緑本などで確認できるので本屋で探して見てください。
下は少し説明します。
MCと言うのは、状態とその遷移を確率的に表した確率遷移行列で表現されていて、Wikipediaのここを覗いてもらいたいのですが、既約なマルコフチェーンは定常過程に収束します。
既約とは、「全ての頂点が互いに繋がっている」つまりステップをたくさん踏めばマルコフ連鎖内にある全ての状態から別の状態にたどり着ける形であり、これは「確率遷移行列が収束する」状態です。
ここで注意しておきたいのが、既約であっても収束しない状態があると言う事です。この場合は「確率線維行列に周期性がある」状態です。
他方、既約でない場合はchainが収束しないので、パラメータの確率分布が一意に定まらないことを示しています。分布が一意に定まらないのに推定しちゃダメってことです。
RStanでは、パラメータ推定時にうまく収束したかどうかについては Rhat が1の付近の値になっているかを見ることで判定できます。
Stanで扱えるモデルは？
結論、どこまで扱えるかわからないですが、
アヒル本の8章までを読んだ結果、とてもフレキシブルで拡張性の高いソフトなので、
一般的な多変量線形/非線形回帰で階層型のモデルなら大抵扱えるんじゃないかなと思ってます。
ただ、モデルが複雑になると、その分計算時間は増えますし、収束もしなくなるので複雑にしすぎるのも悩みどころのようです。
Stanの書式基礎(1)
RStanはRファイルとは別に同じ階層に”.stan”と言う拡張子で作られたテキストファイルを用意します。
例えばここにもっとわかりやすい資料があるのですがStanコードには6つのBlockと呼ばれる箱があり、概ね4つのblockを使うようです。
Data{}ブロック：データを宣言します。C++で書かれているため、INTやREALなど型宣言もここでします
Transformed data{}ブロック:渡したデータを加工する場合はここで宣言します
Parameters{}ブロック：推定したいパラメータの宣言です
Transformed paramaters{}ブロック：パラメータの加工はこちらで宣言します
Model{}ブロック：モデルを記述します
Generated quantities{}ブロック：追加で確認したい統計値を宣言。95%区間を求めるなどができる
必須なのは model ブロックだけのようですが、大抵はtransformed 系ブロック以外は使っているようです。
ちなみに、自作の関数で拡張できる functions{}ブロックもあるみたいです。
結構長く書いてしまったので、ここで終えます。
GW中にもう一本実装踏まえて書いて見たいんですが・・・できるかな？