dlmのGibbsサンプラー、それからStan

Rのdlmパッケージには、Gibbsサンプラーでパラメータ推定するdlmGibbsDIG()もある。

Rによるベイジアン動的線形モデル (統計ライブラリー) の4章を参考にする。スペインの投資額のデータ。

##
## Petris at al. (2009) Dynamic Linear Models with R
## Chapter 4
##
library(dlm)

data.file <- url("http://definetti.uark.edu/~gpetris/dlm/Datasets/invest2.dat")
# data.file <- "Datasets/invest2.dat"
invSpain <- ts(read.table(data.file,
                          colClasses = "numeric")[, 2] / 1000,
               start = 1960)
plot(invSpain, las = 1, xlab = "Year", ylab = "Investiments (×1000)")

MCMC

set.seed(5672)
MCMC <- 12000
gibbsOut <- dlmGibbsDIG(invSpain,             # data vector
                        mod = dlmModPoly(2),  # dlm
                        a.y = 1,              # prior mean of observation precision
                        b.y = 1000,           # prior variance of observation precision
                        a.theta = 10,         # prior mean of system precisions
                        b.theta = 1000,       # prior variance of system precisions
                        n.sample = MCMC,
                        thin = 1,
                        save.states = TRUE)

burn <- 2000
theta <- mcmcMean(cbind(gibbsOut$dV[-(1:burn)],
                        gibbsOut$dW[-(1:burn), ]))

計算はあまりはやくない。結果。

> theta
                W.1        W.2   
  0.012197   0.117391   0.329588 
 (0.000743) (0.007682) (0.007833)

推定された状態をプロット。

y <- apply(gibbsOut$theta[, 1, -(1:burn)], 1, mean)
plot(invSpain, las = 1, xlab = "Year", ylab = "Investiments (×1000)")
lines(start(invSpain)[1]:end(invSpain)[1], dropFirst(y), col = 2, lty = 2)

観測誤差がちいさい。（Investmentのつづりをまちがえているが、きにしない。）

カルマンフィルタにかけてみる。

build <- function(theta) {
  dlmModPoly(order = 2, dV = theta[1], dW = theta[2:3])
}
mod <- build(theta[1, ])
ks <- dlmSmooth(invSpain, mod)
plot(invSpain, las = 1,
     xlab = "Year", ylab = "Investiments (×1000)")
lines(start(invSpain)[1]:end(invSpain)[1],
      dropFirst(ks$s[, 1]), col = "red", lty = 2)

MCMCの部分をStanでやってみる。

library(rstan)

model.text <- "
  data {
  int<lower=1>  N;
  real          y[N];
}
parameters {
  real          mu[N];
  real          mu0;
  real          beta[N];
  real          beta0;
  real<lower=0> sigma[3];
}
model {
  for (t in 1:N) {
    y[t] ~ normal(mu[t], sigma[1]);
  }
  mu[1] ~ normal(mu0, sigma[2]);
  beta[1] ~ normal(beta0, sigma[3]);
  for (t in 2:N) {
    mu[t] ~ normal(mu[t - 1] + beta[t - 1], sigma[2]);
    beta[t] ~ normal(beta[t - 1], sigma[3]);
  }
  sigma[1] ~ uniform(0, 1000);
  sigma[2] ~ uniform(0, 1000);
  sigma[3] ~ uniform(0, 1000);
}
"

model <- stan_model(model_code = model.text)
fit <- sampling(model, data = list(N = length(invSpain),
                                   y = c(invSpain)),
                chains = 3, iter = 6000, warmup = 1000, thin = 5)

結果

> print(fit, pars = "sigma", digits = 2)
Inference for Stan model: model.text.
3 chains, each with iter=6000; warmup=1000; thin=5; 
post-warmup draws per chain=1000, total post-warmup draws=3000.

         mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
sigma[1] 0.11    0.01 0.07 0.01 0.05 0.10 0.16  0.27   136 1.03
sigma[2] 0.30    0.01 0.17 0.05 0.16 0.28 0.41  0.66   303 1.01
sigma[3] 0.57    0.01 0.15 0.22 0.48 0.58 0.66  0.83   474 1.00

Samples were drawn using NUTS(diag_e) at Sat Jan 18 13:08:20 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).

カルマンフィルタへ。

sigma <- apply(extract(fit, "sigma")$sigma, 2, mean)
mod <- build(sigma)
ks <- dlmSmooth(invSpain, mod)
plot(invSpain, las = 1, xlab = "Year", ylab = "Investiments (×1000)")
lines(start(invSpain)[1]:end(invSpain)[1],
      dropFirst(ks$s[, 1]), col = "red")

Stanのモデルだが、increment_log_prob()をつかおうと、下のようにかいたら、Error in function stan::prob::gaussian_dlm_obs_log(l): columns of Fcolumns of F (2) and rows of y (1) must match in sizeRejecting proposed initial value with zero density.というエラーになった。Fは1行2列のはずなのだが。

[追記] Fは転置したものがincrement_log_prob()関数に渡されるので、2行1列が正しい。下のコードは修正ずみ。

data {
  int<lower=0>  N;
  matrix[1, N]  y;
  matrix[2, 1]  F;
  matrix[2, 2]  G;
  vector[2]     m0;
  cov_matrix[2] C0;
}
parameters {
  real<lower=0> sigma[3];
}
transformed parameters {
  vector[1]     V;
  cov_matrix[2] W;

  V[1] <- sigma[1];
  W[1, 1] <- sigma[2];
  W[2, 2] <- sigma[3];
  W[1, 2] <- 0;
  W[2, 1] <- 0;
}
model {
  increment_log_prob(gaussian_dlm_obs_log(y, F, G, V, W, m0, C0));
  sigma ~ uniform(0, 1.0e+3);
}

[追記] 上のコードだが、VとWには分散がはいる。sigmaは分散とおもってください。

このモデルを使用するRStanのコード

model <- stan_model(model_code = model.text)
fit <- sampling(model, data = list(N = length(invSpain),
                                   y = matrix(c(invSpain), nrow = 1),
                                   F = matrix(c(1, 0), nrow = 2, ncol = 1),
                                   G = matrix(c(1, 0, 1, 1), nrow = 2, ncol = 2),
                                   m0 = c(10, 10),
                                   C0 = matrix(c(10, 0, 0, 10), nrow =2, ncol = 2)),
                chains = 3, iter = 6000, warmup = 1000, thin = 5)

タグ：R STAn