以下の書籍
を中心に時系列解析を勉強していきます。
16. 一般状態空間モデルの逐次解法
16.2 粒子フィルタによる状態の推定
ここでは自作で粒子フィルタによる状態の推定を行なう。
library(rstan) # 人工データの生成 library(dlm) # 前処理 set.seed(4521) W <- 1 V <- 2 m0 <- 10 C0 <- 9 mod <- dlmModPoly(order = 1, dW = W, dV = V, m0 = m0, C0 = C0) t_max <- 200 sim_data <- dlmForecast(mod = mod, nAhead = t_max, sampleNew = 1) y <- sim_data$newObs[[1]] # 粒子数 n <- 10000 # データの整形 y <- c(NA_real_, y) # 事前分布を表すダミー分を先頭に追加 # リサンプリング用の系列は全時点で保存しておく k <- matrix(1:n, nrow = n, ncol = t_max+1) # 事前分布の設定 # 粒子(実現値) x <- matrix(NA_real_, nrow = t_max+1, ncol = n) x[1,] <- rnorm(n, mean = mod$m0, sd = sqrt(mod$C0)) # 粒子(重み) w <- matrix(NA_real_, nrow = t_max+1, ncol = n) w[1,] <- 1/n # 時間軸方向の処理 for (t in (1:t_max)+1){ # 状態方程式：粒子(実現値)を生成 x[t,] <- rnorm(n, mean = x[t-1,], sd = sqrt(mod$W)) # 観測方程式：粒子(重み)を更新 w[t,] <- w[t-1,] * dnorm(y[t], mean = x[t,], sd = sqrt(mod$V)) # 重みの規格化 w[t,] <- w[t,] / sum(w[t,]) ### リサンプリング # リサンプリング用のインデックス列 k[,t] <- sample(1:n, prob = w[t,], replace = T, size = n) # 粒子(実現値)：リサンプリング用のインデックス列を新たな番号とする x[t,] <- x[t, k[,t]] # 粒子(重み)：リセット w[t,] <- 1/n } # 結果の整形：事前分布の分を除去するなど y <- ts(y[-1]) k <- k[,-1, drop = F] x <- x[-1,,drop = F] w <- w[-1,,drop = F] # 平均・25%値・75%値を求める scratch_m <- sapply(1:t_max, function(t){ mean(x[t,]) }) scratch_m_quant <- lapply(c(0.05,0.95),function(quant){ sapply(1:t_max, function(t){ quantile(x[t,], probs = quant) }) }) df <- data.frame(t = 1:t_max, var = "観測値", val = as.numeric(y)) # yはTimeSeriesなので型変換 df <- rbind(df,data.frame(t = 1:t_max, var = "平均", val = scratch_m)) df <- rbind(df,data.frame(t = 1:t_max, var = "5%点", val = scratch_m_quant[[1]])) df <- rbind(df,data.frame(t = 1:t_max, var = "95%点", val = scratch_m_quant[[2]])) df$var <- factor(df$var,levels = c("観測値","平均","5%点","95%点")) g1 <- ggplot(df, aes(x = t, y = val, color = var, linetype = var)) g1 <- g1 + geom_line(size = 0.6) g1 <- g1 + theme_bw() g1 <- g1 + labs(title = "粒子フィルタリングによるモデリング", x = "時点",y = "y") g1 <- g1 + theme(axis.text= element_text(size = 12), axis.text.x=element_text(size = 12), legend.title = element_blank(), legend.text = element_text(size = 12), legend.position = "top", axis.title = element_text(size = 12), axis.title.y = element_text(angle = 90,size = 12), plot.title= element_text(size = 16,hjust = 0.5)) plot(g1) ########################## # 予測を実施 # 未来時点のデータ領域を追加 future_time <- 20 x <- rbind(x, matrix(NA_real_, nrow = future_time, ncol = n)) w <- rbind(w, matrix(NA_real_, nrow = future_time, ncol = n)) # 時間軸方向の処理 for(t in t_max + (1:20)){ # 状態方程式：粒子(実現値)を生成 x[t,] <- rnorm(n, mean = x[t-1,], sd = sqrt(mod$W)) # 粒子(重みを更新) w[t,] <- w[t-1,] } scratch_a <- sapply(t_max+(1:future_time), function(t){ mean(x[t,]) }) scratch_a_quant <- lapply(c(0.05,0.95),function(quant){ sapply(t_max+(1:future_time), function(t){ quantile(x[t,], probs = quant) }) }) df2 <- data.frame(t = 1:t_max, var = "観測値", val = as.numeric(y)) # yはTimeSeriesなので型変換 df2 <- rbind(df2, data.frame(t =t_max+(0:future_time), var = "予測(平均)", val = c(as.numeric(y)[length(y)],scratch_a))) df2 <- rbind(df2, data.frame(t =t_max+(0:future_time), var = "予測(5%点)", val = c(as.numeric(y)[length(y)],scratch_a_quant[[1]]))) df2 <- rbind(df2, data.frame(t =t_max+(0:future_time), var = "予測(95%点)", val = c(as.numeric(y)[length(y)],scratch_a_quant[[2]]))) df2$var <- factor(df2$var,levels = c("観測値","予測(平均)","予測(5%点)","予測(95%点)")) g2 <- ggplot(df2, aes(x = t, y = val, color = var, linetype = var)) g2 <- g2 + geom_line(size = 0.6) g2 <- g2 + theme_bw() g2 <- g2 + labs(title = "粒子フィルタリングによるモデリング(予測値)", x = "時点",y = "y") g2 <- g2 + theme(axis.text= element_text(size = 12), axis.text.x=element_text(size = 12), legend.title = element_blank(), legend.text = element_text(size = 12), legend.position = "top", axis.title = element_text(size = 12), axis.title.y = element_text(angle = 90,size = 12), plot.title= element_text(size = 16,hjust = 0.5)) plot(g2)
推定 |
予測 |
 |
 |
