はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

どうでもいい話ですが、この記事は私にとって、記念すべきPositronを使って書いた初めてのアウトプットです！たぶんまだRStudioから完全に卒業することはないと思いますが、Positronもかなり使いやすかったので、気になっている方はぜひ一度試してみてください！

さて、皆さんもベイズ統計の教科書などで、ベイズの強みとして、事前分布の設計による情報の借用やシュリンケージ（縮約）の話を聞いたことがあると思います。

でも、はっきり言って、私自身も最初は「なにそれ？それの何が嬉しいの？」と思っていましたw

そこで本記事では、時系列で平均値を推定するタスクを例に、情報の借用が実際にどのような場面で役に立つのかをシミュレーションで確かめてみます。

想定するのは、普段は十分な量のデータが取得できているのに、インフラの障害や運用ルールの変更など何らかの原因によって、ある期間だけデータが極端に少なくなってしまったというシナリオです。

このような状況でも、時系列の構造を事前分布に組み込んだ階層ベイズモデルであれば、前後の時点から情報を借りることで、わずかなデータからでも真の値に近い推定ができるのでしょうか？

本記事では、時系列の依存関係を組み込んだ階層ベイズモデルと、各時点をより独立に扱うシンプルな階層ベイズモデルを比較します。

さらに、「たまたま1回うまくいっただけでは？」という疑問にも答えるため、同じ実験を1,000回シミュレーションします。

その結果、時系列の依存関係を組み込んだモデルは、シンプルなモデルに対して73%のシミュレーションでRMSEが小さくなることを確認できました。

なお、本記事で扱うシナリオには重要な前提があります。

ここでは、データが取得できなくなること自体と、推定したい値との間に系統的な関係がないことを想定しています。

例えば不動産業界で、これまではネット上に掲載されている物件データをほぼ全件取得できていたものの、担当者の長期休暇や一時的なシステム障害などによって、ある期間だけランダムに5%程度しか取得できなくなった、という状況は今回のシナリオに対応します。

一方で、例えば業績が悪化した物件や企業ほどデータを公開しなくなる、あるいは「見られたくないデータ」を意図的に取得・公開しないといった場合には、欠けているデータそのものに情報が含まれています。このようなケースでは、単純に前後の時点から情報を借りればよいわけではなく、データが欠ける仕組みそのものを考慮したモデル化が必要になります。

したがって、本記事で示す方法は「データが少なくなったらベイズで何でも解決できる」というものではありません。あくまで、データ量の急減が推定対象そのものと系統的に結びついていない状況で、残された少量のデータと時系列構造をどこまで活用できるかを検証するものです。

さて、「情報の借用」とは結局何が嬉しいのか、実際にデータが足りなくなったときにどのような力を発揮するのかを、Stanのコードとシミュレーション結果を見ながら確認していきましょう。

シナリオ

全体像

まずは、今回検証するシナリオを設定します。

ある指標の平均値を時系列で継続的に推定している状況を考えてみましょう。例えば、毎月取得している商品の価格、求人の給与、物件の賃料などをイメージしてください。

今回は500時点分のデータを用意し、それぞれの時点に「本当の平均値」が存在するとします。ただし、この平均値は時点ごとに完全に独立して変化するのではなく、前の時点の値とある程度似た値を取りながら、徐々に変化していくものとします。

つまり、

今月の平均値は、先月の平均値とまったく無関係ではない

という、時系列データでは比較的自然な状況を想定します。

数式でいうと、$t$期の平均$\mu_{t}$はこのように生成されたとします：

$$
\mu_{t} \sim Normal(\mu_{t-1}, \rho)
$$

さらにt期のi番目の観測値$y_{t,i}$が

$$
y_{t,i} \sim Normal(\mu_{t}, \sigma)
$$

でサンプリングされます。

普段は500件、しかし突然1件しか取れなくなる😨

通常時には、それぞれの時点について 500件 のデータを観測できるものとします。

ところが、200〜250時点目だけ、担当者の長期休暇など何らかの理由によって、取得できるデータが突然 1件まで減少したとします。

期間	1時点あたりのデータ数
1〜199	500件
200〜250	1件 😨
251〜500	500件

