階層ベイズモデルの基礎
2026-08-09
久保先生の 階層ベイズモデルの解説 の考え方を, 打率の推定でなぞります.
問題: 20人の打者について, シーズン序盤40〜50打数ぶんの成績があります. ここからシーズン終了時の打率を当てられるでしょうか.
library(data.table)
library(dplyr)
library(ggplot2)
library(rstanarm)
library(knitr)
options(mc.cores = min(4, parallel::detectCores()))
dat <- fread("kubo-data.csv")
dat |> head(5) |> kable(digits = 3)| FULLNAME | ATBAT | HIT | AVG | AVG_SEASON |
|---|---|---|---|---|
| 1 | 43 | 7 | 0.163 | 0.273 |
| 2 | 42 | 11 | 0.262 | 0.280 |
| 3 | 42 | 6 | 0.143 | 0.284 |
| 4 | 50 | 14 | 0.280 | 0.287 |
| 5 | 44 | 10 | 0.227 | 0.289 |
ATBATが打数, HITが安打数,
AVGが序盤の打率,
AVG_SEASONがシーズン終了時の打率です.
序盤の打率はあてにならない
dat |>
ggplot(aes(x = AVG, y = AVG_SEASON)) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
geom_point(size = 3) +
coord_fixed(xlim = c(0.10, 0.42), ylim = c(0.10, 0.42)) +
xlab("序盤の打率") + ylab("シーズン打率") +
ggtitle("点が破線に乗っていれば当たっている")
見ての通り, 縦方向にはほとんど広がっていないのに横方向には大きく散っています.
- 序盤の打率のばらつき: 0.071
- シーズン打率のばらつき: 0.014
観測されたばらつきの大半は実力差ではなく偶然です. 40打数で1本の安打は打率0.025ぶんに相当するので, これくらいは動きます.
3つの予測を比べる
rmse <- function(pred, truth) sqrt(mean((pred - truth)^2))
group_mean <- sum(dat$HIT) / sum(dat$ATBAT)- 最尤推定: 序盤の打率をそのまま使う.
- 全体平均: 誰であろうと全体の打率0.283と答える.
- 階層ベイズ: 1と2の間をデータに決めさせる.
3つ目のモデルはこうです.
\[ \text{HIT}_i \sim \mathrm{Binomial}(\text{ATBAT}_i,\ q_i), \qquad \mathrm{logit}(q_i) = a + b_i, \qquad b_i \sim \mathrm{Normal}(0, \sigma) \]
\(b_i\) が個人差で, そのばらつき \(\sigma\) もデータから推定します. \(\sigma\) が小さいと推定されれば全員が全体平均に寄せられ, 大きいと推定されれば各自の成績が尊重されます.
これは打者を変量効果とするランダム切片ロジットなので,
stan_glmerで書けます.
rstanと違ってrstanarmはコンパイル済みで出荷されているため,
C++ツールチェーンが無くても動きます.
dat[, batter := factor(FULLNAME)]
fit <- stan_glmer(cbind(HIT, ATBAT - HIT) ~ 1 + (1 | batter),
data = dat, family = binomial("logit"),
chains = 4, iter = 4000, seed = 1234, refresh = 0)dat[, 階層ベイズ := colMeans(posterior_epred(fit))]
sigma_post <- sqrt(as.data.frame(fit)[["Sigma[batter:(Intercept),(Intercept)]"]])\(\sigma\) の事後平均は0.149 (95%区間 0.007 〜 0.385) でした. ロジットスケールで小さい値です.
result <- tibble::tibble(
予測 = c("最尤推定 (序盤の打率)", "全体平均", "階層ベイズ"),
RMSE = c(rmse(dat$AVG, dat$AVG_SEASON),
rmse(group_mean, dat$AVG_SEASON),
rmse(dat$階層ベイズ, dat$AVG_SEASON))
) |>
mutate(`最尤推定に対する比` = RMSE / RMSE[1])
kable(result, digits = 4)| 予測 | RMSE | 最尤推定に対する比 |
|---|---|---|
| 最尤推定 (序盤の打率) | 0.0635 | 1.0000 |
| 全体平均 | 0.0219 | 0.3452 |
| 階層ベイズ | 0.0203 | 0.3201 |
順番は予想通りですが, 中身は予想と違うかもしれません.
- 最尤推定は0.0635.
- ただの全体平均が0.0219で, 最尤推定の2.9分の1です.
- 階層ベイズは0.0203. 全体平均より 7%良いだけです.
つまりこの問題の大半は「個人差を無視する」だけで解けます. 階層ベイズの手柄は, 無視すべきかどうかを自分で決めたことです. \(\sigma\) が大きく推定されるデータなら, 同じモデルが個人差を残します.
どれくらい縮んだか
dat |>
select(FULLNAME, 最尤推定 = AVG, 階層ベイズ, 実際 = AVG_SEASON) |>
tidyr::pivot_longer(-FULLNAME, names_to = "種類", values_to = "打率") |>
mutate(種類 = factor(種類, levels = c("最尤推定", "階層ベイズ", "実際"))) |>
ggplot(aes(x = 種類, y = 打率, group = FULLNAME)) +
geom_line(alpha = 0.4) +
geom_point(size = 2, alpha = 0.7) +
geom_hline(yintercept = group_mean, linetype = "dashed") +
ggtitle("縮小のようす (破線は全体平均)")
ばらつきは18%まで縮みました. 82%は偶然だと判断されたわけです.
推定値の幅は最尤推定が0.143〜0.396, 階層ベイズが0.259〜0.305, 実際が0.273〜0.328です.
順位は変わらない
打数がほぼ同じなので, 縮小は全員に同じようにかかります. つまり順位は動きません.
tibble::tibble(
比較 = c("最尤推定 vs 実際", "階層ベイズ vs 実際", "最尤推定 vs 階層ベイズ"),
順位相関 = c(cor(dat$AVG, dat$AVG_SEASON, method = "spearman"),
cor(dat$階層ベイズ, dat$AVG_SEASON, method = "spearman"),
cor(dat$AVG, dat$階層ベイズ, method = "spearman"))
) |> kable(digits = 3)| 比較 | 順位相関 |
|---|---|
| 最尤推定 vs 実際 | 0.621 |
| 階層ベイズ vs 実際 | 0.615 |
| 最尤推定 vs 階層ベイズ | 0.997 |
縮小推定が直すのは値であって, 並び順ではありません. 「誰が一番打つか」を当てたいなら, 縮小しても情報は増えません. 打数が人によって大きく違う場合は話が別で, そのときは打数の少ない人だけが強く縮むので順位も変わります.
残った誤差はバイアス
階層ベイズでも RMSE は0.0203残りました. 中身を見ます.
err <- dat$階層ベイズ - dat$AVG_SEASON
tibble::tibble(
`平均のずれ (バイアス)` = mean(err),
`ばらつき` = sd(err),
`RMSE` = sqrt(mean(err^2)),
`MSEに占めるバイアスの割合` = mean(err)^2 / mean(err^2)
) |> kable(digits = 4)| 平均のずれ (バイアス) | ばらつき | RMSE | MSEに占めるバイアスの割合 |
|---|---|---|---|
| -0.017 | 0.0114 | 0.0203 | 0.7019 |
誤差の70%はバイアスです. 序盤の全体打率が0.283, シーズンの全体打率が0.300で, 0.017ずれています.
序盤のデータだけを見ている限り, このずれは分かりません. 縮小推定は分散を減らしますが, 中心がずれていれば 全員をずれた場所に集めるだけです.
寒い4月は打率が下がる, といった外の知識を入れるか, 過去のシーズン平均を事前分布に使うかしないと直りません.
まとめ
- 40打数の打率のばらつきは, ほとんどが偶然です.
- 個人差を無視して全体平均と答えるだけで, 序盤の打率より 2.9倍正確でした.
- 階層ベイズは縮小の強さをデータから決めます. ここでは ばらつきの82%を捨てました.
- 縮小しても順位は変わりません. 打数が揃っているためです.
- 残った誤差の大半はバイアスで, 縮小推定では直せません.
同じ考え方は 勝率テーブルの集計 でも使っています. あちらは1試合しかない状態の勝率を モデルの予測に寄せています.
注意
- 20人ぶんの教材用データです.
AVG_SEASONが昇順に並んでいることからも 分かる通り, 説明のために整えられています. 実データでの効き方は 4月の首位打者は, 首位打者になるのか を見てください (2013年4月20日までの規定打席到達者). - 事前分布は
rstanarmの弱情報事前分布です. posterior_epredの事後平均を点推定に使っています.