RStan:ローカルレベルモデル+カルマンフィルタ-パラメーター未知

# データの生成
set.seed(seed)
T <- 150 # 観測期間 (日数)
mu_zero <- 5.0 # 初期水位 (メートル)

# パラメータの真の値 (未知として推定する対象)
sigma_w_true <- 0.2 # システムノイズの標準偏差 (日々の水位変動の大きさ)
sigma_v_true <- 0.4 # 観測ノイズの標準偏差 (センサーの測定誤差)

# 状態(mu)と観測値(y)を生成
mu_true <- numeric(T)
y <- numeric(T)

# t=1 の状態
w <- rnorm(1, mean = 0, sd = sigma_w_true)
mu_true[1] <- mu_zero + w
v <- rnorm(1, mean = 0, sd = sigma_v_true)
y[1] <- mu_true[1] + v

# t=2以降の状態
for (t in 2:T) {
  w <- rnorm(1, mean = 0, sd = sigma_w_true)
  mu_true[t] <- mu_true[t - 1] + w
  v <- rnorm(1, mean = 0, sd = sigma_v_true)
  y[t] <- mu_true[t] + v
}

# 生成したデータの可視化
true_data <- tibble::tibble(
  time = 1:T,
  mu_true = mu_true,
  y = y
)

library(ggplot2)
ggplot(true_data, aes(x = time)) +
  geom_line(aes(y = mu_true, color = "真の水位 (μ_t)"), linewidth = 1) +
  geom_point(aes(y = y, color = "観測値 (y_t)"), alpha = 0.7) +
  labs(
    title = "シミュレーションデータ",
    subtitle = "真の水位と観測値の推移",
    x = "日数", y = "水位 (m)", color = "データ系列"
  ) +
  scale_color_manual(values = c("真の水位 (μ_t)" = "dodgerblue", "観測値 (y_t)" = "black")) +
  theme_bw() +
  theme(legend.position = "bottom")

# stanfit オブジェクトの読み込み
setwd(stan_output)
fit <- readRDS("stan_fit.rds")

# 収束診断 (Rhatが1.1未満であれば概ね良好)
print(fit, pars = c("mu_zero", "sigma_v", "sigma_w"))

# MCMCサンプルを抽出
samples <- rstan::extract(fit)

# パラメータの事後分布をプロット
p1 <- ggplot(tibble::tibble(sigma_v = samples$sigma_v), aes(x = sigma_v)) +
  geom_density(fill = "skyblue", alpha = 0.7) +
  geom_vline(xintercept = sigma_v_true, color = "red", linetype = "dashed", linewidth = 1) +
  labs(
    title = "観測ノイズ (sigma_v) の事後分布",
    subtitle = paste0("真値 = ", sigma_v_true), x = "sigma_v"
  ) +
  theme_bw()

p2 <- ggplot(tibble::tibble(sigma_w = samples$sigma_w), aes(x = sigma_w)) +
  geom_density(fill = "lightgreen", alpha = 0.7) +
  geom_vline(xintercept = sigma_w_true, color = "red", linetype = "dashed", linewidth = 1) +
  labs(
    title = "システムノイズ (sigma_w) の事後分布",
    subtitle = paste0("真値 = ", sigma_w_true), x = "sigma_w"
  ) +
  theme_bw()

# パラメータのプロットを並べて表示(Figure 2)
library(patchwork)
p1 + p2

# トレースプロットを表示(Figure 3)
bayesplot::mcmc_trace(fit, pars = c("sigma_v", "sigma_w"))

library(dplyr)

# 推定された状態(mu_filt)の要約統計量を取得
mu_summary <- summary(fit, pars = "mu_filt")$summary %>%
  as_tibble() %>%
  mutate(time = 1:T)

# 可視化
ggplot(true_data, aes(x = time)) +
  # 95%信用区間
  geom_ribbon(data = mu_summary, aes(ymin = `2.5%`, ymax = `97.5%`), fill = "skyblue", alpha = 0.8) +
  # 真の水位
  geom_line(aes(y = mu_true, color = "真の水位 (μ_t)"), linewidth = 1) +
  # 観測値
  geom_point(aes(y = y, color = "観測値 (y_t)"), alpha = 0.6) +
  # 状態の事後中央値
  geom_line(data = mu_summary, aes(y = `50%`, color = "推定された水位 (μ_t)"), linewidth = 1) +
  labs(
    title = "状態推定の結果",
    subtitle = "カルマンフィルタによる水位の推定値と95%信用区間",
    x = "日数", y = "水位 (m)", color = "データ系列"
  ) +
  scale_color_manual(values = c(
    "真の水位 (μ_t)" = "black",
    "観測値 (y_t)" = "gray50",
    "推定された水位 (μ_t)" = "dodgerblue"
  )) +
  theme_bw() +
  theme(legend.position = "bottom")