かなり極端な設定ですが、情報の借用の効果をわかりやすく確認するため、あえてこのような状況を作ります。

普段500件のデータから平均値を推定していたのに、突然1件しか取得できなくなったらどうなるでしょうか。

当然、その1件のデータだけで結果を報告すると、たまたま大きな値や小さな値が観測されただけで、推定結果も大きく振り回されてしまいます。

しかし、ここで一つ使えそうな情報があります：前後の時点の平均値です。

今回のシミュレーションでは、本当の平均値が時間とともに徐々に変化するようにデータを生成しています。そのため、例えば220時点目のデータが1件しかなくても、219時点目や218時点目などの情報は、220時点目の平均値を考えるうえで参考になるはずです。

そこで今回は、次の2つのモデルを比較します。

時点間の依存関係を考慮しないモデル（independent）
各時点の平均値を共通の分布から生成されたものとして扱いますが、隣接する時点同士の関係は明示的には利用しません。
時点間の依存関係を考慮するモデル（dependent）
「今回の平均値は前回の平均値と似ている」という構造を事前分布に組み込み、前の時点から情報を借りられるようにします。

ポイントは、どちらも階層ベイズモデルであることです。

「ベイズ vs. 非ベイズ」を比較したいわけではありません。今回確認したいのは、データの構造に合わせて事前分布を設計し、適切な場所から情報の借用に意味があるのかという点です。

「本当の平均値」がわかる世界を作る

実際のビジネスデータでは、推定したい「本当の平均値」はもちろんわかりません。

そのため、推定結果が滑らかになったとしても、

なんとなくきれいな推定結果になっただけでは？

という疑問が残ります。

そこで今回は、まず「本当の平均値」をこちらで生成し、その値を使って観測データをシミュレーションします。

具体的には、500時点にわたってランダムウォークする潜在的な平均値を生成し、そこにランダムなノイズを加えて観測データを作ります。

coef_master <- 500 |>
  seq_len() |>
  tibble::tibble(
    time = _
  ) |>
  dplyr::mutate(
    beta_init = cumsum(rnorm(dplyr::n()))
  ) |>
  dplyr::mutate(
    # モデル推定の都合上、平均値がゼロを中心に分散が1になるように標準化
    beta = (beta_init - mean(beta_init)) / sd(beta_init)
  ) |>
  dplyr::select(!beta_init)

そして通常時には各時点から500件、200〜250時点目だけは1件の観測値を生成します。

sim_data <- 500 |>
  seq_len() |>
  purrr::map(
    \(i){
      if (dplyr::between(i, 200, 250)){
        rep(i, 1)
      } else {
        rep(i, 500)
      }
    }
  ) |>
  unlist() |>
  tibble::tibble(
    time = _
  ) |>
  dplyr::left_join(
    coef_master, by = "time"
  ) |>
  dplyr::mutate(
    e = rnorm(dplyr::n()),
    y = beta + e
  )

まず、平均値とデータの分布を可視化しましょう：

g_y <- sim_data |> 
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = time, y = y), color = ggplot2::alpha("blue", 0.1))

g_beta <- coef_master |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = time, y = beta))

gridExtra::grid.arrange(g_y, g_beta, nrow = 1)

左側の図でわかるように、データの急な減少があまりにもひどすぎるwww

また、平均値 $\beta$（beta）の変動は時系列的な相関の強い株価のような形をしていることがわかります。

これで、

本当の平均値 $\beta$
実際に観測される値 $y$
どの期間でデータが極端に少なくなったか

のすべてを把握した状態で、モデルの性能を検証できます。

つまり今回は、単に「推定結果がそれっぽく見えるか」ではなく、推定された平均値が、こちらがあらかじめ知っている正解 $\beta$ にどれだけ近づけるかを直接評価できます。

そして1回のシミュレーションだけでは「たまたま今回のデータでうまくいっただけ」という可能性も残ります。

そこで最後には、データ生成からモデル推定までを1,000回繰り返し、どちらのモデルがより安定して正解に近づけるのかも検証します。

果たして、データが 500件からわずか1件まで減ってしまったとき、前後の時点から情報を借りることでどこまで推定を救えるのでしょうか。

モデル設計

