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

About connecting the dots.

statistics/machine learning adversaria.

中古マンション売買データを分析してみた(2.2) - BUGSのモデルをStanに置き換えてみた

ごぶさたしています.シリーズの途中にも関わらず,1ヶ月半ほどブログを放置していました.仕事にかまけて更新が億劫になったのと,Pythonでせこせこbot作ったりしてたのが主な要因ではありますが.気がついたら第3回BUGS/Stan勉強会が開催されていて*1,そうだそうだいい加減やらなきゃということで,再会しました.ちょうど先日開発用にMac miniを購入したこともあり,i7クアッドコアでなら計算も捗ること間違いなしです.

今回も引き続きテクニカルな話になるので,Stanに興味がある人のみお読みください.

BUGSからStanへ

Stanの導入については,TJOさんの記事など参考にインストールすればよいのではないでしょうか.私の環境はMacなので,公式のガイドをみながら普通にインストールしました*2

Stanコード

とりあえず今回はコードを動かしてみようということで,同じくTJOさんの記事を参考に,前回記事のlm() in BUGSをそのままStanに移植してみました.特にひねりもないので,Stanコードと,それをkickするRコードは以下の通りです.

data {
    int<lower=1> N;
    int<lower=1> M;
    matrix[N, M] X;
    int<lower=0, upper=200> Y[N];
}
parameters {
    real alpha;
    vector[M] beta;
    real<lower=0> sigma;
}
model {
    for (i in 1:N)
        Y[i] ~ normal(alpha + dot_product(beta, X[i]), sigma);
}
# load library
library('ggplot2')
library('rstan')

# load data
d <- read.delim('data/mantions.csv', header=T, sep=',')
d = na.omit(d)
attach(d)

# package data for stan
X = t(rbind(distance, from, room, space))
Y = price
d.stan = list(N=nrow(X), M=ncol(X), X=X, Y=Y)
basic_lm.fit<-stan(file="script/mcmc_basic_lm.stan",
                   data=d.stan,
                   iter=10000,
                   chains=3)

結果

ということで,結果は以下の通りです.実施条件は,10000 iteration, 5000 burn-in, 1 thinでした*3.lm()の結果とほぼ同じ値が得られていることがわかるかと思います.

> basic_lm.fit
Inference for Stan model: mcmc_basic_lm.
3 chains, each with iter=10000; warmup=5000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=15000.

             mean se_mean    sd      2.5%       25%       50%       75%     97.5% n_eff Rhat
alpha    -1873.43    0.36 26.97  -1926.05  -1891.81  -1873.25  -1855.43  -1820.64  5767    1
beta[1]     -4.23    0.00  0.12     -4.47     -4.31     -4.23     -4.15     -4.00  9278    1
beta[2]      0.96    0.00  0.01      0.94      0.95      0.96      0.97      0.99  5765    1
beta[3]     -2.60    0.00  0.25     -3.10     -2.77     -2.60     -2.43     -2.11  7327    1
beta[4]      2.79    0.00  0.20      2.40      2.66      2.79      2.92      3.18  6920    1
sigma       10.07    0.00  0.10      9.87     10.00     10.07     10.14     10.27 10643    1
lp__    -13553.97    0.02  1.73 -13558.15 -13554.89 -13553.67 -13552.71 -13551.58  5047    1

Samples were drawn using NUTS(diag_e) at Sun Jul 13 20:00:27 2014.
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).

計算速度について

何も考えずにデフォルトの値を弄らず実行したら,今回の計算はやたらに時間がかかってしまいました*4.また,その際にいくつかの計算式の書き方を試していたら,微妙に時間のかかり方が違っていたため,その辺りを検証するコードを書いて,改めてエントリとしてあげたいと思います*5.あと,せっかくクアッドコアになったので,beroberoさんの記事を参考に,3 chainをパラで回せるようにしたいとも思ってます.これもまた改めて...

SAMPLING FOR MODEL 'mcmc_basic_lm' NOW (CHAIN 1).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 5701.48 seconds (Warm-up)
#                7073.33 seconds (Sampling)
#                12774.8 seconds (Total)


SAMPLING FOR MODEL 'mcmc_basic_lm' NOW (CHAIN 2).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 5417.52 seconds (Warm-up)
#                7073.67 seconds (Sampling)
#                12491.2 seconds (Total)


SAMPLING FOR MODEL 'mcmc_basic_lm' NOW (CHAIN 3).

Iteration:    1 / 10000 [  0%]  (Warmup)
Iteration: 1000 / 10000 [ 10%]  (Warmup)
Iteration: 2000 / 10000 [ 20%]  (Warmup)
Iteration: 3000 / 10000 [ 30%]  (Warmup)
Iteration: 4000 / 10000 [ 40%]  (Warmup)
Iteration: 5000 / 10000 [ 50%]  (Warmup)
Iteration: 5001 / 10000 [ 50%]  (Sampling)
Iteration: 6000 / 10000 [ 60%]  (Sampling)
Iteration: 7000 / 10000 [ 70%]  (Sampling)
Iteration: 8000 / 10000 [ 80%]  (Sampling)
Iteration: 9000 / 10000 [ 90%]  (Sampling)
Iteration: 10000 / 10000 [100%]  (Sampling)
#  Elapsed Time: 6178.64 seconds (Warm-up)
#                6914.37 seconds (Sampling)
#                13093 seconds (Total)

*1:勉強会の名前とは裏腹に,ほぼStanの話だけだったと評判ですが.まぁBUGSは開発が止まっているみたいなので,仕方はないかもしれないですが...

*2:ちなみに,R3.0.2+Xcode+Commandline toolsでRStanを入れようとしたら,g++あたりでエラーを吐いてしまったので,Google groupのディスカッションをみて,R3.1.1の最新版にあげたら正常にインストールできました.

*3:何も考えずにデフォルトの値を変えずに実施したら,1 thinだったため非常に実行に時間がかかってしまいました...

*4:おかしーなー,Stanの長所は収束が早いことのはずなのに,こんなに計算自体が重たいんじゃ結局意味ないじゃんー,とか勝手に勘違いしていたことは内緒です.

*5:具体的には,arrayを使うか,vectorを使うか,matrixを使うかという選択の問題です.正解自体は,Stan manualのp141以降に書いてあるのですが,いちおう自分でも試してみたいと思います.