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

About connecting the dots.

statistics/machine learning adversaria.

中古マンション売買データを分析してみた(2.3) - Stanの行列演算速度を軽く試す

とりあえずStanになれるというのも兼ねて,軽くマニュアル読みながら*1,どういう書き方が良い書き方なのかを学ぶことにしました.そうしたら,マニュアルのp141 Vectorizationのところで,基本的にvectorやmatrix使った方が計算早いよって書いてあるので,実際にコード書いてみて試してみました.

サンプルコード

ベタな回帰式の書き方

data {
    int<lower=1> N;
    int<lower=1, upper=5> d_d[N];
    int<lower=1900, upper=2014> d_f[N];
    int<lower=1, upper=20> d_r[N];
    int<lower=1, upper=20> d_s[N];
    int<lower=0, upper=200> d_p[N];
}
parameters {
    real a;
    real b_d;
    real b_f;
    real b_r;
    real b_s;
    real b_p;
    real<lower=0> s;
}
model {
    for (i in 1:N)
        d_p[i] ~ normal(a + b_d*d_d[i] + b_f*d_f[i] + b_r*d_r[i] + b_s*d_s[i], s);
}

マトリックス演算

data {
    int<lower=1> N;
    int<lower=1> M;
    matrix[N, M] X;
    vector[N]    Y;
}
parameters {
    real          alpha;
    vector[M]     beta;
    real<lower=0> sigma;
}
model {
    for(i in 1:N)
        Y[i] ~ normal(alpha + X[i] * beta, sigma);
}

Rコード

キックする側のコードはこちら

# load library
library('rstan')

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

# store elapsed time
time_matrix = array()
time_linear = array()

# matrix operation
X = t(rbind(distance, from, room, space))
Y = price
d.stan = list(N=nrow(X), M=ncol(X), X=X, Y=Y)
## matrix operation
fit.matrix = stan(file="script/mcmc_basic_lm_vector.stan",
                  data=d.stan,
                  iter=10,
                  chains=1)
for (i in 1:100) {
  time_matrix[i] = system.time(
    basic_lm.fit<-stan(fit=fit.matrix,
                       data=d.stan,
                       iter=100,
                       chains=1,
                       thin=10,
                       warmup=20)  
  )
}

# usual calculation using linear expression
d.stan = list(N=nrow(d),
              d_d=d$distance,
              d_f=d$from,
              d_r=d$room,
              d_s=d$space,
              d_p=d$price)
fit.linear<-stan(file="script/mcmc_basic_lm.stan",
                 data=d.stan,
                 iter=10,
                 chains=1)
for (i in 1:100) {
  time_linear[i] = system.time(
    basic_lm.fit<-stan(fit=fit.linear,
                       data=d.stan,
                       iter=100,
                       chains=1,
                       thin=10,
                       warmup=20)  
  )
}

結果

結論からいうと,全くそんなことはなかったので,どういうことなのか未だに良くわかっていないんですよね... まぁある程度変数が増えると,マトリックス使わないと変数管理が面倒でやってられない感じもするので,基本的にはマトリックスでいった方がいいのかなぁとは思っていますが...

> summary(time_matrix)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.526   1.542   2.222  42.330  40.490 254.600 
> summary(time_linear)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.219   0.614   0.933  22.770  46.710 112.900 

追記

このエントリについて@beroberoさんから,以下のようなご指摘をいただきました.おそらく,この程度の行列演算だと高速化の恩恵はほとんどなくて,むしろ通常計算時に値の範囲を絞っているほうが効果がある可能性はたしかにありそうかなぁと思われます.

そして,何より階層ベイズで個体差を考えた瞬間に,そもそも行列演算はできないので,結局このエントリの趣旨とは何だったのか状態だったり...