時点間の依存関係を考慮しないモデル

まず、時点間の依存関係を考慮しないモデルからはじめましょう。

まず、各時点の平均値の全体的なばらつきパラメータと平均値をサンプリングします：

$$
\rho \sim InvGamma(1, 1)
$$

$$
\beta \sim Normal(0, \rho)
$$

次に、観測値の分散と切片をサンプリングし：

$$
\sigma \sim InvGamma(1, 1)
$$

$$
\alpha \sim Normal(0, 5)
$$

最後に観測値をサンプリングします：

$$
y_{i} \sim Normal(\alpha + \beta_{t_{i}}, \sigma)
$$

Stan言語で記述するとこうなります：

time_series_independent.stan

data {
  int time_type;

  int N;
  array[N] int time;
  array[N] real y;
}
parameters {
  real<lower=0> rho;
  sum_to_zero_vector[time_type] beta;
  real<lower=0> sigma;
  real intercept;
}
model {
  rho ~ inv_gamma(1, 1);
  
  beta ~ normal(0, rho);
  sigma ~ inv_gamma(1, 1);
  intercept ~ normal(0, 5);
  y ~ normal(intercept + beta[time], sigma);
}

時点間の依存関係を考慮するモデル

次に、本記事の主役である、時点間の依存関係を考慮するモデルを紹介します。

まず、一期目の平均値のばらつきと一期目以降の平均値のブレ幅パラメータをサンプリングします：

$$
\rho_{1} \sim InvGamma(1, 1)
$$

$$
\rho_{2} \sim InvGamma(1, 1)
$$

次に、一期目の平均値をサンプリングします：

$$
\beta_{1} \sim Normal(0, \rho_{1})
$$

次に、t期目の平均値はこのようにサンプリングされます：

$$
\beta_{t} \sim Normal(\beta_{t - 1}, \rho_{2})
$$

残りは時点間の依存関係を考慮しないモデルと一緒ですが、観測値の分散と切片をサンプリングし：

$$
\sigma \sim InvGamma(1, 1)
$$

$$
\alpha \sim Normal(0, 5)
$$

最後に観測値をサンプリングします：

$$
y_{i} \sim Normal(\alpha + \beta_{t_{i}}, \sigma)
$$

Stan言語で記述するとこうなります：

time_series_dependent.stan

data {
  int time_type;

  int N;
  array[N] int time;
  array[N] real y;
}
parameters {
  vector<lower=0>[2] rho;
  sum_to_zero_vector[time_type] beta;
  real<lower=0> sigma;
  real intercept;
}
model {
  rho ~ inv_gamma(1, 1);
  
  beta[1] ~ normal(0, rho[1]);
  for (i in 2:time_type){
    beta[i] ~ normal(beta[i - 1], rho[2]);
  }
  sigma ~ inv_gamma(1, 1);
  intercept ~ normal(0, 5);
  y ~ normal(intercept + beta[time], sigma);
}

モデル推定結果

では早速モデル推定を実施しましょう：

> m_d_init <- cmdstanr::cmdstan_model("time_series_dependent.stan") 
> m_d_estimate <- m_d_init$variational(
   seed = 1,
   data = list(
     time_type = nrow(coef_master),
 
     N = nrow(sim_data),
     time = sim_data$time,
     y = sim_data$y
   )
 )
> m_d_summary <- m_d_estimate$summary() 
> m_id_init <- cmdstanr::cmdstan_model("time_series_independent.stan") 
> m_id_estimate <- m_id_init$variational(
   seed = 1,
   data = list(
     time_type = nrow(coef_master),
 
     N = nrow(sim_data),
     time = sim_data$time,
     y = sim_data$y
   )
 )
> m_id_summary <- m_id_estimate$summary()
Model executable is up to date!
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.004311 seconds 
1000 transitions using 10 leapfrog steps per transition would take 43.11 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Success! Found best value [eta = 1] earlier than expected. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100      -343235.013             1.000            1.000 
   200      -320768.300             0.535            1.000 
   300      -320367.410             0.357            0.070 
   400      -320051.585             0.268            0.070 
   500      -320191.943             0.215            0.001   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  4.2 seconds.
