【R】DLM for Nile via RStan【demo】

◆NileデータにDLM
====== dlm.stan =====
data {
  int<lower=0>  N;
  real          y[N];
}
parameters {
  real          theta[N];
  real<lower=0> sigma[2];
}
model {
  //観測モデル
  for (t in 1:N) {
    y[t] ~ normal(theta[t], sigma[1]);
  }
  //状態モデル
  for (t in 2:N) {
    theta[t] ~ normal(theta[t - 1], sigma[2]);
  }
}
=====

===== dlm_stan.R =====
library(rstan)
data(Nile)
fit <- stan("dlm.stan",
            data = list(y = as.vector(Nile),
                        N = length(Nile)),
            pars = c("sigma", "theta"),
            chains = 3,
            iter = 1500,
            warmup = 500,
            thin = 1)

tmp <- apply(as.matrix(fit), 2, mean)
theta <- tmp[grep(pattern = "theta", x = names(tmp))]
s <- start(Nile)[1]
e <- end(Nile)[1]
plot(Nile, type="o")
lines(s:e, theta, col="red", lwd=2)
=====

◆状態の初期値を考慮する
===== dlm_m0.stan =====
data {
  int<lower=0>  N;
  real          y[N];
}
parameters {
  real          theta[N];
  real<lower=0> sigma[2];
  real          m0;
}
model {
  //観測モデル
  for (t in 1:N) {
    y[t] ~ normal(theta[t], sigma[1]);
  }
  //状態モデル
  theta[1] ~ normal(m0, sigma[2]);
  for (t in 2:N) {
    theta[t] ~ normal(theta[t - 1], sigma[2]);
  }
}
=====

===== dlm_m0_stan.R =====
library(rstan)
data(Nile)
fit <- stan("dlm_m0.stan",
            data = list(y = as.vector(Nile),
                        N = length(Nile)),
            pars = c("sigma", "theta", "m0"),
            chains = 3,
            iter = 1500,
            warmup = 500,
            thin = 1)

tmp <- apply(as.matrix(fit), 2, mean)
theta <- tmp[grep(pattern = "theta", x = names(tmp))]
s <- start(Nile)[1]
e <- end(Nile)[1]
plot(Nile, type="o")
lines(s:e, theta, col="red", lwd=2)
=====

◆ダム建設を考慮する
===== dlm_dam.stan =====
data {
  int<lower=0>  N;
  real          y[N];
  int<lower=0, upper=1> x[N];
}
parameters {
  real          theta[N];
  real<lower=0> sigma[2];
  real          lambda;
  real          m0;
}
model {
  //観測モデル
  for (t in 1:N) {
    y[t] ~ normal(theta[t] + lambda * x[t], sigma[1]);
  }
  //状態モデル
  theta[1] ~ normal(m0, sigma[2]);
  for (t in 2:N) {
    theta[t] ~ normal(theta[t - 1], sigma[2]);
  }
}
=====
 
===== dlm_dam_stan.R =====
library(rstan)
data(Nile)

dat <- list(y = as.vector(Nile),
            N = length(Nile),
            x = as.vector(c(rep(0, 27), rep(1, length(Nile)-27))))

fit <- stan("dlm_dam.stan",
            data = dat,
            pars = c("sigma", "theta", "lambda", "m0"),
            chains = 3,
            iter = 1500,
            warmup = 500,
            thin = 1)

tmp <- apply(as.matrix(fit), 2, mean)
theta <- tmp[grep(pattern = "theta", x = names(tmp))]
lambda <- tmp[grep(pattern = "lambda", x = names(tmp))]

s <- start(Nile)[1]
e <- end(Nile)[1]

plot(Nile, type="o")
lines(s:e, theta+lambda*dat$x, col="red", lwd=2)
=====

【R】RStanのインストール

=====
Sys.setenv(MAKEFLAGS = "-j4")
source('http://mc-stan.org/rstan/install.R', echo = TRUE, max.deparse.length = 2000)
install_rstan()
library(rstan)
=====

◆参考
How to Install RStan

【Python】numpyのインストール for Mac

Macでnumpyを使いたい

sudo pip install numpy
とすると、
RuntimeError: Broken toolchain: cannot link a simple C program
と言われるので、

sudo ARCHFLAGS=-Wno-error=unused-command-line-argument-hard-error-in-future pip install --upgrade numpy
とする。


※参考
Marvericksにnumpyをインストールする