1. No category

2014年3月16日 第61回日本生態学会大会（広島）
T13 生態学における状態空間モデルの利用
状態空間モデルの
実行方法と
実行環境の比較
森林総合研究所
伊東宏樹
本日とりあつかうソフトウェア
•
Rパッケージ
• dlm
• KFAS
•
MCMC
• BUGS言語
• Stan
サンプルコードなどの置き場所
http://www001.upp.so-net.ne.jp/ito-hi/stat/2014ESJ/
Statistical
Software for
State Space
Models
Commandeur et al. (2011)
Journal of Statistical
Software 41(1)
State Space
Models in R
Petris & Petrone (2011)
Journal of Statistical
Software 41(4)
dlm
dlm
•
Dynamic Linear Model （動的線形モデル）
•
線形+正規分布
•
カルマンフィルタ
•
パラメータ推定
•
最尤推定／ベイズ推定
dlmの記法
データモデル
yt = F t ✓ t + vt ,
プロセスモデル
✓t = Gt ✓t
1
vt ⇠ N (0, Vt )
+ wt , wt ⇠ N (0, Wt )
t = 1, ... , n
✓0 ⇠ N (m0 , C0 )
ナイル川の流量の変化
data(Nile)
dlmによるLocal Level Model
Petris and Petrone (2011)より
## build functionの定義!
BuildLLM <- function(theta) {!
dlmModPoly(order = 1,!
dV = theta[1],!
dW = theta[2])!
}
このような関数を定義しておく。
dlmによるLocal Level Model
## パラメーターの最尤推定!
fit.llm <- dlmMLE(Nile, parm = c(100, 2),!
build = BuildLLM,!
lower = rep(1e-4, 2))!
!
## 推定したパラメーターをbuild functionで使用!
model.llm <- BuildLLM(fit.llm$par)!
!
## 平滑化!
smooth.llm <- dlmSmooth(Nile, model.llm)
平滑化
dlmSmooth()
アスワンダム着工
ナイル川の流量の変化
data(Nile)
dlmによる回帰モデル
# アスワンダム着工の前後を変数に!
x <- matrix(c(rep(0, 27),!
rep(1, length(Nile) - 27)),!
ncol = 1)
dlmによる回帰モデル
## モデル定義!
model.reg <- dlmModReg(x, dW = c(1, 0))!
BuildReg <- function(theta) {!
V(model.reg) <- exp(theta[1])!
diag(W(model.reg))[1] <- exp(theta[2])!
return(model.reg)!
}
dlmによる回帰モデル
## 最尤推定!
fit.reg <- dlmMLE(Nile,!
parm = rep(0, 2),!
build = BuildReg)!
model.reg <- BuildReg(fit.reg$par)!
smooth.reg <- dlmSmooth(Nile,!
mod = model.reg)
アスワンダム着工
ナイル川の流量の変化
data(Nile)
dlmの文献
•
Petris G, Petrone S, Campagnoli (2009)
Dynamic Linear Models with R Springer
•
•
和合肇（監訳）・萩原淳一郎（訳）(2013)「R
によるベイジアン動的線形モデル」朝倉書店
Petris G (2010) An R package for dynamic
linear models. Journal of Statistical
Software 36(12)
KFAS
KFAS
•
Kalman Filter and Smoother for Exponential
Family State Space Models
•
正規分布以外の分布（ポアソン分布など）を扱
える
•
最尤推定
KFASの記法
データモデル
yt = Z t ↵ t + ✏ t ,
プロセスモデル
✏t ⇠ N (0, Ht )
↵t+1 = Tt ↵t + Rt ⌘t , ⌘t ⇠ N (0, Qt )
t = 1, ..., n
↵1 ⇠ N (a1 , P1 )
イギリスのバン運転手の死者・重傷者数
data(Seatbelts)
KFASによるポアソン分布の状態空間モデル
help(KFAS)より
model.van <- SSModel(VanKilled ~ law +!
SSMtrend(degree = 1,!
Q = list(matrix(NA))) +!
SSMseasonal(period = 12,!
sea.type = “dummy",!
Q = matrix(NA)),!
data = Seatbelts,!
distribution = "poisson")
KFASによるポアソン分布の
状態空間モデル
fit.van <- fitSSM(inits = c(-4, -7, 2),!
model = model.van,!
method = “BFGS")!
!
pred.van <- predict(fit.van$model,!
states = 1:2)
lawとSSMtrend()のみをつかう
シートベルト着用義務化
季節変化をのぞいた予測値
BUGS
WinBUGS, OpenBUGS, JAGS
BUGS
•
MCMCによるベイズ推定
•
柔軟なモデリング
•
Rパッケージでは対応できないモデル
例題
•
ある生物の個体数を推定する。
•
一定の発見確率にしたがって発見される。
Kéry & Schaub (2011) Bayesian Population
Analysis using WinBUGS: A hierarchical
perspective Chapter 5を参考にした。
データ生成
set.seed(1234)!
n.t <- 50
N.lat <- rep(50, n.t)
p <- 0.7
N.obs <- rbinom(n.t, N.lat, p)
#
#
#
#
観察回数!
真の個体数!
発見確率!
観察個体数!
真の個体数
観測された個体数
生成されたデータ
Binomial(50, 0.7)
BUGSモデル
var!
N,
y[N],
y_hat[N],
lambda[N],
p,
tau, sigma;
#
#
#
#
#
観察回数!
観察された個体数!
「真の個体数」の推定値!
log(y_hat)!
発見確率!
BUGSモデル
model {!
## データモデル!
for (t in 1:N) {!
y[t] ~ dbin(p, y_hat[t]);!
y_hat[t] <- trunc(exp(lambda[t]));!
}!
## プロセスモデル!
for (t in 2:N) {!
lambda[t] ~ dnorm(lambda[t - 1], tau);!
}!
## 事前分布!
lambda[1] ~ dnorm(0, 1.0E-4);!
p ~ dbeta(2, 2);!
sigma ~ dunif(0, 100);!
tau <- 1 / (sigma * sigma);!
}
JAGSによる実行
inits <- list()!
inits[[1]] <- list(p = 0.9,
lambda =
inits[[2]] <- list(p = 0.7,
lambda =
inits[[3]] <- list(p = 0.8,
lambda =
sigma = 1,!
rep(log(max(N.obs) + 1), n.t))!
sigma = 3,!
rep(log(max(N.obs) + 1), n.t))!
sigma = 5,!
rep(log(max(N.obs) + 1), n.t))!
!
model <- jags.model("ks51.bug.txt",!
data = list(N = n.t, y = N.obs),!
inits = inits, n.chains = 3,!
n.adapt = 100000)!
samp <- coda.samples(model,!
variable.names = c("y_hat", “sigma",!
"p"),!
n.iter = 3000000, thin = 3000)!
真の個体数
「真の個体数」の推定値
観測された個体数
推定結果
Stan
http://mc-stan.org/
Stan
•
•
MCMCによるベイズ推定
•
Hamiltonian Monte Carlo (HMC)
•
No U-Turn Sampling (NUTS)
Stan → C++ → ネイティブバイナリ
Stan
•
CmdStan
• コマンドラインから
•
RStan
• Rから
•
PyStan
• Pythonから
StanによるDLM
data(Nile)を使用
StanによるDLM
データ
data {!
int<lower=0> N;!
matrix[1, N] y;!
}!
transformed data {!
matrix[1, 1] F;!
matrix[1, 1] G;!
vector[1]
m0;!
cov_matrix[1] C0;!
!
}
F[1, 1] <- 1;!
dlmと同様の
G[1, 1] <- 1;!
データを用意
m0[1] <- 0;!
C0[1, 1] <- 1.0e+6;!
StanによるDLM
パラメータ
parameters {!
real<lower=0> sigma[2];!
}!
transformed parameters {!
vector[1]
V;!
cov_matrix[1] W;!
!
dlmと同様の
V[1] <- sigma[1] * sigma[1];!
W[1, 1] <- sigma[2] * sigma[2];!
パラメータを
}!
用意
StanによるDLM
モデル
model {!
y ~ gaussian_dlm_obs(F, G, V, W, m0, C0);!
sigma ~ uniform(0, 1.0e+6);!
}
StanによるDLM
library(rstan)!
!
model <- stan("kalman.stan",!
data = list(y = matrix(c(Nile),!
nrow = 1),!
N = length(Nile)),!
pars = c("sigma"),!
chains = 3,!
iter = 1500, warmup = 500,!
thin = 1)
MCMCの軌跡
traceplot(fit, pars = "sigma", inc_warmup = FALSE)
StanによるDLM
> print(fit)!
Inference for Stan model: kalman.!
3 chains, each with iter=1500; warmup=500; thin=1; !
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] 121.2
0.5 13.8
92.6 112.7 121.5 130.3 148.4
889
1!
sigma[2]
45.5
0.6 17.6
18.3
32.7
43.2
55.7
85.2
833
1!
lp__
-541.6
0.0 1.1 -544.6 -542.0 -541.3 -540.9 -540.6
904
1!
!
Samples were drawn using NUTS(diag_e) at Sun Feb 9 06:06:42 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).!
StanによるDLM
ベイズ推定されたパラメータをdlmで使用
sigma <- apply(extract(fit, "sigma")$sigma, 2, mean)!
!
library(dlm)!
!
buildNile <- function(theta) {!
dlmModPoly(order = 1, dV = theta[1], dW = theta[2])!
}!
modNile <- buildNile(sigma^2)!
smoothNile <- dlmSmooth(Nile, modNile)
平滑化
Stanでベイズ推定されたパラメータをdlmで使用
Stanによる状態空間モデルの解析
•
gaussian_dlm_obs()でうまくいかないことも
•
自分でモデルを記述することも当然可能
Stanによる状態空間モデルの解析
data {!
int<lower=0>
real
}!
parameters {!
real
real<lower=0>
}!
N;!
y[N];!
theta[N];!
sigma[2];!
Stanによる状態空間モデルの解析
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]);!
}!
// 事前分布!
theta[1] ~ normal(0, 1.0e+4);!
sigma ~ uniform(0, 1.0e+6);!
まとめ
状態空間モデルをあつかえるソフトウェア
•
Rパッケージ: dlm, KFAS
•
•
•
関数に与える引数の意味を理解する。
ベイズ推定: BUGS, Stan
•
柔軟なモデリングが可能。
•
計算時間はかかる。
上記以外のソフトウェアもある。