Model executable is up to date!
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.00448 seconds 
1000 transitions using 10 leapfrog steps per transition would take 44.8 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Success! Found best value [eta = 1] earlier than expected. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100      -501449.802             1.000            1.000 
   200      -325824.844             0.770            1.000 
   300      -321836.324             0.517            0.539 
   400      -320611.530             0.389            0.539 
   500      -320829.530             0.311            0.012 
   600      -320616.121             0.259            0.012 
   700      -320529.419             0.222            0.004   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  5.2 seconds.

Positronは、コードをペーストして実行した際に、複数行のコードをまとめて実行するなど、RStudioとは挙動が若干異なるように感じました。

推定結果を見てみる

結果をこのように可視化します：

g_d_series <- m_d_summary |>
  dplyr::filter(stringr::str_detect(variable, "beta")) |>
  dplyr::mutate(
    time = dplyr::row_number()
  ) |>
  dplyr::bind_cols(
    answer = coef_master$beta
  ) |>
  dplyr::mutate(
    status = dplyr::case_when(
      dplyr::between(time, 200, 250) ~ "few data",
      TRUE ~ "normal"
    )
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = time, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_point(ggplot2::aes(x = time, y = answer, color = status), alpha = 0.3) + 
  ggplot2::geom_vline(xintercept = c(200, 250), linetype = "dashed", color = "red", linewidth = 1) +
  ggplot2::labs(
    title = "dependent"
  )

g_id_series <- m_id_summary |>
  dplyr::filter(stringr::str_detect(variable, "beta")) |>
  dplyr::mutate(
    time = dplyr::row_number()
  ) |>
  dplyr::bind_cols(
    answer = coef_master$beta
  ) |>
  dplyr::mutate(
    status = dplyr::case_when(
      dplyr::between(time, 200, 250) ~ "few data",
      TRUE ~ "normal"
    )
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = time, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::geom_point(ggplot2::aes(x = time, y = answer, color = status), alpha = 0.3) + 
  ggplot2::geom_vline(xintercept = c(200, 250), linetype = "dashed", color = "red", linewidth = 1) +
  ggplot2::labs(
    title = "independent"
  )

g_d_error <- m_d_summary |>
  dplyr::filter(stringr::str_detect(variable, "beta")) |>
  dplyr::mutate(
    time = dplyr::row_number()
  ) |>
  dplyr::bind_cols(
    answer = coef_master$beta
  ) |>
  dplyr::mutate(
    error = answer - mean,
    status = dplyr::case_when(
      dplyr::between(time, 200, 250) ~ "few data",
      TRUE ~ "normal"
    )
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = time, y = error, color = status)) + 
  ggplot2::geom_vline(xintercept = c(200, 250), linetype = "dashed", color = "red", linewidth = 1) +
  ggplot2::ylim(-1.5, 1.5) +
  ggplot2::labs(
    title = "dependent"
  )

g_id_error <- m_id_summary |>
  dplyr::filter(stringr::str_detect(variable, "beta")) |>
  dplyr::mutate(
    time = dplyr::row_number()
  ) |>
  dplyr::bind_cols(
    answer = coef_master$beta
  ) |>
  dplyr::mutate(
    error = answer - mean,
    status = dplyr::case_when(
      dplyr::between(time, 200, 250) ~ "few data",
      TRUE ~ "normal"
    )
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = time, y = error, color = status)) + 
  ggplot2::geom_vline(xintercept = c(200, 250), linetype = "dashed", color = "red", linewidth = 1) +
  ggplot2::ylim(-1.5, 1.5) +
  ggplot2::labs(
    title = "independent"
  )

library(patchwork)

(g_d_series | g_d_error)/(g_id_series | g_id_error)

では、結果を見てみましょう。

上段が時点間の依存関係を考慮したモデル （dependent）、下段が各時点を独立に扱ったモデル （independent）です。

また、左側の図では黒線がモデルによる beta の推定値、点がシミュレーション時に生成した「正解」の beta を表しています。赤い点と2本の赤い破線で囲まれた200〜250時点目が、データ数を500件からわずか1件まで減らした期間です。

右側には、それぞれの時点について

正解 - 推定値

