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

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Stan コードをベクトル化したら速くなるか:AR(1) モデルを例に

ちょっと速くなる。

こういうモデルのパラメータを推定しました。

 X_t = \beta_0 + \beta_1 X_{t-1}+\varepsilon_t
\varepsilon_t は平均 0、分散 \sigma^2正規分布
\beta_0, \beta_1 の事前分布は一様分布です。

ベクトル化前のコードはこちら。

// AR1_1.stan
data {
  int<lower=1> T;
  real y[T];
}
parameters {
  real beta0;
  real beta1;
  real<lower=0> sigma;
}
model{
  for(t in 2:T){
   y[t] ~ normal(beta0+beta1*y[t-1],sigma); 
  }
}

ベクトル化後のコードはこちらです。

//AR1_2.stan
data {
  int<lower=1> T;
  vector[T] y;
}
parameters {
  real beta0;
  real beta1;
  real<lower=0> sigma;
}
model{
    y[2:T] ~ normal(rep_vector(beta0,T-1)+beta1*y[1:(T-1)],sigma); 
}

R でシミュレーションデータを生成して、パラメータを推定します。

set.seed(10)
y <- numeric(10000)
y[1] <- 2
for(i in 2:10000)
  y[i]<-rnorm(1,0.5+0.8*y[i-1])
dat<-list(T=10000,y=y)
plot(y,type="l")
library(rstan)
fit1 <-stan("AR1_1.stan",data=dat,chain=1)
#8.14929 seconds (Total)
fit2 <-stan("AR1_2.stan",data=dat,chain=1)
#6.44406 seconds (Total)

計算させたPCはMacBook Airで、
プロセッサ:2.2 GHz Intel Core i7
メモリ:8 GB
です。