StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

不等間隔の状態空間モデル

Stan

日付単位とかでデータを取ることが多いこのご時世、等間隔の状態空間モデルを使うことが多いと思います。しかし、ふと不等間隔の状態空間モデルってどうやるんだろーとつぶやいたところ、ご指導いただきました。いつも大変感謝です。

今回は不等間隔の時系列データに対し、ざっくり小数点第1位で等間隔のメッシュを作って、データのない所は欠損値として扱う、通常の状態空間モデルに落とす方法でやってみました。Stanコードは以下になります。

data {
   int<lower=1> N;
   int<lower=1> N_obs;
   int<lower=1, upper=N> Idx_obs[N_obs];
   real Y_obs[N_obs];
}

parameters {
   real r[N];
   real<lower=0> s_r;
   real<lower=0> s_y;
}

model {
   for (n in 3:N)
      r[n] ~ normal(2*r[n-1] - r[n-2], s_r);
   for(i in 1:N_obs)
      Y_obs[i] ~ normal(r[Idx_obs[i]], s_y);
}
  • 9行目:背後のrはすべてのメッシュで存在します。
  • 16行目:今回は2階差分のCAR modelを使っています。
  • 3～5行目:あとのRコードの方で詳しく見ますが、観測したところのデータとメッシュのインデックスを渡します。
  • 17～18行目:観測したデータのところの尤度を計算しています。

キックするRコードは以下になります。

library(rstan)

d <- read.delim("input/GPbook-Fig2_05.data.txt", sep="\t")
X_obs <- d$X
X_mesh <- -7.0 + 0:140 * 0.1
N <- length(X_mesh)
N_obs <- nrow(d)
X_mod <- round(X_obs, 1)
Idx_obs <- as.integer(round((X_mod + 7.0)/0.1))
Y_obs <- d$Y
data <- list(N=N, N_obs=N_obs, Idx_obs=Idx_obs, Y_obs=Y_obs)

stan_fit <- stan(file="model/model.stan", chains=0)
fit <- stan(fit=stan_fit,　data=data,　iter=1000, warmup=400, thin=1,　seed=123,　chains=3)
  • 5,8,9行目:ここでX_obsを小数点第1位のデータX_modに落として、そのインデックス（X_modがメッシュで左から何番目の位置か）を計算しています。なお計算時間はSurface Pro 3（core i5）で2次のCAR modelが1chainあたり20秒ぐらいでした。意外と早い。

結果は以下になります。想定よりうまくいっているのではないかと思います。

f:id:StatModeling:20201107072907p:plain

白抜き点はデータ、青線はrのMCMCサンプルの中央値、薄い灰帯は同じく90%信用区間、濃い灰帯は同じく50%信用区間です。

f:id:StatModeling:20201107072902p:plain

こちらは、灰帯だけ変えて、1800個のMCMCサンプルを全てプロットしたものです。

最後に上の図を描いたRコードを載せておきます。

library(rstan)
library(ggplot2)
library(reshape2)

...(Stanで計算)...

d_obs <- data.frame(X=X_obs, Y=Y_obs)

la <- extract(fit)
N_mcmc <- length(la$lp__)

qua <- apply(la$r, 2, quantile, prob=c(0.05, 0.25, 0.5, 0.75, 0.95))
d_qua <- data.frame(X=X_mesh, t(qua))
colnames(d_qua) <- c('X', 'p5', 'p25', 'p50', 'p75', 'p95')
p <- ggplot()
p <- p + geom_ribbon(data=d_qua, aes(x=X, y=p50, ymax=p95, ymin=p5), alpha=1/4)
p <- p + geom_ribbon(data=d_qua, aes(x=X, y=p50, ymax=p75, ymin=p25), alpha=2/4)
p <- p + geom_line(data=d_qua, aes(x=X, y=p50), color='blue')
p <- p + geom_point(data=d_obs, aes(x=X, y=Y), shape=as.factor(1), size=2)
p <- p + coord_cartesian(ylim=c(-3.5, 3.5))
p <- p + xlab('X') + ylab('Y')
p <- p + theme(legend.position='none')
ggsave(file='output/r_est_quantile.png', plot=p, dpi=300, width=4, height=3)


d_mcmc <- data.frame(mcmc_sample=1:N_mcmc, la$r)
colnames(d_mcmc) <- c('mcmc_sample', X_mesh)
d_melt <- melt(d_mcmc, id='mcmc_sample', variable.name='X', value.name='Y')
d_melt$X <- as.numeric(as.character(d_melt$X))
p <- ggplot()
p <- p + geom_line(data=d_melt, aes(x=X, y=Y, group=mcmc_sample), color='darkorange3', alpha=1/100)
p <- p + geom_line(data=d_qua, aes(x=X, y=p50), color='blue')
p <- p + geom_point(data=d_obs, aes(x=X, y=Y), shape=as.factor(1), size=2)
p <- p + coord_cartesian(ylim=c(-3.5, 3.5))
p <- p + xlab('X') + ylab('Y')
p <- p + theme(legend.position='none')
ggsave(file='output/r_est_MCMC.png', plot=p, dpi=300, width=4 , height=3)