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)
y t = F t t + v t, v t N (0, V t ) t = G t t 1 + w t, w t N (0, W t ) t = 1,..., n 0 N (m 0, C 0 )
## build function! BuildLLM <- function(theta) {! dlmmodpoly(order = 1,! dv = theta[1],! dw = theta[2])! }
##! 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)
#! x <- matrix(c(rep(0, 27),! rep(1, length(nile) - 27)),! ncol = 1)
##! 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)! }
##! fit.reg <- dlmmle(nile,! parm = rep(0, 2),! build = BuildReg)! model.reg <- BuildReg(fit.reg$par)! smooth.reg <- dlmsmooth(nile,! mod = model.reg)
y t = Z t t + t, t N (0, H t ) t+1 = T t t + R t t, t N (0, Q t ) t = 1,..., n 1 N (a 1, P 1 )
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")
fit.van <- fitssm(inits = c(-4, -7, 2),! model = model.van,! method = BFGS")!! pred.van <- predict(fit.van$model,! states = 1:2)
set.seed(1234)! n.t <- 50 #! N.lat <- rep(50, n.t) #! p <- 0.7 #! N.obs <- rbinom(n.t, N.lat, p) #!
var! N, #! y[n], #! y_hat[n], #! lambda[n], # log(y_hat)! p, #! tau, sigma;
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);! }
inits <- list()! inits[[1]] <- list(p = 0.9, sigma = 1,! lambda = rep(log(max(n.obs) + 1), n.t))! inits[[2]] <- list(p = 0.7, sigma = 3,! lambda = rep(log(max(n.obs) + 1), n.t))! inits[[3]] <- list(p = 0.8, sigma = 5,! lambda = 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)!
http://mc-stan.org/
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;! G[1, 1] <- 1;! m0[1] <- 0;! C0[1, 1] <- 1.0e+6;!
parameters {! real<lower=0> sigma[2];! }! transformed parameters {! vector[1] V;! cov_matrix[1] W;!! V[1] <- sigma[1] * sigma[1];! W[1, 1] <- sigma[2] * sigma[2];! }!
model {! y ~ gaussian_dlm_obs(f, G, V, W, m0, C0);! sigma ~ uniform(0, 1.0e+6);! }
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)
traceplot(fit, pars = "sigma", inc_warmup = FALSE)
> 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).!
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)
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]);! }!! //! theta[1] ~ normal(0, 1.0e+4);! sigma ~ uniform(0, 1.0e+6);! }