読者です 読者をやめる 読者になる 読者になる

Knowledge As Practice

JAIST(東京)でサービス経営の研究をしている社会人大学院生の研究・勉強メモ(統計分析多め)。

非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()

f:id:hikaru1122:20150815181025j:plain

 
では、いよいよ Stan を使って単回帰分析を行います。まず分析するためのデータを作ります。必要なのは、サンプルサイズ、最高気温、客数です。リスト形式でまとめて、dat という名前をつけておきます。Stan で使うデータはリスト形式にしておくのがお約束のようです。

dat <- list(N = 20, kion = icecream$kion, kyaku = icecream$kyaku)

 
次に、Stan のコードを作ります。切片を beta0、最高気温(kion)の回帰係数を beta1 とします。要は「 { \displaystyle
客数 = \beta_{0} + \beta_{1}*気温
} 」という式を作るわけです。

 
作ったコードは 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で何ができるか」を眺めてみる - 銀座で働くデータサイエンティストのブログ

4.4 回帰直線を求めてみよう

クリエイティブ・コモンズ・ライセンス
この 作品 は クリエイティブ・コモンズ 表示 - 継承 4.0 国際 ライセンスの下に提供されています。