stanの導入方法とベイズ分析の例 - バナナでもわかる話

前回まででベイズがどういうものかっていう話をしました。
bananarian.hatenablog.com

ベイズを簡単に試すにはStanというソフトを使うのが手っ取り早いです。
Stan モデリング言語: ユーザーガイド・リファレンスマニュアル

そこで、本記事ではStanをR上で利用するRstanの導入方法と、Stanを使ってどんなことが出来るのかという話をします。

ちなみにPythonユーザーの方であれば、PyStanというものがあるので、そちらで試してみると良いかと。

一応、Rはインストール済みという前提で導入方法の話をしていきます。

あ、Rをインストールしていない方はこちらでどうぞ。
R のインストール - RjpWiki

準備

Rでstanを使うには、R以外にC++コンパイラが必要なので、インストールします。

Windowsの方はRtoolsを(インストールは下のサイトから、自分のRのバージョンに合わせて行う)
https://cran.r-project.org/bin/windows/Rtools/

【注意】
(2018年9月時点での情報)
Pathを通す必要があるので、あまりコンピュータ関係の話について詳しくない人は、とりあえず表示されるチェックボックスは全てチェックして進んでください。

Macの方はXcodeを
Xcode - Apple Developer

※私はWindowsユーザーなので、Macの方はあまり詳しくありませんが、上のXcodeというもので大丈夫らしいです。

これで準備はオッケーです。Rでrstanパッケージを入れましょう。

install.packages("rstan")

Stanを試してみる

Stanコードの説明は次回に回すとして、今回はどういうことが出来るのかをざっとお見せします。

コードの説明は次回行うので、「へー、こんなことが出来るんだー」程度でごらんください。

前記事のベイズ流線形回帰を見てみます。
ベイズ統計学入門の入門 - バナナでもわかる話

$y=aX+\epsilon$

データとしてXとyが与えられているとします。この時のパラメータは $a$ ですと前回説明しました。

今回は更に、 $\epsilon$ の分散を $\sigma$ と置き、 $\sigma$ というパラメータも増やして考えてみます。

では、試しに回帰するデータを適当に作っておきます。

X=seq(1,10,0.1)
Y=5*X+rnorm(1)
data_list=list(N=length(Y),X=X,Y=Y)

Rstanで処理する場合は、データをリスト型にして一つにまとめておくのが良いです。
リストの中にdataブロックに書き込むものを入れておきます(dataブロックについては後述)

それではモデルをコードに直していきましょう。
今回は事前分布として一様分布を用いることとします。

library(rstan)
localModel_text="
data {
int N;
real X[N];
real Y[N];
}
parameters {
real a;      
real<lower=0> sigma;
}
model {
for(i in 1:N) {
Y[i] ~ normal(a*X[i], sqrt(sigma));
}
}
"

これでモデルの作成は終了です。コードがRっぽくないですね。これは要はStanのコードを文字列のように保存するという処理です。ここで、今回行う分析のモデルを書きます。

書いたモデルを実行(MCMCで計算)して、その出力をModel1.1という形で保持します。

Model1.1=stan(model_code=localModel_text,data=data_list,iter=2000,thin=1,warmup=1000,chains=3)

結果はこんな感じ

> Model1.1
Inference for Stan model: 3acd581e1fc643a3dd2e27505a5a9dd7.
3 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
a      5.08    0.00 0.00  5.07  5.07  5.08  5.08  5.08  2782    1
sigma  0.05    0.00 0.01  0.04  0.05  0.05  0.06  0.07  1757    1
lp__  87.06    0.03 1.04 84.41 86.64 87.35 87.79 88.08  1463    1

Samples were drawn using NUTS(diag_e) at Thu Sep 27 22:01:34 2018.
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).

何やらズラズラ書いていますが、今回注目すべきはココです。

       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
a      5.08    0.00 0.00  5.07  5.07  5.08  5.08  5.08  2782    1
sigma  0.05    0.00 0.01  0.04  0.05  0.05  0.06  0.07  1757    1

パラメータ $a$ と $\sigma$ をベイズ分析にかけた結果である、分布の要約です。

$a$ の平均・中央値が5.08
$\sigma$ の平均が0.05、中央値が0.06の分布が得られました。

実際のデータは $a=5$ , $\sigma=1$ の設定だったので、まずまず近い位置が平均・中央値として得られていますね。

stanでは、この結果を計算するにあたってMCMCを利用しています。
MCMCについて詳しくはこちら
マルコフ連鎖の不変確率分布とMCMCとの関係 - バナナでもわかる話
※続きものですが、一応この記事の真ん中らへんにMCMCの説明があります。

MCMCがしっかり機能したかどうかをプロットによって確認する必要があるので、その辺も確認しておきます。
※詳しくはまた次回

traceplot(Model1.1)

f:id:bananarian:20180927221433p:plain

おおざっぱですが、こういう感じで、一定の範囲にまとまって分布している場合は大体大丈夫です。

更にggmcmcパッケージというパッケージを使って、出力結果とMCMCの結果に関する一連の情報をpdfにまとめることもできます。

save.image(file='Model11.RData')
install.packages("ggmcmc") 
library(ggmcmc) 
 
load('Model11.RData') 
 
write.table(data.frame(summary(Model1.1)$summary),file='fit-summary.txt',sep='/t',quote=FALSE, col.names=NA) 
 
ggmcmc(ggs(Model1.1, inc_warmup=TRUE, stan_include_auxiliar=TRUE),   file='fit-traceplot.pdf', plot='traceplot') 
 
ggmcmc(ggs(Model1.1), file='fit-ggmcmc.pdf')