をプロットしています。そのため、0に近いほど正解に近く、0から離れるほど推定を外していると解釈できます。

データが500件ある期間では、どちらもかなり当たっている

まず注目したいのは、データが通常どおり500件存在する期間です。

この期間では、dependent・independentのどちらのモデルでも、推定値は正解にかなり近い値になっています。右側の誤差を見ても、200〜250時点目以外ではほとんど0付近に集中しています。

これはある意味当然です。

500件ものデータがあれば、その時点のデータだけでも平均値についてかなり多くの情報を得られます。そのため、わざわざ前後の時点から情報を借りなくても、どちらのモデルもそれなりに正確な推定ができます。

今回確認したいのは、ここではありません。

500件から1件になると、モデルの違いが一気に現れる

問題は、赤い破線で囲まれた200〜250時点目です。

この期間では、それまで500件あったデータをわずか1件まで減らしています。

下段のindependentモデルを見ると、この瞬間から推定値が大きく暴れ始めています。

これは、各時点を独立に扱っているためです。

例えば220時点目に1件しかデータがなければ、219時点目や221時点目にどれだけ情報が存在していても、それを220時点目の推定に直接利用することができません。その結果、たまたま観測された1件の値に推定結果が強く左右されます。

右下の誤差を見ると、この違いはさらに明確です。

通常期間では0付近に集中していた誤差が、データが1件しかない期間では一気に大きくなり、正解から大きく外れる時点が多数発生しています。

一方、dependentモデルは前の時点から情報を借りられる

それに対して、上段のdependentモデルでは様子がかなり異なります。

こちらのモデルでは、

beta[i] ~ normal(beta[i - 1], rho[2]);

とすることで、

今回の beta は、前回の beta とある程度似ているはず

という時系列の構造を事前分布に組み込んでいます。

そのため、ある時点のデータが1件しかなくなったとしても、その1件だけを頼りに推定する必要はありません。

前の時点までに得られた情報を利用しながら、現在の beta を推定できます。

実際、左上を見ると、200〜250時点目でも推定値が正解の動きに比較的よく追従しています。

右上の誤差を見ても、データ不足期間では通常期間より誤差が大きくなっているものの、independentモデルほど大きくは暴れていません。

これが、今回見たかった情報の借用の効果です。

データが大量に存在するときには、それぞれの時点のデータだけでも十分に推定できます。しかし、何らかのトラブルによってデータが急激に少なくなったとき、時系列の構造を事前分布に組み込んでおけば、データが豊富だった時点から情報を借りることで、少量のデータしかない時点の推定を補うことができます。

でも、きれいになっただけでは？

ここまでの結果を見ると、

やっぱり時系列構造を入れた階層ベイズのほうが強い！

と言いたくなります。

しかし、まだそう結論づけるのは早いです。

今回見ているのは、set.seed(123) で生成したたった1回のシミュレーション結果にすぎません。

たまたま今回生成されたデータでdependentモデルがうまくいっただけかもしれません。

さらに、推定値が滑らかになったからといって、それだけで良いモデルとは限りません。強い事前分布を設定すれば推定値を滑らかにすること自体はできますが、重要なのは見た目がきれいになることではなく、本当に正解に近づいていることです。

幸い、今回はシミュレーションなので「正解」の beta を私たち自身が知っています。

だったら、1回だけではなく何度も世界を作り直して、本当にdependentモデルのほうが正解に近づきやすいのかを確認すればよいでしょう。

ということで次は、データ生成からモデル推定までを1,000回繰り返して検証してみます。

大規模シミュレーション検証

ここで、たまたまうまくいったシード値だけをチェリーピックした結果ではなく、時点間の依存関係を考慮したdependentモデルを、各時点を独立に扱うindependentモデルよりも本当に優先して使ってよいのかをより体系的に判断するため、シミュレーションを1,000回繰り返して性能を比較します。

ちなみに余談ですが、私は普段かなりベイズモデルを使っているので「ベイジアンなのかな？」と思われがちですが、違います。ベイズという都合のいい推定手法を酷使している頻度論者ですwww

今回のように、データ生成そのものを何度も繰り返し、それぞれのデータに対する推定量のパフォーマンスを比較するという発想は、まさに頻度論的な考え方です。

