中古マンション売買データを分析してみた(2.2) - BUGSのモデルをStanに置き換えてみた
ごぶさたしています.シリーズの途中にも関わらず,1ヶ月半ほどブログを放置していました.仕事にかまけて更新が億劫になったのと,Pythonでせこせこbot作ったりしてたのが主な要因ではありますが.気がついたら第3回BUGS/Stan勉強会が開催されていて*1,そうだそうだいい加減やらなきゃということで,再会しました.ちょうど先日開発用にMac miniを購入したこともあり,i7クアッドコアでなら計算も捗ること間違いなしです.
今回も引き続きテクニカルな話になるので,Stanに興味がある人のみお読みください.
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以降に書いてあるのですが,いちおう自分でも試してみたいと思います.