Knowledge As Practice

JAIST(東京)で Transformative Service Research に取り組んでる社会人大学院生の研究・勉強メモ

非ITエンジニアのための Stan と Rstan でベイズな重回帰分析

前回に続いて、Stan と Rstan で重回帰分析を行います。単回帰分析ができれば、重回帰分析はすんなりできます。っていうか、できました。変数が増えただけ、だからです。ほんのちょっぴり、Stan コードを洗練させました(IT系の方には「どこが?」と言われそうだけど)。

 
Stan で単回帰分析をしたことがない方は、前回のエントリーを見てください。先に見たほうがこのエントリーを理解しやすいと思います。 hikaru1122.hatenadiary.jp

 
データは再び「アイスクリーム統計学」からお借りしました。ほんとにいい本ですよね。 http://kogolab.chillout.jp/elearn/icecream/chap5/sec1.html

 
まず分析するデータを上記「アイスクリーム統計学」のページから、Excel にコピペ。タブ区切りのテキスト形式 icecream2.txt で保存して、Rstudio に読み込みます。何も考えなくていいので、楽チンです。でも、きっとRだけで行える方法があるんですよね…。誰が教えていただけないでしょうか?

 
なお、データ番号、最高気温、最低気温、客数はそれぞれ id、saiko_kion、saitei_kion、kyaku という列名にしています。

> icecream2
   id saiko saitei kyaku
 1    33     22   382
 2    33     26   324
 3    34     27   338
 4    34     28   317
 5    35     28   341
 6    35     27   360
 7    34     28   339
 8    32     25   329
 9    28     24   218
10    35     24   402
11    33     26   342
12    28     25   205
13    32     23   368
14    25     22   196
15    28     21   304
16    30     23   294
17    29     23   275
18    32     25   336
19    34     26   384
20    35     27   368

 
今回の分析で行うのは、次の重回帰式の { \displaystyle \beta} (ベータ)を推定することです。

{ \displaystyle
客数 = \beta_{1} + \beta_{2}*最高気温 + \beta_{3}*最低気温
}

 
推定する { \displaystyle \beta} は3つあるので、この数を p として、Stan で分析するデータをリスト形式で dat にセットします。Stan で分析するデータはリスト形式はお約束(たぶん)。

dat <- list(
    N = nrow(icecream2),  #サンプルサイズ、20。つまりデータフレーム icecream2 の行数。
    saiko_kion = icecream2$saiko, #最高気温。
    saitei_kion = icecream2$saitei, #最低気温。
    kyaku = icecream2$kyaku, #客数。
    p = 3 #推定するベータの数。今回は3。
    )

 
次はStan コード。単回帰分析の拡張版なので、難しくありません。icecream2.stan というファイル名で R の作業ディレクトリに保存します。作業ディレクトリの場所は getwd() でわかります。

data {
   int<lower=0> N; #サンプルサイズの数。
   real<lower=0> kyaku[N]; #継続したお客さんの数は20回。
   real<lower=0>  saiko_kion[N]; #計測した最高気温は20個。
   real<lower=0>  saitei_kion[N]; #計測した最低気温は20個。
   int<lower=0> p; #切片および偏回帰係数の合計の数。
}
 
parameters {
   real beta[p]; #切片beta1、最高気温の偏回帰係数beta2、最低気温の偏回帰係数beta3。
   real<lower=0> sigma; #客数が従う正規分布の標準偏差。
}
 
transformed parameters  {
   real mu[N]; #mu は data と parameters セクションで定義した変数を使って定義する。
   for (i in 1:N) {
     mu[i] <- beta[1] + beta[2]*saiko_kion[i] + beta[3]*saitei_kion[i]; #重回帰の式をここで書いて mu に代入。
   }
}
 
model { 
     kyaku ~ normal(mu, sigma);  #客数は正規分布に従う。
}

 
あとは分析実行。私の環境では35秒ほどかかりました。結果もいっしょに出力します。サンプリング回数などは『基礎からベイズ統計学』の第7章章末問題そのままです。適切なサンプリング回数などはまだよくわからないんですよね…。

fit<-stan(file = scr, model_name = scr, data = dat, chains = cha, warmup = war, iter=ite)
print(fit, pars = c("beta"))

 
結果は次のとおりです。上が Stan によるベイズ推定、下が一般の重回帰分析 lm(kyaku ~ saiko + saitei, data = icecream2) の結果です。

          mean se_mean    sd    2.5%     25%    50%    75%  97.5% n_eff Rhat
beta[1] -90.57    0.21 41.39 -172.32 -117.40 -90.53 -63.87  -8.07 39158    1
beta[2]  25.96    0.01  1.63   22.72   24.90  25.96  27.01  29.20 33787    1
beta[3] -16.71    0.01  2.21  -21.07  -18.13 -16.71 -15.29 -12.30 33320    1
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -90.640     37.718  -2.403   0.0279 *  
saiko         25.957      1.486  17.463 2.71e-12 ***
saitei       -16.703      2.012  -8.301 2.20e-07 ***

 
だいたい同じになりました。めでたし、めでたし。今回参考にしたページは次のところです。ありがとうございました。

RPubs - Multiple linear regression with Stan

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

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