「この1回のデータでうまくいったか」ではなく、同じデータ生成過程から何度もデータが得られたとしたら、どちらの推定方法が安定して良い結果を出すのかを見ていきましょう。

ただ、1,000回のシミュレーションをやると、一回一回の結果を図で可視化するのが難しくなるため、二乗平均平方根誤差（RMSE）を記録して全体的な比較を実施します：

future::plan(future::multisession(workers = 8))
seed_simulation_df <- 1000 |>
  seq_len() |>
  furrr::future_map(
    \(simulation_id){
      set.seed(simulation_id)

      coef_master <- 500 |>
        seq_len() |>
        tibble::tibble(
          time = _
        ) |>
        dplyr::mutate(
          beta_init = cumsum(rnorm(dplyr::n()))
        ) |>
        dplyr::mutate(
          beta = (beta_init - mean(beta_init))/sd(beta_init)
        ) |>
        dplyr::select(!beta_init)

      sim_data <- 500 |>
        seq_len() |>
        purrr::map(
          \(i){
            if (dplyr::between(i, 200, 250)){
              rep(i, 1)
            } else {
              rep(i, 500)
            }
          }
        ) |>
        unlist() |>
        tibble::tibble(
          time = _
        ) |>
        dplyr::left_join(
          coef_master, by = "time"
        ) |>
        dplyr::mutate(
          e = rnorm(dplyr::n()),
          y = beta + e
        )

      m_d_init <- cmdstanr::cmdstan_model("time_series_dependent.stan")

      m_d_estimate <- m_d_init$variational(
        seed = 1,
        data = list(
          time_type = nrow(coef_master),

          N = nrow(sim_data),
          time = sim_data$time,
          y = sim_data$y
        )
      )

      m_d_summary <- m_d_estimate$summary()

      m_id_init <- cmdstanr::cmdstan_model("time_series_independent.stan")

      m_id_estimate <- m_id_init$variational(
        seed = 1,
        data = list(
          time_type = nrow(coef_master),

          N = nrow(sim_data),
          time = sim_data$time,
          y = sim_data$y
        )
      )

      m_id_summary <- m_id_estimate$summary()

      tibble::tibble(
        id = simulation_id,
        d_rmse = m_d_summary |>
          dplyr::filter(stringr::str_detect(variable, "beta")) |>
          dplyr::bind_cols(answer = coef_master$beta) |>
          dplyr::mutate(e2 = (answer - mean)^2) |>
          dplyr::pull(e2) |>
          mean() |>
          sqrt(),
        id_rmse = m_id_summary |>
          dplyr::filter(stringr::str_detect(variable, "beta")) |>
          dplyr::bind_cols(answer = coef_master$beta) |>
          dplyr::mutate(e2 = (answer - mean)^2) |>
          dplyr::pull(e2) |>
          mean() |>
          sqrt()
      )
    },
    .progress = TRUE,
    .options = furrr::furrr_options(seed = 1)
  ) |>
  dplyr::bind_rows()

かなり時間がかかる処理でした！では、結果を可視化しましょう：

g_density <- seed_simulation_df |> 
  tidyr::pivot_longer(!id, names_to = "model", values_to = "rmse") |>
  ggplot2::ggplot() + 
  ggplot2::geom_density(ggplot2::aes(x = rmse, fill = model), alpha = 0.3)

g_scatter <- seed_simulation_df |>
  dplyr::mutate(
    dependent_win = d_rmse < id_rmse
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = d_rmse, y = id_rmse, color = dependent_win), alpha = 0.6)

gridExtra::grid.arrange(g_density, g_scatter)

1,000回やってみると、正しいモデルでも100%は勝てない

結果を見ると、時点間の依存関係を考慮したdependentモデルのRMSE分布は、independentモデルよりも全体として小さい方向に寄っていることがわかります。

実際に、それぞれのシミュレーションでどちらのRMSEが小さかったかを確認すると、dependentモデルの勝率は約73%となりました。

> seed_simulation_df |>
   dplyr::mutate(
     dependent_win = d_rmse < id_rmse
   ) |>
   dplyr::pull(dependent_win) |>
   mean()
[1] 0.73

