はじめに
こんにちは、事業会社で働いているデータサイエンティストです。
本記事は、前回の記事の続編です。
前回の記事では、指数二乗カーネルを用いたガウス過程が、ランダムウォークで生成されたデータのデータ不足区間を推定するタスクにおいて、より単純なモデルにボロ負けするという結果を、シミュレーションを通じてお見せしました。
この結果については、
「いや、指数二乗カーネルは滑らかな関数を想定しているんだから、ギザギザしたランダムウォークを突っ込んだら壊れるのは当たり前では？」
と思った方もいるでしょう。
はい。その通りです。
前回の記事を取り上げていただいたYouTube動画でも、まさにこの点を強調していただきました。
では、なぜ私はそんな「相性の悪い組み合わせ」をあえて試したのでしょうか。
理由は単純です。
現実のデータでは、モデルを推定する前に「真の潜在関数」がどんな形をしているのか分からないからです。
たとえば、求人の採用数を予測するときに、「賃金と採用数の間にどのような非線形な関係があるのか」をガウス過程で表現したいとします。
あるいは、不動産価格を予測するときに、緯度・経度から潜在的な地理効果を推定したいとします。
こうした問題で、モデルを推定する前に真の潜在関数を眺めて、
「あ、こいつギザギザしてるから指数二乗カーネルは無理ンゴ」
と判断することはできません。
そもそも、その関数が見えないから推定したいのです。
もちろん、ドメイン知識を使ったり、複数のカーネルを比較したり、事前予測チェックを行ったりすることで、より妥当なモデルを選ぶことはできます。
しかし、指数二乗カーネルはガウス過程で非常によく使われる代表的なカーネルの一つです。
だから前回の記事では、あえてこの定番カーネルを「モデルにとって都合の悪い世界」に放り込み、どこまで耐えられるのかをストレステストしました。
そして今回は、そこからさらに一歩進めます。
前回は、
ランダムウォーク vs. 指数二乗カーネル
という、かなり極端な対決でした。
そこで今回は、ランダムウォークを少しずつ、少しずつ、どんどん滑らかにしていきます。
そして、
真の潜在関数がどの程度滑らかになれば、指数二乗カーネルを使ったガウス過程は安定して高い精度を出せるようになるのか？
を、大量のシミュレーションで確認します。
言い換えれば今回は、ガウス過程にいきなり苦手な問題を投げつけて「爆死した！」と笑う記事ではありません。
ガウス過程が苦手な世界から得意な世界へ、データ生成過程を少しずつ変化させながら、その性能がどう変わっていくのかを観察する実験です。
結果を先にお見せすると、かなり面白いことになりました。
潜在関数を滑らかにしていくにつれて、ガウス過程の典型的な予測誤差は確かに小さくなっていきます。
ここまでは、ある意味で期待通りです。
ところが――
滑らかにすれば、いつでも安定して高精度になるわけではありません。
大量にシミュレーションを繰り返して誤差の「平均」だけではなく分布そのものを見てみると、典型的にはかなり高精度になっている条件でも、一部のシードでは予測誤差が大きく跳ね上がります。
つまり、
「平均的には強い」と「安心して任せられる」は、同じ意味ではない
ということです。
もちろん、
「ハイパーパラメータをもっとちゃんとチューニングすればいいのでは？」
「別のカーネルも比較すべきでは？」
「推論方法も最適化すべきでは？」
という指摘は当然あり得ます。
そして、それらはすべて正しいです。
ただし、今回はあえてやりません。
今回確認したいのは「ガウス過程を限界までチューニングすればどこまで強くなるか」ではなく、
代表的な指数二乗カーネルを使ったガウス過程に、深く考えずにデータを渡したら、どの程度安心して任せられるのか？
という問題だからです。
実務では残念ながら、LightGBM、ガウス過程、深層学習などの柔軟なモデルにデータを渡せば、モデル側がいい感じに何とかしてくれるだろう、という「ディープでポン」的なスタンスがゼロではないわけなので、モデルが柔軟だから、何も考えずに使ってもロバストであるなのかを丁寧に検証するのが非常に大事になります。
モデルの名前が強そうだからといって、データ生成過程までモデルに忖度してくれるわけではありません。
というわけで今回は、ガウス過程くんに大量のデータ生成過程をぶつけます。
どこまで滑らかにすれば強くなるのか。
強くなったように見えても、どの程度「事故」が残るのか。
そして、「柔軟なモデルに任せておけば何とかなる」という期待は、どこまで正しいのか。
大量のシミュレーションで確認していきましょう。
ガウス過程5000回酷使劇場、開幕です。
劇場セッティング
今回の記事は前回の記事の数学的な説明とStanファイルを前提としています。数式やStan言語のモデルコードを確認したい方は前回の記事をご覧ください：
では、今回のシミュレーションで何をするのかを説明します。
前回の記事では、真の潜在的な平均値としてランダムウォークを使いました。
set.seed(1)
500 |>
rnorm() |>
cumsum() |>
tibble::tibble(
x = _
) |>
dplyr::mutate(
time = dplyr::row_number()
) |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = x))
株化の動きみたいですね。ランダムウォークは、前の時点から少しずつ値が変化していくため時点間の依存関係は強い一方で、指数二乗カーネルが好むような非常に滑らかな関数ではありません。
そこで今回は、このランダムウォークを出発点として、少しずつ滑らかな関数へ変化させてみます。
まず、何をしているのかを簡単に可視化してみましょう。
set.seed(1)
500 |>
rnorm() |>
cumsum() |>
tibble::tibble(
x = _
) |>
dplyr::mutate(
smooth_20_x = slider::slide_dbl(x, mean, .before = 20, .after = 20)
) |>
dplyr::mutate(
time = dplyr::row_number()
) |>
tidyr::pivot_longer(!time) |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = value, color = name))
ここで行っている操作は非常に単純です。
まず、
rnorm(500) |> cumsum()
によって、500時点のランダムウォークを生成しています。
そして、
slider::slide_dbl(x, mean, .before = 20, .after = 20)
によって、それぞれの時点について前後20時点の値を使った移動平均を計算しています。
詳細はこちらの記事をご参考ください：
要するに、中央付近の時点であれば、
「その時点の値をそのまま使うのではなく、前20時点＋現在＋後20時点、最大41個の値の平均に置き換える」
という操作です。
当然ながら、こうするとランダムウォークに存在していた細かなギザギザが消えていきます。
元のランダムウォークと比較すると、大きな上昇・下降といった長期的な形状はある程度残っている一方で、短期的な細かい変動が平均化され、かなり滑らかな曲線になっていることが分かると思います。
ここが今回の実験のポイントです。
前回は、
ギザギザしたランダムウォークを指数二乗カーネルのガウス過程に突っ込む
という、ガウス過程にとってかなり厳しい条件でした。
今回はそこから、移動平均に使う範囲を徐々に大きくしていきます。
たとえば、
- 平滑化なし：元のランダムウォーク
- 前後10時点で平滑化
- 前後20時点で平滑化
- 前後30時点で平滑化
- …
- 前後90時点で平滑化
というように、同じランダムウォークから滑らかさの異なる潜在関数を作ります。
重要なのは、完全に別々のデータ生成過程を適当に作って比較するのではなく、同じランダムウォークを出発点として、平滑化の強さだけを段階的に変えていることです。
これによって、
潜在関数を滑らかにして、指数二乗カーネルが想定する世界に少しずつ近づけていったとき、ガウス過程の推定精度はどう変化するのか？
を観察できます。
もちろん、slider による移動平均の窓幅そのものを、ガウス過程における数学的な「滑らかさ」の尺度だと解釈することはできません。
今回操作しているのは、あくまでデータ生成過程側の平滑化の強さです。
したがって、本記事で知りたいのは、
「smooth = 30ならガウス過程にとってこの程度の滑らかさである」
といった普遍的な対応関係ではありません。
むしろ、
同じ系列を人為的にどんどん滑らかにしていったとき、指数二乗カーネルを使ったガウス過程の誤差分布がどのように変化していくのか
というストレステストです。
もし前回の惨敗が本当に「指数二乗カーネルに対して潜在関数がギザギザすぎた」ことによって起きているのであれば、平滑化を強くするにつれてRMSEは下がっていくはずです。
逆に、十分滑らかにしても大きな誤差が頻繁に発生するのであれば、
「指数二乗カーネルが得意そうな世界に近づければ、それだけで安心して任せられる」
というほど単純な話でもなさそうです。
では、ランダムウォークをどんどん滑らかにしながら、ガウス過程くんを大量投入してみましょう。
一回のシミュレーション
では、大量のシミュレーションを回す前に、まずは一つの乱数シードを使って、
本当に潜在関数を滑らかにしていけば、ガウス過程の推定精度も改善していきそうなのか？
を目で確認してみましょう。
いきなり500シード × 10種類を回して集計結果だけを見るよりも、まずはガウス過程が実際にどのような関数を推定しているのかを確認した方が、後ほど出てくる結果も理解しやすくなります。
滑らかさの異なる10本の潜在関数を作る
まずはランダムウォークを1本生成し、それを段階的に平滑化したデータを作ります：
set.seed(111)
coef_master <- 200 |>
rnorm() |>
cumsum() |>
tibble::tibble(
beta_0 = _
) |>
purrr::reduce(
seq(10,90,length.out = 9),
\(df, smooth_id) {
df |>
dplyr::mutate(
!!stringr::str_c("beta_", smooth_id) :=
slider::slide_dbl(
beta_0,
mean,
.before = smooth_id,
.after = smooth_id,
.complete = FALSE
)
)
},
.init = _
) |>
dplyr::mutate(
dplyr::across(
dplyr::starts_with("beta_"),
\(x) (x - mean(x)) / sd(x)
)
)
ここでは、まずbeta_0として元のランダムウォークを作り、そこから移動平均の範囲を10、20、30、……、90と広げながら、9種類の平滑化された系列を作っています。
最後に、それぞれの系列を平均0、標準偏差1に標準化しています。
これは重要な処理です。平滑化を強くすると系列の分散そのものも小さくなってしまうため、そのまま比較すると「滑らかになったから推定しやすくなった」のか、「単に変動幅が小さくなったからRMSEが小さくなった」のかが分かりにくくなります。
そこで、すべての系列のスケールを揃えた上で、主に関数形の違いを比較できるようにしています。
こんな感じのデータが出来上がります：
> coef_master
# A tibble: 200 × 10
beta_0 beta_10 beta_20 beta_30 beta_40 beta_50 beta_60 beta_70 beta_80 beta_90
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0.397 -0.387 -0.752 -0.905 -1.32 -1.85 -1.98 -1.94 -1.82 -1.77
2 0.330 -0.459 -0.770 -0.943 -1.38 -1.85 -1.96 -1.90 -1.78 -1.74
3 0.267 -0.488 -0.786 -0.967 -1.44 -1.84 -1.93 -1.86 -1.75 -1.71
4 -0.201 -0.507 -0.801 -1.00 -1.50 -1.83 -1.91 -1.82 -1.72 -1.68
5 -0.236 -0.511 -0.812 -1.05 -1.55 -1.83 -1.87 -1.79 -1.70 -1.65
6 -0.208 -0.536 -0.816 -1.09 -1.59 -1.82 -1.84 -1.75 -1.67 -1.63
7 -0.512 -0.560 -0.838 -1.13 -1.63 -1.81 -1.81 -1.71 -1.64 -1.60
8 -0.717 -0.586 -0.833 -1.16 -1.67 -1.80 -1.78 -1.68 -1.62 -1.57
9 -0.910 -0.623 -0.827 -1.20 -1.69 -1.78 -1.75 -1.64 -1.60 -1.55
10 -1.01 -0.653 -0.832 -1.24 -1.70 -1.77 -1.72 -1.60 -1.58 -1.53
# ℹ 190 more rows
# ℹ Use `print(n = ...)` to see more rows
つまり、beta_0からbeta_90に向かうにつれて、元は同じランダムウォークだった系列がどんどん滑らかになっていきます。
10本すべてにガウス過程を当てる
では次に、それぞれの系列を「真の潜在平均」とみなし、各時点につき100件の観測値を生成します。
そして、10種類のデータそれぞれに対して、前回の記事と同じ指数二乗カーネルを使ったガウス過程を推定します。
future::plan(future::multisession(workers = 10))
furrr::future_map(
1:ncol(coef_master),
\(this_smooth){
set.seed(123)
sim_data <- tibble::tibble(
time = rep(
seq_len(nrow(coef_master)),
each = 100
),
y = rnorm(
nrow(coef_master) * 100,
mean = rep(coef_master[[this_smooth]], each = 100)
)
)
m_gp_init <- cmdstanr::cmdstan_model("time_series_gp.stan")
m_gp_estimate <- m_gp_init$variational(
seed = 1,
data = list(
time_type = nrow(coef_master),
time_seq = (1:nrow(coef_master) - mean(1:nrow(coef_master)))/sd(1:nrow(coef_master)),
N = nrow(sim_data),
time = sim_data$time,
y = sim_data$y
)
)
m_gp_summary <- m_gp_estimate$summary()
m_gp_summary |>
dplyr::filter(stringr::str_detect(variable, "^f_trend\\[")) |>
dplyr::bind_cols(y = coef_master[[this_smooth]]) |>
dplyr::mutate(time = dplyr::row_number()) |>
ggplot2::ggplot() +
ggplot2::geom_point(ggplot2::aes(x = time, y = y)) +
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::labs(
title = stringr::str_c("smooth level: ", this_smooth)
)
},
.progress = TRUE,
.options = furrr::furrr_options(seed = TRUE)
) |>
gridExtra::grid.arrange(
grobs = _,
ncol = 5,
nrow = 2
)
結果はこちらです。
おおおおおおおお！！！！
かなり分かりやすい結果になりました。
黒い点が真の潜在平均、青い線がガウス過程による推定値、青い帯が推定された信用区間です。
まずsmooth level: 1、つまり平滑化していない元のランダムウォークを見てください。
真の系列は細かく上下しているのに対して、ガウス過程はその変動を追い切れず、かなり滑らかな曲線を引いています。前回の記事で見た現象とほぼ同じです。
ところが、平滑化を強くしていくと様子が変わります。
smooth level: 2、3あたりではまだ局所的なズレが目立ちますが、4、5、6と進むにつれて、青い推定線が黒い真値にかなり近づいていきます。
そして後半になると、
おお！お前、急に仕事できるようになったな 😎
と言いたくなるくらい、ガウス過程が真の潜在関数をきれいに追跡しています。
これはまさに、今回確認したかった現象です。
指数二乗カーネルを使ったガウス過程が常に「ランダムウォークを推定できない」のではありません。
データ生成過程を徐々に滑らかにし、指数二乗カーネルが得意とする世界へ近づけていくと、推定結果も明らかに改善しているように見えます。
少なくともこのシードだけを見る限り、
「前回ガウス過程がボロ負けしたのは、ガウス過程そのものが弱いからではなく、指数二乗カーネルが想定する滑らかさとデータ生成過程が合っていなかったからではないか」
という説明には、かなり説得力がありそうです。
しかし、、、
ここで勝利宣言するわけにはいきません。
前回の記事で散々確認したように、一つの乱数シードだけを見てモデルの性能を語るのは危険です。
今回たまたま、非常にきれいな結果が出るランダムウォークを引いただけかもしれません。
別のランダムウォークを生成したら、十分に平滑化したにもかかわらずガウス過程が盛大に外すかもしれません。
というわけで、一回のシミュレーションで喜ぶのはここまでです。
次は乱数シードを変えながら同じ実験を繰り返します。
500種類のランダムウォーク × 10段階の滑らかさ。
つまり、、、
5000本のガウス過程をリングに投入します。
ここからが、本当の酷使劇場です。
大規模シミュレーション劇場
では、先ほどのセッティングで乱数シードを500回変えながら、大規模にシミュレーションを回してみましょう。
各シードでは、まず200時点のランダムウォークを1本生成します。そこから移動平均の範囲を段階的に広げることで、元のランダムウォークを含む滑らかさの異なる10本の潜在関数を作ります。
つまり、
- 500種類のランダムウォーク
- 各ランダムウォークについて10段階の滑らかさ
- 合計 500 × 10 = 5000回のガウス過程推定
を行うことになります。
5000本です。
さすがにガウス過程くんも「ごっくーん（筆者のあだ名）！そんなんの聞いてないよ！」と言いたくなる仕事量ですが、タイトル通り酷使していきましょう。
コードはこちらです：
future::plan(future::multisession(workers = 20))
m_gp_init <- cmdstanr::cmdstan_model("time_series_gp.stan")
gaussian_process_simulation_df <- 500 |>
seq_len() |>
furrr::future_map(
\(this_seed){
set.seed(this_seed)
coef_master <- 200 |>
rnorm() |>
cumsum() |>
tibble::tibble(
beta_0 = _
) |>
purrr::reduce(
seq(10,90,length.out = 9),
\(df, smooth_id) {
df |>
dplyr::mutate(
!!stringr::str_c("beta_", smooth_id) :=
slider::slide_dbl(
beta_0,
mean,
.before = smooth_id,
.after = smooth_id,
.complete = FALSE
)
)
},
.init = _
) |>
dplyr::mutate(
dplyr::across(
dplyr::starts_with("beta_"),
\(x) (x - mean(x)) / sd(x)
)
)
purrr::map(
1:ncol(coef_master),
\(this_smooth){
set.seed(123)
sim_data <- tibble::tibble(
time = rep(
seq_len(nrow(coef_master)),
each = 100
),
y = rnorm(
nrow(coef_master) * 100,
mean = rep(coef_master[[this_smooth]], each = 100)
)
)
gp_result <- tryCatch(
{
m_gp_estimate <- m_gp_init$variational(
seed = 1,
data = list(
time_type = nrow(coef_master),
time_seq = (1:nrow(coef_master) - mean(1:nrow(coef_master)))/sd(1:nrow(coef_master)),
N = nrow(sim_data),
time = sim_data$time,
y = sim_data$y
)
)
m_gp_summary <- m_gp_estimate$summary()
list(
success = TRUE,
summary = m_gp_summary
)
},
error = function(e){
message(
"GP failed",
conditionMessage(e)
)
list(
success = FALSE,
summary = NULL
)
}
)
if (gp_result$success == TRUE){
gp_result$summary |>
dplyr::filter(stringr::str_detect(variable, "^f_trend\\[")) |>
dplyr::bind_cols(y = coef_master[[this_smooth]], seed = this_seed, smooth = this_smooth) |>
dplyr::mutate(time = dplyr::row_number()) |>
dplyr::select(seed, smooth, time, y, mean, q5, q95)
} else {
tibble::tibble(
y = coef_master[[this_smooth]]
) |>
dplyr::bind_cols(seed = this_seed, smooth = this_smooth) |>
dplyr::mutate(
time = dplyr::row_number(),
mean = NA_real_,
q5 = NA_real_,
q95 = NA_real_
) |>
dplyr::select(seed, smooth, time, y, mean, q5, q95)
}
}
) |>
dplyr::bind_rows()
},
.progress = TRUE,
.options = furrr::furrr_options(seed = TRUE)
) |>
dplyr::bind_rows()
まずは推定失敗率を確認する
RMSEを見る前に、まず確認しておきたいのがそもそもガウス過程の推定が最後まで成功したのかという点です。
今回のシミュレーションでは、推定に失敗したケースについてmeanをNAとして保存しているため、各seed × smoothについてすべてのmeanがNAになっているかどうかを確認すれば、推定失敗を判定できます。
> gaussian_process_simulation_df |>
dplyr::summarise(fail = all(is.na(mean)), .by = c(seed, smooth)) |>
dplyr::summarise(fail_rate = mean(fail), .by = smooth)
# A tibble: 10 × 2
smooth fail_rate
<int> <dbl>
1 1 0.08
2 2 0.08
3 3 0.084
4 4 0.076
5 5 0.084
6 6 0.108
7 7 0.062
8 8 0.046
9 9 0.034
10 10 0.038
これ、かなり面白い結果です。
まず大きな傾向としては、潜在関数を強く平滑化していくと、後半では推定失敗率が低下していることが分かります。
smooth = 1では8.0%、smooth = 6では10.8%もの推定が失敗していますが、smooth = 8では4.6%、smooth = 9では3.4%、smooth = 10では3.8%まで下がっています。
つまり、指数二乗カーネルが想定するような滑らかな世界に近づいていくと、精度だけではなく、推定そのものの安定性も改善している可能性があります。
ただし、ここで重要なのは、
滑らかにすればするほど、失敗率が一直線に下がっていくわけではない
ということです。
実際、
-
smooth = 1：8.0% -
smooth = 2：8.0% -
smooth = 3：8.4% -
smooth = 4：7.6% -
smooth = 5：8.4% -
smooth = 6：10.8%
となっており、前半から中盤にかけてはほとんど改善していません。
むしろsmooth = 6では、今回の10段階の中で**最も高い10.8%**まで悪化しています。
この点はかなり重要です。
もし単純に、
「潜在関数が滑らかになれば、指数二乗カーネルとの相性が良くなって、推定失敗も順調に減っていくだろう」
という話であれば、もっと素直な右肩下がりを期待したくなります。
しかし実際にはそうなっていません。
明確に5%を下回るのはsmooth = 8以降です。つまり、かなり強く平滑化して、ようやく推定失敗率が5%未満まで下がるという結果になっています。
言い換えると、少しくらい潜在関数を滑らかにしただけでは、
「これでもう安心！」
とは到底言えません。
指数二乗カーネルがかなり得意そうな世界に近づけても、ある程度の平滑化までは数%〜10%程度の推定失敗が普通に残っています。
可視化すると、さらに分かりやすいです：
gaussian_process_simulation_df |>
dplyr::summarise(fail = all(is.na(mean)), .by = c(seed, smooth)) |>
dplyr::summarise(fail_rate = mean(fail), .by = smooth) |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = smooth, y = fail_rate))
このグラフからも、失敗率が単調に減少しているわけではなく、途中で上下しながら、かなり後半になってようやく明確に低い水準へ落ちていくことが分かります。
もちろん、この推定失敗がすべて「指数二乗カーネルのモデルとしての弱さ」によって起きていると断定することはできません。
今回使っているのはStanの変分推論なので、最適化の不安定性や共分散行列の数値計算なども関係している可能性があります。
したがって、ここから言えるのは、
今回の「指数二乗カーネル × この事前分布 × 変分推論」という実装では、潜在関数をかなり滑らかにするまで、数%規模の推定失敗が残り続けた
ということです。
そして実務的には、これは無視しづらいポイントです。
モデルの平均的な精度がどれだけ高くても、100回動かしたら5回や10回推定に失敗するのであれば、そのまま安定運用できるとは言いづらいでしょう。
一方で、後半では失敗率が3〜4%程度まで下がっていることから、潜在関数とカーネルの相性が数値的な安定性にも関係している可能性は十分にありそうです。
ただし、まだこれは「推定できたか、できなかったか」の話にすぎません。
次に見たいのは当然、
推定に成功したガウス過程は、本当に滑らかになるほど高精度になっているのか？
という問題です。
ここから、5000本のガウス過程のRMSE分布を見ていきましょう。
滑らかになれば、本当にガウス過程は強くなるのか？
では、いよいよ今回の本題です。
先ほどは「そもそも推定に成功するのか」を確認しましたが、ここからは推定に成功したガウス過程が、真の潜在関数をどの程度正確に復元できているのかを見ていきます。
各シード・各滑らかさについてRMSEを計算し、その分布を可視化します。
gaussian_process_simulation_df |>
tidyr::drop_na() |>
dplyr::summarise(
rmse = sqrt(mean((y - mean)^2, na.rm = TRUE)),
.by = c(seed, smooth)
) |>
dplyr::summarise(
median = median(rmse),
q5 = quantile(rmse, 0.05),
q95 = quantile(rmse, 0.95),
q65 = quantile(rmse, 0.65),
q35 = quantile(rmse, 0.35),
q75 = quantile(rmse, 0.75),
q25 = quantile(rmse, 0.25),
.by = smooth
) |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = smooth, y = median)) +
ggplot2::geom_ribbon(ggplot2::aes(x = smooth, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) +
ggplot2::geom_ribbon(ggplot2::aes(x = smooth, ymin = q35, ymax = q65), fill = ggplot2::alpha("blue", 0.3)) +
ggplot2::geom_ribbon(ggplot2::aes(x = smooth, ymin = q25, ymax = q75), fill = ggplot2::alpha("blue", 0.3))
この図は、各 smooth における500回のシミュレーションのRMSE分布をまとめたものです。
中央の線はRMSEの中央値を表しています。そして、その周囲にある3層の帯は、それぞれRMSE分布の異なる範囲を表しています。
- 最も内側の帯：35%点〜65%点
- 中間の帯：25%点〜75%点
- 最も外側の帯：5%点〜95%点
したがって、中央の線が低いほど典型的なシミュレーションでは精度が高いと解釈できます。一方、帯が広いほど、シードによってRMSEが大きく変動しており、推定精度の安定性が低いことを意味します。
つまり、この図では単に「平均的にどれくらい当たるか」だけではなく、500回繰り返して運用したとき、精度がどの程度安定しているのか、そして悪いケースではどこまで悪化しうるのかを同時に見ています。
この結果から、まず非常に分かりやすい傾向が確認できます。
はい。潜在関数が滑らかになるほど、ガウス過程の推定精度は全体として改善しています。
これはある意味では期待通りです。
今回使用している指数二乗カーネルは非常に滑らかな潜在関数を想定しています。そのため、ギザギザしたランダムウォークそのものから、移動平均によってどんどん滑らかな関数へ変化させていけば、データ生成過程がガウス過程の想定する世界に近づいていきます。
実際、RMSEの中央値を見ると、滑らかさが高まるにつれて明確に低下していきます。
前回の記事で見た「ガウス過程がランダムウォークをうまく捉えられない」という結果は、ガウス過程そのものがポンコツだったからではありません。
モデルが想定する世界と、実際に与えられた世界が合っていなかった。
そして今回、その世界を徐々にガウス過程好みに変えていくことで、ちゃんと性能が改善していくことが確認できました。
めでたしめでたし。
、、、で終われば非常にきれいなのですが、この図にはもう一つ、私がかなり気になった特徴があります。
中央値は優秀。でも、たまに盛大に外す
中央値の線だけを見れば、かなり安心感があります。
ところが、RMSEの分布全体を見ると話が少し変わってきます。
今回の図では、中央値だけではなく、25〜75%、35〜65%、さらに5〜95%のRSME分布の範囲まで表示しています。
中央付近のRMSEは滑らかさとともにかなり小さくなっています。つまり、典型的なケースでは、確かにガウス過程は非常によく働いています。
しかし、外側の5〜95%区間を見ると、かなり滑らかな潜在関数を与えた後でも、上側には無視できない長い裾が残っています。
言い換えると、
普段はかなり当てる。でも、一定の割合で突然かなり大きく外す。
という挙動です。
これは平均的・典型的な精度だけを見ていたら、かなり見落としやすいポイントだと思います。
「RMSEの中央値が低いから、このモデルは安定して高精度だ」
とは必ずしも言えません。
500回も違う世界を作って同じモデルを走らせてみると、その中にはかなり高いRMSEを出してしまうケースが少数ながら確実に存在することが分かります。
しかも前節で確認したように、これは推定そのものに失敗したケースを除外した後の話です。
つまり今回観察しているのは、
推定に成功したうえで、それでも盛大に外しているケース
です。
「柔軟性」はタダではないのかもしれない
ここから先は、数学的な定理というより、今回の実験結果を見た私なりの実務的な解釈です。
ガウス過程の大きな魅力の一つは、関数形を最初から
$$
y = \beta_0 + \beta_1 x
$$
のように強く決め打ちしなくてもよいことです。LightGBMや深層学習などの柔軟なモデルにも共通して言えるメリットです。
「どんな形なのかよく分からないから、ガウス過程に柔軟に推定してもらおう」
という使い方ができます。
これはものすごく強力です。
しかし、少し意地悪な言い方をすれば、これは関数形をどう設計するのかというモデリング作業の一部を、ガウス過程にアウトソースしているとも考えられます。
人間が、
「この変数とはこういう関係があるはずだ」
と細かく構造を指定する代わりに、
「そこはよしなにやってください」
とモデルに任せているわけです。
もちろん、本当に何でも勝手にやってくれるわけではありません。カーネルや事前分布を通じて、こちら側は依然としてかなり重要な仮定を置いています。
それでも、パラメトリックなモデルと比較すれば、潜在関数の形についてモデル側に任せている部分は大きいでしょう。
そして今回の結果を見ると、その自由には代金があるのかもしれません。
典型的なケースでは非常に高精度。
潜在関数が滑らかになれば、中央値のRMSEも明確に改善する。
それでも、何百回も繰り返してみると、一定割合でかなり悪い推定結果が現れる。そして前節で見たように、そもそも推定自体が失敗するケースも残っています。
ここで私が気にしているのは、「平均的にどれだけ賢いか」だけではありません。
何度も運用したときに、どれだけ安定して同じ品質を出してくれるのか。
いわば、モデルの運用上の頻度論的性能です。
実務では、
「95%くらいのケースではめちゃくちゃ強いです。でも残りでは何が起きるか分かりません」
では困る場面があります。
100回、1000回とモデルを動かすのであれば、平均的なRMSEだけではなく、悪いケースがどの程度の頻度で発生し、そのときどこまで壊れるのかも性能の一部です。
No Free Lunch
もちろん、このシミュレーションだけから、
「ノンパラメトリックモデルは不安定である」
などと一般化するつもりはありません。
今回の結果は、あくまでこのデータ生成過程、この指数二乗カーネル、この事前分布、そしてStanの変分推論という組み合わせについての結果です。
しかし、今回の5000回の推定を見ていると、少なくとも一つの実務的な教訓は感じます。
No Free Lunchです。
柔軟性は魔法ではありません。
関数形を人間が細かく設計する作業を減らし、その一部を柔軟なモデルに任せられることには、とてつもなく大きな価値があります。
しかし、その代わりに、「ほとんどの場合かなりうまくいく」ことと、「何度使っても安定してうまくいく」ことは別問題なのかもしれません。
モデルに自由を与えれば、その自由をうまく使って驚くほどきれいな推定をしてくれることがあります。
一方で、その自由を与えた以上、人間が構造を強く指定したモデルとは違う種類の不確実性や不安定性、もしくはその不確実性・不安定性を吸収する運用ルールの策定タスクを引き受ける必要があるかもしれません。
これはガウス過程への批判ではありません。
むしろ、柔軟なモデルを使うなら、その柔軟性から得られる利益だけではなく、その代金まで測ってから本番投入したいという話です。
5000本も走らせてみると、ガウス過程の表の顔――「滑らかな関数なら非常に高精度」――だけではなく、
「普段は優秀。でも、ごく一部でかなり派手に外す」
という裏の顔まで見えてきました。
そして、この裏の顔こそ、最先端モデルを教科書から実運用に連れていくときに確認しておきたいものではないでしょうか。
おわりに
今回、500種類のランダムウォークを10段階に平滑化し、合計5000回のガウス過程を推定してみました。
結果はある意味では期待通りでした。潜在関数を滑らかにして指数二乗カーネルが想定する世界に近づけるほど、典型的なRMSEは明確に改善しました。 一方で、推定失敗率は単調には低下せず、推定に成功したケースだけを見ても、少数ながらRMSEが大きく跳ねるケースが残りました。
ここで私が言いたいのは、
「だからガウス過程は使うな」
ではありません。
むしろ逆です。ガウス過程をはじめとする柔軟なノンパラメトリックモデルは、分からない関数形をモデル側に任せられる、非常に強力な道具です。
ただし、それはある意味で、本来人間が行うはずだったモデリング作業の一部をモデルにアウトソースするショートカットでもあります。そしてショートカットにはリスクがあります。今回の結果が示しているように、「典型的には高精度」と「何度動かしても安定して高精度」は同じではありません。
だから私は、可能なのであれば、ドメイン知識を学び、
なぜこの変数が効くのか？
どのような関係があるはずなのか？
どの構造なら現実の現象を表現できるのか？
を考え、その知識を統計モデルの構造として書いてしまう方が好きです。
そして、いろいろな構造を一つの巨大なモデルの中につなぎ合わせられることこそ、私がベイズモデリングを好きな最大の理由の一つです。
一方で、私は研究のためだけにモデルを作っているわけではありません。実際にオンラインで動くデータサイエンスシステムを持つデータサイエンスマネジャーとして、モデルが「理論的に美しい」「平均的に強い」だけでは困ります。
- 100回動かしたら何回壊れるのか。
- 悪いときにはどこまで悪くなるのか。
- 本番環境で繰り返し使っても、安心して任せられるのか。
だからこそ私は、ベイズモデルが大好きであると同時に、そのモデルを評価するときには非常に頻度論的です。
柔軟性はタダではありません。最先端のモデルも魔法ではありません。
モデルに任せる前に、自分たちが知っていることをモデルに教えられないか考える。そして任せた部分については、何百回でも酷使して、本当に安心して任せられるのかを確認する。
それが今回、5000本のガウス過程を酷使して私が改めて感じたことです。
最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：