RStan:ローカルレベルモデル-パラメーター既知

1. ストーリーとシナリオ

シナリオ：『健康志向のAさん、日々の体重測定日記』

健康を気遣うAさんは、毎朝同じ時間に体重を測り、記録することにしました。しかし、体重は日々の食事や体内の水分量によって細かく変動するため、測定した数値がそのまま「真の体重」の増減を表しているわけではないことに気づきます。

ここで、Aさんの体重推移を状態空間モデルで考えてみましょう。

• 状態 (潜在変数) α_t：時刻 t における「真の体重」
  • これは、Aさんの基礎的な体重のことで、日々の食生活や運動習慣によって少しずつ、しかし連続的に変化していくと考えられます。この変化は、前日の「真の体重」 α_{t-1} を中心としたランダムな動き（ランダムウォーク）で表現できます。この変動の大きさをシステムノイズと呼びます。
• 観測値 y_t：時刻 t における「観測体重」
  • これは、Aさんが毎朝体重計で測定する実際の数値です。この数値は、その日の「真の体重」 α_t に、体内の水分量の一時的な増減や、体重計の測定誤差などが加わったものと考えられます。この誤差の大きさを観測ノイズと呼びます。

モデルの数式表現

このシナリオは、典型的なローカルレベルモデルで表現できます。

1. 状態方程式（真の体重の推移）

   α_t = α_{t-1} + w_t , w_t ~ Normal(0, σ_w^2)

   • 今日の真の体重(α_t)は、昨日の真の体重(α_{t-1})に、日々の生活習慣による変動(w_t)が加わったもの。
   • w_tは、平均0、分散σ_w^2の正規分布に従う（システムノイズ）。

2. 観測方程式（体重計で測る値）

   y_t = α_t + v_t , v_t ~ Normal(0, σ_v^2)

   • 今日測定した体重(y_t)は、今日の真の体重(α_t)に、測定時の誤差(v_t)が加わったもの。
   • v_tは、平均0、分散σ_v^2の正規分布に従う（観測ノイズ）。

シミュレーションのパラメータ設定（既知とする値）

今回はシミュレーションなので、これらのパラメータを事前に決めておきます。

• 観測期間 T: 100日間
• 初期の真の体重 α_1: 65.0 kg
• システムノイズの標準偏差 σ_w: 0.1 kg
  • 真の体重が1日で変動する幅。比較的小さな値で、ゆっくりとした変化を表します。
• 観測ノイズの標準偏差 σ_v: 0.5 kg
  • 測定ごとのばらつき。水分量などで日々の測定値が大きくぶれることを表現するため、システムノイズより大きな値にします。

この設定で、「真の体重」がどのように推移し、それに対して「観測体重」がどのように得られるのかをシミュレートします。

2. RStanによるシミュレーションの実装

それでは、上記のシナリオをRとRStanで実装します。

RStanでパラメータを推定するのではなく、generated quantitiesブロックを使って、モデルからデータを生成（シミュレート）します。

library(rstan)

sapply(X = c("rstan"), packageVersion)

stan_output <- "D:/stan_output"

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

$rstan
[1]  2 32  7

stan_code <- "
data {
  int<lower=1> T;          // データ期間
  real alpha1;             // 初期状態
  real<lower=0> sigma_w;   // システムノイズの標準偏差
  real<lower=0> sigma_v;   // 観測ノイズの標準偏差
}

parameters {
  // パラメータを推定しないので、このブロックは空
}

model {
  // モデルの尤度計算も不要なので、このブロックも空
}

generated quantities {
  vector[T] alpha; // 真の体重（状態）
  vector[T] y;     // 観測体重（観測値）

  // t=1（初日）の状態
  alpha[1] = alpha1;

  // t=2以降の状態を生成
  for (t in 2:T) {
    alpha[t] = normal_rng(alpha[t-1], sigma_w);
  }

  // 各時点の観測値を生成
  for (t in 1:T) {
    y[t] = normal_rng(alpha[t], sigma_v);
  }
}
"

ステップ3：Rでのシミュレーション実行と結果の可視化

RスクリプトでStanコードを呼び出し、シミュレーションを実行します。

seed <- 20250626

# 1. パラメータの設定
T <- 100
alpha1 <- 65.0
sigma_w <- 0.1
sigma_v <- 0.5

# 2. Stanに渡すためのデータリストを作成
data_list <- list(
  T = T,
  alpha1 = alpha1,
  sigma_w = sigma_w,
  sigma_v = sigma_v
)

# 3. Stanでシミュレーションを実行
# algorithm="Fixed_param" を使うと、サンプリングを行わず
# generated quantitiesブロックのみを実行するため高速です。
sim_fit <- stan(
  model_code = stan_code,
  data = data_list,
  algorithm = "Fixed_param", # シミュレーションのための指定
  chains = 1,
  iter = 1,
  seed = seed
)

# stanfit オブジェクトの保存
setwd(stan_output)
saveRDS(object = sim_fit, file = "stan_fit.rds")

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

# 4. シミュレーション結果の抽出
sim_data <- rstan::extract(sim_fit)
alpha_sim <- as.vector(sim_data$alpha)
y_sim <- as.vector(sim_data$y)

# 5. 結果をデータフレームにまとめる
df_sim <- data.frame(
  time = 1:T,
  true_weight = alpha_sim,
  observed_weight = y_sim
)

# 6. 結果の可視化
library(ggplot2)
ggplot(df_sim, aes(x = time)) +
  geom_line(aes(y = true_weight, color = "真の体重 (α_t)"), linewidth = 1) +
  geom_point(aes(y = observed_weight, color = "観測体重 (y_t)"), size = 2, alpha = 0.7) +
  scale_color_manual(
    name = "系列",
    values = c("真の体重 (α_t)" = "darkred", "観測体重 (y_t)" = "steelblue")
  ) +
  labs(
    title = "ローカルレベルモデルによる体重データのシミュレーション",
    subtitle = "Aさんの100日間の体重測定日記",
    x = "日数",
    y = "体重 (kg)"
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

実行結果と解釈

上記のコードを実行すると、以下のようなプロットが生成されます。

• 赤い線（真の体重 α_t）: Aさんの「真の体重」の推移を表します。σ_w = 0.1と小さく設定したため、滑らかにゆっくりと変動（ランダムウォーク）しているのがわかります。
• 青い点（観測体重 y_t）: Aさんが毎日体重計で測った値です。「真の体重」の周りに、σ_v = 0.5という比較的大きなばらつきを持ってプロットされています。日々の測定値だけを見ると体重が大きく増減しているように見えますが、その背後には滑らかに変動する「真の体重」がある、という状態空間モデルの考え方が視覚的に理解できます。

以上です。