つまり、時系列の依存関係をモデルに組み込むことで、平均的には推定精度を改善できるものの、毎回必ず勝てるわけではないという結果になりました。

「正しいモデル」なのに100%勝てない

ここで個人的にかなり面白いと思うポイントがあります。

今回のシミュレーションでは、真の beta を以下のようなランダムウォークとして生成しています。

beta_init = cumsum(rnorm(dplyr::n()))

一方、dependentモデルでは、

beta[i] ~ normal(beta[i - 1], rho[2]);

としているため、少なくとも beta の時間的な変化については、データ生成過程と整合的な構造をモデル側にも与えています。

それでも勝率は100%ではなく、約73%です。

これは重要なポイントだと思います。

データ生成過程をより正しく表現したモデルを使ったからといって、有限のデータから得られる一回一回の推定結果で、必ず他のモデルより高い精度を出せるわけではありません。

毎回生成される観測値にはランダムなノイズがあります。たまたまindependentモデルにとって推定しやすいデータが生成されることもあれば、dependentモデルが時系列方向に情報を借りることで、逆にその回の正解から少し遠ざかることもあります。

つまり、

正しいモデルを使う＝毎回必ず勝つ

ではなく、

正しい構造を利用することで、データ生成を何度も繰り返したときの性能を改善できる

と考えるほうが適切です。

まさにこの違いを見るために、1つのシード値だけではなく1,000回シミュレーションする意味があります。

dependentモデルは「大外し」を抑えているように見える

RMSEの分布を見ると、もう一つ興味深い特徴があります。

今回の1,000回のシミュレーションでは、independentモデルで観測された最大RMSEのほうが、dependentモデルで観測された最大RMSEよりも大きくなっています。

逆に分布の小さい側を見ると、dependentモデルでは、independentモデルではほとんど到達していないような低いRMSEも観測されています。

もちろん、これは今回の1,000回のシミュレーションで観測された最大値・最小値についての話であり、理論的にdependentモデルのRMSEに上限があるという意味ではありません。

しかし少なくとも今回の実験では、

dependentモデルは非常に良い推定結果を出せるケースを増やしつつ、independentモデルで発生するような大きな推定誤差もある程度抑えている

ように見えます。

これは今回のビジネス上のシナリオを考えると、平均RMSEの改善とは別の意味でも重要です。

データが突然500件から1件になったとき、「普段より多少精度が落ちる」こと以上に怖いのは、たまたま取得できた1件に引っ張られて、とんでもなく外れた数字を出してしまうことだからです。

時系列方向への情報の借用は、そのような極端な推定を抑える保険としても機能している可能性があります。

2つのモデルのRMSEは、意外と一緒には動かない

散布図も面白い結果になっています。

横軸がdependentモデルのRMSE、縦軸がindependentモデルのRMSEですが、両者の間には、少なくとも図を見る限り非常に強い相関があるようには見えません。

もし、

「難しいデータが生成された回では、どちらのモデルも同じように悪くなる」

だけであれば、散布図にはもっと明確な右上がりの関係が現れるはずです。

しかし今回の結果では、dependentモデルのRMSEが比較的小さいのにindependentモデルでは大きなRMSEになっているケースもあれば、その逆も存在します。

つまり、シミュレーションごとの「難しさ」だけで両モデルの性能差が決まっているわけではなく、生成されたデータに対して、それぞれのモデルが持つ情報の借用の仕組みが異なる形で作用していることが示唆されます。

相関係数を計算すると：

> cor.test(seed_simulation_df$d_rmse, seed_simulation_df$id_rmse, method = "pearson")

	Pearson's product-moment correlation

data:  seed_simulation_df$d_rmse and seed_simulation_df$id_rmse
t = 6.7504, df = 998, p-value = 2.498e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1488976 0.2674898
sample estimates:
      cor 
0.2089618

相関係数は0.2くらいで、統計的に有意ですが、値としては低いと言えるでしょう。

だから「1回当たった」ではなく、繰り返したときの性能を見る

以上をまとめると、今回のシミュレーションから得られるメッセージは、

時系列構造を正しくモデル化すれば必ず勝てる

というものではありません。

むしろ、

