中古マンション売買データを分析してみた(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)