非ITエンジニアのための Stan と Rstan でベイズな単回帰分析
『基礎からのベイズ統計学』には単回帰分析の例がなかったので、Stan による単回帰分析の練習をしました。今回のエントリーはR および統計の初心者は抜けだした、『基礎からベイズ統計学』を読んでだいたいわかった社会科学系の方が対象です。ITエンジニアやデータサイエンティストの方々のブログ、難しいんですよね…。なので、シンプルな分析をやってみました。
データは「アイスクリーム統計学」からお借りしました。ハンバーガー統計学とともに、わかりやすくていい本です。データはこちらのページを参照してください。
http://kogolab.chillout.jp/elearn/icecream/chap4/sec1.html
データを Excel にコピペして、タブ区切りのCSV形式で出力。それを Rstudio で読み込みました。ええ、このへんが素人感丸出しです。
> icecream id kion kyaku 1 33 382 2 33 324 3 34 338 4 34 317 5 35 341 6 35 360 7 34 339 8 32 329 9 28 218 10 35 402 11 33 342 12 28 205 13 32 368 14 25 196 15 28 304 16 30 294 17 29 275 18 32 336 19 34 384 20 35 368
最高気温（kion）を説明変数、お客さん（kyaku）を目的変数として、Stan で単回帰分析を行います。サンプルサイズは20です。「お客さんの数はマイナスにはならないじゃん！」とか無粋なことは言わないでください…。あくまで Stan に慣れるためなので。その前に、まず見える化。
icecream %>% ggvis(~kion, ~kyaku) %>% layer_points()
では、いよいよ Stan を使って単回帰分析を行います。まず分析するためのデータを作ります。必要なのは、サンプルサイズ、最高気温、客数です。リスト形式でまとめて、dat という名前をつけておきます。Stan で使うデータはリスト形式にしておくのがお約束のようです。
dat <- list(N = 20, kion = icecream$kion, kyaku = icecream$kyaku)
次に、Stan のコードを作ります。切片を beta0、最高気温（kion）の回帰係数を beta1 とします。要は「
」という式を作るわけです。
作ったコードは icecream.stan という名前を付けて、R の作業ディレクトリ（ getwd() で表示される場所）に保存します。お客さんの数（kyaku）が正規分布に従っているとして、その母数（平均と標準偏差）を推定するためのコードです。平均を mu、標準偏差を sigma とします。
data { int<lower=0> N; #サンプルサイズ。0より大きい。 real kion[N]; #最高気温のデータは20。 real<lower=0> kyaku[N]; #データの客数は20人。0より大きい。 } parameters { real beta0; #切片。 real beta1; #最高気温の回帰係数。 real<lower=0> sigma; #標準偏差。 } transformed parameters { #客数の平均は切片、気温の回帰係数、気温を使って定義していく。 real mu[N]; for (i in 1:N) { mu[i] <- beta0 + beta1*kion[i]; #客数の平均を定義。 } } model { kyaku ~ normal(mu, sigma); #客数は平均 mu、標準偏差 sigma の正規分布に従う。 }
では分析実行。分析結果を fit に保存します。サンプリング回数などは『基礎からベイズ統計学』の第7章章末問題のコードにならいました。
library(rstan) scr<-"icecream.stan" par = c("beta0", "beta1") war<-5000 ite<-100000 cha<-1 fit<-stan(file = scr, model_name = scr, data = dat, pars = par, chains = cha, warmup = war, iter=ite)
じっと待つこと10秒くらい。無事、分析が終わりました。結果を出力します。
print(fit,pars=par) mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat beta0 -230.75 0.57 80.23 -389.66 -282.82 -230.41 -178.89 -72.06 19789 1 beta1 17.27 0.02 2.50 12.32 15.66 17.26 18.90 22.23 19771 1
普通の回帰分析 lm(kyaku ~ kion, data = icecream) の結果も載せます。
Estimate Std. Error t value Pr(>|t|) (Intercept) -229.98 73.79 -3.117 0.00596 ** kion 17.25 2.30 7.499 6.08e-07 ***
ほぼ、同じ結果になりました。だから何だと言われても困ります。とりあえず練習してみたかったので…。今回の分析から「最高気温が1度上がるごとに増えるお客さんの数は平均は17人ほどです。また、最高気温1度が上がるごとに増えるお客さんの数は、95％の確率で12～22人の範囲に収まります」と言えます。2文目が言えることがベイズの利点ですね。
今回のエントリーで参考にしたのは、次のページです。ありがとうございました。
RPubs - Multiple linear regression with Stan
Stanで統計モデリングを学ぶ(3): ざっと「Stanで何ができるか」を眺めてみる - 銀座で働くデータサイエンティストのブログ