時系列構造が実際に存在する状況では、その構造を利用して情報を借用することで、確率的にいうと推定性能を改善できる可能性が高い。ただし、個々のデータセットではシンプルなモデルが勝つこともある。

という結果です。

個人的には、dependentモデルが100%勝たなかったことも含めて、かなり現実的な結果だと思います。

モデルを正しく設計する目的は、未来のあらゆるデータセットで絶対に勝つことではありません。

同じような問題に何度も直面したときに、より高い確率で、より小さな誤差を出せる意思決定を選ぶことにあります。

今回の結果では、その選択がdependentモデルだった、ということになります。

結論

いかがでしたか？

今回は、「情報プーリングって結局、何が嬉しいの？」という疑問から出発して、普段500件取得できていたデータが、ある期間だけ突然1件になってしまうという極端なシナリオを使って検証してみました。

比較したのは、どちらも階層ベイズモデルです。

各時点のパラメータを独立に扱うモデル
前の時点との依存関係を事前分布に組み込み、時系列方向に情報をプーリングするモデル

データが十分に存在する期間では、どちらのモデルも正解にかなり近い値を推定できました。

一方、データが500件から1件まで減少すると、各時点を独立に扱うモデルの推定値は大きく不安定になりました。それに対して、時系列の依存関係を組み込んだモデルでは、前後の時点から情報を借りることで、データがほとんどない期間でも比較的安定した推定を行うことができました。

さらに、「たまたま都合のよいシード値だっただけでは？」という可能性を排除するため、データ生成からモデル推定までを1,000回繰り返しました。

その結果、時系列の依存関係を考慮したモデルは、73%のシミュレーションで独立モデルより小さなRMSEを記録しました。

興味深いのは、今回のdependentモデルがデータ生成過程と整合的な構造を持っているにもかかわらず、勝率が100%ではなかったことです。

正しいモデルを使えば、どんなデータでも必ず勝てるわけではありません。

有限のデータには偶然性があります。ある1回のデータセットでは、よりシンプルなモデルのほうが正解に近づくこともあります。

それでも、同じデータ生成過程から何度もデータを取得すると考えたとき、適切な構造をモデルに組み込むことで、より高い確率で良い推定結果を得られる。今回の1,000回のシミュレーションは、そのことを非常にわかりやすく示してくれたと思います。

情報の借用は「データを増やす魔法」ではない

もちろん、今回の結果は、

データが足りなくなったら、とりあえずベイズを使えば何とかなる！

という話ではありません。

情報の借用によって、存在しないデータを新しく作り出せるわけではありません。借用しているだけです。

今回のdependentモデルが情報を借りることができたのは、「隣り合う時点の値は似ている」という構造が実際のデータ生成過程に存在していたからです。

もし現実には時点間の関係がほとんどなかったり、途中で大きな構造変化が発生していたりすれば、過去の情報を積極的に借りることが、かえって推定を悪化させる可能性もあります。

また、本記事の冒頭でも触れたように、「どのデータが欠けるか」が推定対象そのものと関係している場合には、今回とは別の問題になります。

重要なのは、何でもプーリングすることではなく、現実に存在すると考えられる構造をモデルに組み込み、そこから情報を借りることです。

「事前分布を設計する」の意味

ベイズ統計を勉強していると、

事前分布を適切に設計しましょう
情報を借用できます
シュリンケージが働きます

といった説明をよく目にします。

私自身、最初にこれらの説明を読んだときは、

で、それの何が嬉しいの？

と思っていました。

しかし今回のように、昨日まで500件あったデータが突然1件になった世界を考えると、その意味はかなり具体的になります。

手元の1件だけで戦う必要はありません。

「昨日と今日はある程度似ている」「同じ店舗同士は似ている」「同じ地域の顧客は似ている」といった、データの背後にある構造そのものを利用して、情報が少ない場所に情報が豊富な場所の知識を分けてもらう。

それが、階層ベイズにおける情報プーリングの大きな魅力の一つだと思います。

そして、もし明日の会議までに数字を出さなければならないのに、担当者がまだ戻ってこなかったら、担当者の帰りを待ちながら、ベイズモデルを推定してみるのも一つの手かもしれませんw

最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：