「店舗別・顧客別に効果を見たいが、個別のデータは少ない」というマーケの根源的ジレンマを、partial pooling（部分プーリング） で解く。complete pooling（全体一律）と no pooling（個別推定）の中間だ。
shrinkage の公式 \(\hat\theta_i = B_i\,\bar\theta + (1-B_i)\,y_i\) 、\(B_i = \dfrac{\sigma_i^2}{\sigma_i^2 + \tau^2}\) を、正規-正規モデルから導く。データが少ない店ほど全体平均に強く引っ張られる 。
この \(\tau^2\) と \(\bar\theta\) をデータから推定するのが経験ベイズ 。ただし、店舗が交換可能で正規分布から生まれるという仮定までデータが自動で正当化するわけではない。ここで本コース初めてベイズ推定 を正面から導入する。
事後分布を数値的に得る道具として Gibbs サンプラー を手で実装する。「条件付き分布から交互に乱数を引くだけ」で、計量Iの「モデル＝条件付き分布の束」がそのまま計算アルゴリズムになる。これがこの回の第二の見せ場だ。
mixed logit（第4回）、hierarchical Bayes conjoint（第6回）、購買履歴モデル（第10回）に共通する「異質性を分布で表し、情報を部分的に共有する」発想をつかむ。CATE は方法によって階層 prior を使わない点も区別する。
前回までの道のりと、今日の位置
第6回でコンジョイントの AMCE（属性効果の平均）を、第7回でA/Bテストの平均処置効果を、第8回でDiD・RDDの因果効果を学んだ。ここまで我々が測ってきたのは、一貫して「平均 」だ。この属性の平均的な部分効用、この施策の平均的な効果、この制度の平均的なインパクト。
だが実務で本当に効くのは、たいてい平均の一歩先にある。第6回では「消費者ごとに部分効用が違う（個人異質性）」ことを指摘し、その本格的な推定は第9回=今日に予告した。第8回では「サブグループに切って多重検定するとまずい」という問題に触れた。今日はこの2つの宿題に、正面から答える。
問いはこうだ。
同じ価格・広告・商品属性に対して、店舗ごと・顧客ごとに反応が違う。この異質性そのもの を推定したい。しかも各店舗・各顧客のデータは少ない。どうやって「個別に見つつ、ノイズに溺れない」を両立するか。
答えの核心は一言、partial pooling（部分プーリング） だ。そしてそれを支える数学が shrinkage（縮小推定） 、それをモデルとして組み立てる枠組みが階層ベイズ・経験ベイズ である。今日はここで、本コースで初めてベイズ推定 を正面から扱う。計量I・IIではずっと頻度論（標本分布・信頼区間）でやってきたが、階層モデルは「ベイズで考えると圧倒的に自然」な対象なので、ここが導入の好機だ。身構えなくていい。ベイズは新しい宗教ではなく、計量Iでやった「モデル＝条件付き分布の束」の素直な延長として入る。
70%くらい分かったら先へ進もう。この回はコース後半の考え方の土台になるので、一度で全部飲み込めなくて当然だ。
1. 実務の痛みから：ダッシュボードの罠
具体例から始める。あなたは全国50店舗を持つ小売チェーンのマーケ担当で、先月、全店で新しいクーポン施策を打った。店舗別の効果（クーポンによる購入額の増分）をダッシュボードに並べ、「今月の効果トップ店舗」ランキングを経営会議で発表した。トップ5店に来月さらに追加予算を投じることに決まった。
ところが翌月、そのトップ5店の効果を見ると、軒並み凡庸 になっている。逆に先月ドベだった店の一部が今月は上位に来ている。何が起きたのか。
これは平均への回帰（regression to the mean） だ。店舗ごとの「観測された効果」は、真の効果 \(\theta_i\) に測定ノイズが乗ったものだ。たまたまノイズが上振れした店が先月のトップに来る。ノイズは翌月には別方向に振れるので、トップ店の効果は真の水準へ引き戻される。「トップ」という選抜そのものが、上振れした店を拾う仕掛けになっている 。
まずこの現象を、データを作って体感しよう。真の効果 \(\theta_i \sim N(2, 1)\) 、各店の観測は \(\theta_i\) にノイズが乗る。とくに小さい店ほど観測が暴れる （サンプルが少ないから）。
# 真の店舗効果 theta_i ~ N(2, tau2=1)、店ごとにノイズ分散が違う <- 50 <- 2.0 <- 1.0 <- 4.0 # 1観測あたりの分散（既知とする） <- sample (5 : 30 , N, replace = TRUE ) # 店ごとの観測数（不均等・小さめ） <- rnorm (N, mu_true, sqrt (tau2_true))<- sigma2 / n_i # 店舗平均 ybar_i の分散 = sigma2/n_i # 先月(month1)と翌月(month2)の観測を、同じ theta から独立に生成 <- rnorm (N, theta_true, sqrt (s2_i))<- rnorm (N, theta_true, sqrt (s2_i))# 先月トップ5を選び、その店の「先月・真値・翌月」を比べる <- order (month1, decreasing = TRUE )[1 : 5 ]data.frame (metric = c ("先月の観測（ダッシュボード表示）" , "その店の真の効果" , "翌月の観測" ),value = round (c (mean (month1[top5]), mean (theta_true[top5]), mean (month2[top5])), 3 )|> kable (caption = "先月トップ5店の、先月・真値・翌月の平均効果" )
先月トップ5店の、先月・真値・翌月の平均効果
先月の観測（ダッシュボード表示）
4.100
その店の真の効果
3.753
翌月の観測
3.408
先月トップ5の「先月の値」は派手に高いが、その店の真の効果 はそれより低く、翌月の観測 は真値のあたりまで下がる。ダッシュボードの数字は上振れを見ていただけだった。これは施策が失敗したのでも、店長が油断したのでもない。ただの統計的な必然だ。
小さいセグメントほど、この暴れは激しくなる。店ごとの観測 \(\bar y_i\) の散らばりは、真の \(\theta_i\) の散らばりより大きくなる （ノイズが上乗せされるから）。図で見よう。
<- data.frame (value = c (theta_true, month1),kind = rep (c ("真の効果 theta_i" , "no pooling 推定 (ybar_i)" ), each = N)ggplot (df_disp, aes (value, fill = kind)) + geom_histogram (alpha = 0.5 , position = "identity" , bins = 20 ) + geom_vline (data = data.frame (kind = c ("真の効果 theta_i" , "no pooling 推定 (ybar_i)" ),m = c (mean (theta_true), mean (month1))),aes (xintercept = m, color = kind), linetype = "dashed" , show.legend = FALSE ) + labs (title = "真の効果の分布 vs 個別推定の分布" ,x = "店舗効果" , y = "店舗数" , fill = NULL ) + theme (legend.position = "top" )
cat ("真の効果の標準偏差 :" , round (sd (theta_true), 3 ), " \n " )cat ("no pooling 推定の標準偏差:" , round (sd (month1), 3 ), " <- 過分散 \n " )
no pooling 推定の標準偏差: 1.091 <- 過分散
no pooling 推定の広がりが真の広がりより大きい。だから「トップとドベの差」を過大評価し、極端な店にばかり目が行く。個別に見たいが、個別のデータは少ない。 これがマーケの根源的ジレンマだ。今日はこれを解く。
同じ罠は人事評価にもある。営業担当を成約率でランキングし、下位の人を叱責する会社は多い。だが担当ごとの案件数はバラバラで、10件しか回していない担当の成約率は激しく暴れる。たまたま下振れした少数案件の担当が「成績不良」として吊るし上げられ、たまたま上振れした担当が表彰される。今日学ぶ shrinkage は、「案件数の少ない担当は全体平均に寄せる」ことで測定ノイズによる評価差を減らす 。ただし、公平性には案件難度・担当割当・評価目的も関わるので、shrinkage だけで公平性が保証されるわけではない。
2. 3つの pooling：complete / no / partial
異質性のあるデータをどう扱うか。極端な2つと、その中間を並べる。店舗 \(i\) の効果 \(\theta_i\) を推定したい。店舗 \(i\) の観測は \(n_i\) 個あり、その標本平均を \(y_i \equiv \bar y_i\) と書く（これが1店舗ぶんの生の推定値だ）。
complete pooling（完全プーリング） ：全店を1つに混ぜ、共通効果 \(\theta\) を全店に使う。共通な個人レベル分散のもとでは \(\hat\theta_{pool}=\sum_i n_i y_i/\sum_i n_i\) 。異質性が完全に死ぬ 。「どの店も同じ効果」と決めつける。no pooling（プーリングなし） ：各店を完全に独立に推定する。\(\hat\theta_i = y_i\) 。ノイズまみれ 。店ごとのデータが少ないと暴れる（§1で見た）。partial pooling（部分プーリング） ：両者の間を取る。データが厚い店は \(y_i\) を信じ、薄い店は \(\bar\theta\) に寄せる。いいとこ取り 。
まず、3つの選択肢から「何を失い、何を得るか」を一段ずつ整理しよう。
complete pooling
全店に同じ精度加重平均
店舗差を無視し、真の異質性を潰す
店舗別の情報も使う
no pooling
店舗平均 \(y_i\)
店舗間で情報を共有せず、小標本のノイズを拾う
推定誤差の大きさを測る
partial pooling
全体平均と \(y_i\) の加重平均
交換可能性と母集団分布を仮定する
ノイズと店舗差から重みを決める
理解の小さな階段は、(1) 店舗ごとの推定対象 \(\theta_i\) を決める、(2) \(y_i\) の測定分散 \(\sigma_i^2\) を求める、(3) 評価基準を全店舗の合計二乗誤差に定める、(4) 正規-正規モデルから事後平均を出す、の4段だ。いきなり「ベイズ」に飛ばず、この順に進む。
言葉だけだとピンとこないので、シミュレーションで3つを比べる。真の分布・no pooling・complete pooling を1枚に描く。
<- weighted.mean (month1, w = n_i) # theta_iが共通なら精度はn_i/sigma2 <- data.frame (store = 1 : N,truth = theta_true,nopool = month1ggplot (df3, aes (x = truth)) + geom_point (aes (y = nopool, color = "no pooling (ybar_i)" ), size = 1.6 , alpha = 0.8 ) + geom_hline (aes (yintercept = pooled_est, color = "complete pooling" ), linetype = "dashed" ) + geom_abline (slope = 1 , intercept = 0 , color = "grey60" ) + scale_color_manual (values = c ("no pooling (ybar_i)" = "#1f77b4" ,"complete pooling" = "#d62728" )) + labs (title = "真の効果 vs 2つの極端な推定" ,x = "真の効果 theta_i" , y = "推定値" , color = NULL ,subtitle = "灰色の対角線が理想（推定=真値）。青は散らばりすぎ、赤は潰れすぎ" ) + theme (legend.position = "top" )
no pooling（青）は対角線のまわりに大きく散らばる——真値をとらえてはいるが、ノイズで上下に暴れる。complete pooling（赤い水平線）は全店を1つの値に潰す——安定だが、店ごとの違いを完全に無視する。理想は対角線 （推定=真値）だが、そこにピタリと乗る推定は作れない。我々にできるのは、青と赤のあいだのどこか に推定を置くことだ。それが partial pooling である。
complete pooling はバイアス大・分散小（全店同じと決めつける）。no pooling はバイアス小・分散大（ノイズを丸ごと拾う）。partial pooling は、このバイアス-分散トレードオフ を、置いた階層モデルのもとで店ごとのデータ量に応じて調整する。データが多い店は no pooling 寄り、少ない店は complete pooling 寄りに、連続的に ブレンドする。
「どれくらいブレンドするか」を決めるのが、次節の shrinkage 係数 \(B_i\) だ。
3. shrinkage の公式（この回の核心）
partial pooling を具体的な式にする。ここが今日の数学的な核だ。設定を正規-正規モデルで書く。
モデル（条件付き分布の束） 。店舗効果 \(\theta_i\) は、共通の分布から生まれる。各店の観測は、その店の効果のまわりにばらつく。
\[
\theta_i \sim N(\bar\theta,\ \tau^2), \qquad
y_i \mid \theta_i \sim N(\theta_i,\ \sigma_i^2).
\]
上の段（\(\theta_i\) の分布）が母集団のばらつき \(\tau^2\) 、下の段（\(y_i\) の分布）が測定のノイズ \(\sigma_i^2\) だ。\(y_i\) は店舗 \(i\) の標本平均なので、\(\sigma_i^2 = \sigma^2/n_i\) ——サンプルが多い店ほど \(\sigma_i^2\) は小さい （測定が正確）。ここが後で効く。
いま \(\bar\theta\) と \(\tau^2\) を既知としよう（後でデータから推定する）。すると店舗 \(i\) の効果 \(\theta_i\) の事後平均 （二乗誤差損失のもとでのベイズ推定量）は、次の形になる。
\[
\hat\theta_i \;=\; B_i\,\bar\theta \;+\; (1 - B_i)\,y_i,
\qquad
B_i \;=\; \frac{\sigma_i^2}{\sigma_i^2 + \tau^2}.
\]
これが shrinkage（縮小）の公式 だ。読み方はこう。
\(B_i\) は縮小の強さ 。\(B_i\) が大きいほど全体平均 \(\bar\theta\) に強く引っ張られる。\(B_i = \dfrac{\sigma_i^2}{\sigma_i^2+\tau^2}\) なので、測定ノイズ \(\sigma_i^2\) が大きい店ほど \(B_i\) が大きい ＝強く縮む。\(\sigma_i^2 = \sigma^2/n_i\) だったから、これは \(n_i\) が小さい（データが少ない）店ほど強く全体平均に寄せる\(y_i\) より、みんなの平均 \(\bar\theta\) を信じた方が当たる。逆に \(n_i \to \infty\) なら \(\sigma_i^2 \to 0\) 、\(B_i \to 0\) で \(\hat\theta_i \to y_i\) 。データが十分な店は no pooling に一致する。
\(\theta_i \sim N(\bar\theta, \tau^2)\) （prior）と \(y_i \mid \theta_i \sim N(\theta_i, \sigma_i^2)\) （尤度）を掛けて、\(\theta_i\) の事後分布を求める。正規×正規は正規になる。指数の肩を \(\theta_i\) について平方完成すればよい。事後分布は
\[
\theta_i \mid y_i \sim N\!\left( \hat\theta_i,\ v_i \right),\qquad
\frac{1}{v_i} = \frac{1}{\sigma_i^2} + \frac{1}{\tau^2},
\]
事後平均は精度（分散の逆数）で重みづけた平均 ：
\[
\hat\theta_i
= v_i\left( \frac{y_i}{\sigma_i^2} + \frac{\bar\theta}{\tau^2} \right)
= \frac{\tau^2}{\sigma_i^2 + \tau^2}\, y_i + \frac{\sigma_i^2}{\sigma_i^2 + \tau^2}\, \bar\theta.
\]
\(B_i = \sigma_i^2/(\sigma_i^2+\tau^2)\) と置けば本文の式そのものだ。「データの精度 \(1/\sigma_i^2\) と prior の精度 \(1/\tau^2\) を足して、それぞれの重みで平均する」——精度が高い方を信じる、というだけの話である。事後分散 \(v_i\) は \(\sigma_i^2\) より必ず小さい（prior が情報を足すから）。この \(v_i\) が後で信用区間を作るのに効く。
3.1 経験ベイズ：\(\bar\theta\) と \(\tau^2\) をデータから決める
\(\bar\theta\) と \(\tau^2\) を既知として話を進めたが、現実には分からない。ここで経験ベイズ（empirical Bayes） の出番だ。「上の段のパラメータ（\(\bar\theta, \tau^2\) ）を、全店のデータから推定してしまう」。モーメント法（積率法）で簡単に出せる。
観測 \(y_i\) の周辺分布を考える。\(y_i = \theta_i + (\text{ノイズ})\) で、\(\theta_i \sim N(\bar\theta, \tau^2)\) 、ノイズの分散は \(\sigma_i^2\) だから、
\[
\mathbb{E}[y_i] = \bar\theta, \qquad
\mathrm{Var}(y_i) = \tau^2 + \sigma_i^2.
\]
したがって、\(\bar\theta\) は \(y_i\) の平均で、\(\tau^2\) は「\(y_i\) の標本分散から、平均的な測定ノイズ \(\overline{\sigma_i^2}\) を引いた残り」で推定できる。
\[
\hat{\bar\theta} = \frac{1}{N}\sum_i y_i, \qquad
\hat\tau^2 = \underbrace{\widehat{\mathrm{Var}}(y_i)}_{\text{観測の散らばり}} - \underbrace{\overline{\sigma_i^2}}_{\text{ノイズの分}}.
\]
「観測の散らばりから、ノイズで説明できる分を差し引いた残りが、真の店舗差 \(\tau^2\) 」という、非常に直感的な式だ。§1で「no pooling は過分散」と言ったのは、まさに \(\mathrm{Var}(y_i) = \tau^2 + \overline{\sigma_i^2} > \tau^2\) だったからで、その過剰分 \(\overline{\sigma_i^2}\) を引き算するのが経験ベイズなのだ。
コードで、真値を回収できるか確認する。
# 経験ベイズ：モーメント法で mu と tau2 を推定 <- mean (month1)<- var (month1) # R の var() は n-1 で割る <- max (var_y - mean (s2_i), 1e-6 ) # 負にならないよう下限を置く cat ("mu_hat =" , round (mu_hat, 3 ), " （真値" , mu_true, "） \n " )cat ("tau2_hat =" , round (tau2_hat, 3 ), " （真値" , tau2_true, "） \n " )cat ("Var(y) =" , round (var_y, 3 ), " のうち、ノイズ分 mean(s2_i) =" ,round (mean (s2_i), 3 ), " を引いた残りが tau2_hat \n " )
Var(y) = 1.191 のうち、ノイズ分 mean(s2_i) = 0.31 を引いた残りが tau2_hat
\(\bar\theta\) も \(\tau^2\) も、真値の近くに戻ってきた。あとはこれを shrinkage の公式に入れるだけだ。
ここで mu_hat = mean(y_i) は「各店舗を1単位として同じ重みで見る」単純なモーメント推定である。全店舗に効果が1つしかないと仮定する complete-pooling モデルの共通効果 MLE は、測定精度 \(1/s_i^2=n_i/\sigma^2\) で重みづけた weighted.mean(y_i, n_i) で、両者は目的も仮定も違う。より効率的な経験ベイズでは周辺尤度を最大化し、\(1/(\tau^2+s_i^2)\) 型の重みを使う方法もある。ここでは縮小公式を見通しよくするためモーメント法を使う。
3.2 shrinkage を可視化する（見せ場）
経験ベイズの推定値 \(\hat\theta_i = B_i\hat{\bar\theta} + (1-B_i)y_i\) を計算し、古典的な矢印プロット で「no pooling → shrinkage 後」の動きを描く。
<- s2_i / (s2_i + tau2_hat) # 縮小係数（店ごと） <- B_i * mu_hat + (1 - B_i) * month1 # 経験ベイズ推定
<- data.frame (store = 1 : N,nopool = month1,eb = eb_est,n_i = n_iggplot (df_arrow) + geom_segment (aes (x = 1 , xend = 2 , y = nopool, yend = eb, color = n_i),arrow = arrow (length = unit (0.12 , "cm" )), alpha = 0.8 ) + geom_hline (yintercept = mu_hat, linetype = "dashed" , color = "grey40" ) + annotate ("text" , x = 2.05 , y = mu_hat, label = "全体平均" , hjust = 0 , color = "grey40" ) + scale_x_continuous (breaks = c (1 , 2 ), labels = c ("no pooling" , "shrinkage 後" ),limits = c (0.9 , 2.4 )) + scale_color_viridis_c (name = "店の観測数 n_i" ) + labs (title = "経験ベイズによる縮小" , x = NULL , y = "推定された店舗効果" ) + theme (legend.position = "right" )
矢印の色（＝観測数 \(n_i\) ）に注目してほしい。濃い色（\(n_i\) が小さい店）ほど矢印が長い ——つまり大きく全体平均に引き寄せられている。データが薄い店の暴れた推定値ほど、強く補正される。データが厚い店（明るい色）はほとんど動かない。\(n_i\) と縮小の強さ \(B_i\) の関係も直接見よう。
ggplot (data.frame (n_i = n_i, B_i = B_i), aes (n_i, B_i)) + geom_point (color = "#1f77b4" , size = 2 , alpha = 0.8 ) + labs (title = "データが少ない店ほど強く縮む" ,x = "店の観測数 n_i" , y = expression ("縮小係数 " * B[i]))
きれいな右下がりだ。\(n_i\) が小さい店ほど \(B_i\) （縮小の強さ）が大きい。これが shrinkage の心臓部——情報量に応じて、信じる度合いを連続的に変える 。
3.3 MSE：本当に改善しているのか
「縮めた方が当たる」を、数字で確かめる。真値 \(\theta_i\) を知っているシミュレーションだからこそできる検算だ。3つの推定量の平均二乗誤差（MSE）\(= \frac{1}{N}\sum_i(\hat\theta_i - \theta_i)^2\) を比べる。
<- function (est) mean ((est - theta_true)^ 2 )data.frame (method = c ("complete pooling" , "no pooling" , "empirical Bayes" ),MSE = round (c (mse (pooled_est), mse (month1), mse (eb_est)), 4 )|> kable (caption = "3つの推定量の MSE（真値との二乗誤差の平均）" )
3つの推定量の MSE（真値との二乗誤差の平均）
complete pooling
1.0937
no pooling
0.3077
empirical Bayes
0.2860
経験ベイズの MSE が最小だ。no pooling より小さく、complete pooling よりずっと小さい。ただし、これは「たまたまこのデータで」かもしれない。データ生成を何度も繰り返して、平均的に 改善するかを確かめよう（計量IIの哲学：モデルを書けばシミュレーションできる、シミュレーションできれば推定量の性質を実験で確かめられる）。
<- function () {<- sample (5 : 30 , N, replace = TRUE )<- sigma2 / n_i<- rnorm (N, mu_true, sqrt (tau2_true))<- rnorm (N, th, sqrt (s2)) # ybar_i ~ N(theta_i, s2_i) <- mean (y)<- weighted.mean (y, w = n_i)<- max (var (y) - mean (s2), 1e-6 )<- s2 / (s2 + t2h)<- B * muh + (1 - B) * yc (complete = mean ((pool - th)^ 2 ),nopool = mean ((y - th)^ 2 ),eb = mean ((eb - th)^ 2 ))set.seed (20269 )<- replicate (1000 , one_rep ()) # 1000 データセット round (rowMeans (sim), 4 )
complete nopool eb
0.9912 0.2969 0.2205
cat (" \n EB が no pooling に勝ったデータセットの割合:" ,round (mean (sim["eb" , ] < sim["nopool" , ]), 3 ), " \n " )
EB が no pooling に勝ったデータセットの割合: 0.943
この DGP では、1000 個のデータセットで平均した MSE は経験ベイズが最小になり、多くのデータセットで no pooling に勝つ。「縮める」のは、平均への回帰を逆手に取った合理的な戦略になりうる のだ。ただし、これはこの正規階層 DGP での平均的な結果であり、経験ベイズが任意のモデル・任意の実現データで必ず勝つという定理ではない。
4. James-Stein の衝撃
前節の「縮めると平均的な合計精度が上がりうる」は、統計学史に残る Stein のパラドクス （Stein 1956、James and Stein 1961）と地続きだ。ただし、定理が保証する対象と、前節の経験ベイズを混同しないように正確に述べる。
\(Y\sim N_p(\theta,\sigma^2I)\) で \(\sigma^2\) が既知、損失が全成分の合計二乗誤差 \(L(\theta,\delta)=\sum_j(\delta_j-\theta_j)^2\) とする。データと独立に決めた固定ターゲット \(a\) へ縮める正部分 James–Stein 推定量は
\[
\delta_{JS}^+(Y)=a+\left[1-\frac{(p-2)\sigma^2}{\lVert Y-a\rVert^2}\right]_+(Y-a),
\qquad [x]_+=\max(x,0).
\]
\(p\ge3\) なら\(Y\) より期待リスク \(R(\theta,\delta)=\mathbb E_\theta[L(\theta,\delta(Y))]\) が小さい。ここで「勝つ」とは、各真値 \(\theta\) を固定してデータを繰り返し生成したときの期待合計二乗誤差 が小さいという意味である。すべての実現データ、すべての成分で誤差が小さくなる、という意味ではない。
一方、縮小先にデータから計算した標本総平均 \(\bar Y\) を使うと、係数は
\[
\delta_L^+(Y)=\bar Y\mathbf 1+
\left[1-\frac{(p-3)\sigma^2}{\sum_j(Y_j-\bar Y)^2}\right]_+
(Y-\bar Y\mathbf 1)
\]
となる（Lindley 型推定量）。こちらで非自明な縮小が始まるのは \(p\ge4\) \(p\ge3\) 」「標本総平均へ縮めるなら \(p\ge4\) 」を分けて覚えよう。
一番有名な実例が、Efron and Morris による打率 の話だ。
問い ：メジャーリーグの選手について、シーズン序盤のわずかな打席の打率から、シーズン全体（残り）の打率を予測したい。各選手の序盤打率をそのまま使うのと、全体平均に縮めるのと、どちらが当たるか。
データ ：1970年シーズン、序盤に約45打席を終えた18人の打者。序盤の打率を「no pooling 推定」とし、その後のシーズン全体の打率を「真値」とみなす。
識別戦略・モデル ：各選手の序盤打率を、全選手の平均打率に向かって James-Stein 型に縮小する。打率は二項なので、分散を安定化する変換（arcsin√）をかけてから正規-正規の枠組みに載せる。縮小の強さは、選手間のばらつきとサンプル誤差の比から自動で決まる（まさに経験ベイズ）。
主要な発見 ：論文の野球例では、Stein 型ルールの18人ぶんの平均二乗予測誤差は、素朴な序盤打率を使う場合の半分未満 だった。これは全選手の損失をまとめた比較であり、「各選手について必ず縮小の方が近い」という主張ではない。
なぜこの回と繋がるか ：これは shrinkage の鮮烈な実例だ。ここでの縮小先は標本総平均なので、対応する非自明な条件は \(p\ge4\) （18人なら満たす）。ただし打率データは正規平均問題そのものではなく、分散安定化変換を使った近似・経験的比較である。序盤の派手な打率に飛びつくと外す、という教訓が §1 のダッシュボードの罠につながる。
パラドクスの「なぜ」を、シミュレーションで体感しておこう。18個の真の値を用意し、少ないサンプルで観測し、縮小あり／なしの合計誤差を比べる。
set.seed (20269 )<- 18 ; nn_js <- 45 <- pmin (pmax (0.265 + 0.06 * rnorm (k_js), 0.15 ), 0.40 ) # 真値ベクトルを固定 <- function () {<- rbinom (k_js, nn_js, p_true_js)<- asin (sqrt ((hits + 0.375 ) / (nn_js + 0.75 ))) # 分散安定化変換 <- 1 / (4 * nn_js) # 変換後の近似分散 <- mean (y)<- sum ((y - gm)^ 2 )<- max (0 , 1 - (k_js - 3 ) * v / S) # 標本総平均へ縮める Lindley 型係数 <- gm + Bl * (y - gm)<- asin (sqrt (p_true_js))c (raw = mean ((y - mu_t)^ 2 ), lindley = mean ((lindley - mu_t)^ 2 ))set.seed (20270 )<- replicate (2000 , js_rep ())cat ("素朴推定の平均MSE :" , round (mean (js_sim["raw" , ]), 5 ), " \n " )cat ("Lindley 型の平均MSE :" , round (mean (js_sim["lindley" , ]), 5 ), " \n " )
Lindley 型の平均MSE : 0.00279
cat ("誤差の比 (縮小/raw) :" , round (mean (js_sim["lindley" , ]) / mean (js_sim["raw" , ]), 3 ), " \n " )cat ("縮小が勝った割合 :" , round (mean (js_sim["lindley" , ] < js_sim["raw" , ]), 3 ), " \n " )
真の18次元ベクトルを固定して観測だけを2000回引き直すと、この設定では縮小推定の平均 MSEが小さくなる一方、縮小が勝った割合 は1ではない。これが「期待リスクで優位」と「毎回勝つ」の違いである。なお二項打率に分散安定化変換を使うここでの計算は、正規平均定理の近似的な実演である。
5. ベイズ推定の導入（コース初）
ここまで「事後平均」「prior」という言葉を使ってきた。この回は本コースでベイズ推定 を正面から扱う初めての回なので、ここで立ち止まって整理する。計量I・IIでずっとやってきた頻度論との関係から入ろう。身構えなくていい。
ベイズの基本等式 はこれだけだ。
\[
\underbrace{p(\theta \mid \text{data})}_{\text{事後分布 posterior}}
\;\propto\;
\underbrace{p(\text{data} \mid \theta)}_{\text{尤度 likelihood}}
\;\times\;
\underbrace{p(\theta)}_{\text{事前分布 prior}}.
\]
「データを見る前の信念（prior）に、データが持つ情報（尤度）を掛けると、データを見た後の信念（posterior）になる」。第2回でやった尤度 は、ここにそのまま登場する。ベイズは尤度を捨てるのではなく、尤度に prior を掛ける だけだ。頻度論の最尤法は「尤度を最大化する1点」を返すが、ベイズは「尤度×prior という分布まるごと」を返す。点でなく分布で答える のがベイズの特徴だ。
「prior なんて恣意的な思い込みじゃないか」という反発は健全だ。階層モデルは判断の一部をデータで学べるようにするが、恣意性を消してくれるわけではない。
\[
\theta_i \sim N(\bar\theta, \tau^2) \quad\Longleftarrow\quad \text{これが }\theta_i\text{ の prior}
\]
店舗効果 \(\theta_i\) の prior は \(N(\bar\theta,\tau^2)\) で、経験ベイズでは \(\bar\theta\) と \(\tau^2\) を全店のデータから推定する。したがって、縮小先と縮小の強さの一部は「他店の実績」から学ばれる。
しかし、(i) 店舗ラベルを入れ替えても同じモデルが妥当という交換可能性 、(ii) 店舗効果の分布が正規、(iii) 店舗間・測定誤差の独立性、などは分析者が置く仮定だ。経験ベイズは同じデータで hyperparameter を推定するため、その推定不確実性を plug-in 区間が十分に反映しないこともある。フルベイズは \(\bar\theta,\tau^2\) に hyperprior を置いて不確実性を伝播させるが、hyperprior の選択という判断は残る。だからモデル診断と prior 感度分析が必要になる。
頻度論とベイズの用語対応も表にしておく。計量Iでやった信頼区間（CI）と、ベイズの信用区間（credible interval）は、似て非なるものだ。
パラメータ \(\theta\) の扱い
固定された未知の定数
確率変数（分布を持つ）
データの扱い
確率変数（標本抽出でばらつく）
観測された定数
推定の答え
点推定＋標準誤差
事後分布まるごと
区間の意味
信頼区間 ：真値を固定して手続きを繰り返せば95%が真値を含む信用区間 ：指定したモデル・prior と観測データのもとで、真値が区間にある事後確率が95%
区間の解釈
「区間がランダム、\(\theta\) は固定」
「\(\theta\) がランダム、区間は固定」
信用区間の方が、素朴な直感（「95%の確率で真値はこの中」）に近い。ただし必ず「置いたモデルと prior のもとで 」という条件が付く。頻度論の信頼区間では、観測後の1本についてこの言い方は厳密には許されない（\(\theta\) は固定定数だから）。マーケの現場で「この店の効果が1.5を超える事後確率は80%です」と言えるのがベイズの利点で、§8で使う。
6. Gibbs サンプラー手実装（第二の見せ場）
経験ベイズは \(\bar\theta, \tau^2\) を点推定して固定した。だが本当は \(\tau^2\) にも不確実性がある（\(\tau^2\) もデータから推定した値だから）。hyperprior のもとでその不確実性も伝播させるのがフルベイズ だ。問題は、\(\theta_1,\ldots,\theta_N, \bar\theta, \tau^2\) という多数のパラメータの同時事後分布が、閉じた形で書けないこと。
ここで Gibbs サンプラー が効く。アイデアは拍子抜けするほど単純だ。
同時分布から直接引くのは難しい。でも「他のパラメータを全部固定したときの、1つのパラメータの条件付き分布」なら簡単に引ける。だったら、パラメータを1つずつ、条件付き分布から順番に引く のを繰り返せばいい。これを大量に回すと、引いた値の集まりが同時事後分布からのサンプルになる。
「条件付き分布から交互に乱数を引くだけ」——これが MCMC（マルコフ連鎖モンテカルロ）の一種、Gibbs サンプリングだ。計量Iで「モデル＝条件付き分布の束」と言ったのを覚えているだろうか。Gibbs は、その条件付き分布の束を、そのまま計算アルゴリズムにする 。モデルの書き方が、そのまま推定の手続きになる。これがこの回の第二の見せ場だ。
6.1 完全条件付き分布を導出する
正規階層モデル（分散 \(\sigma^2\) は既知の単純版）：
\[
y_{ij} \mid \theta_i \sim N(\theta_i, \sigma^2), \qquad
\theta_i \mid \bar\theta, \tau^2 \sim N(\bar\theta, \tau^2), \qquad
\bar\theta \sim N(\mu_0, \kappa^2), \qquad
\tau^2 \sim \text{Inv-Gamma}(\alpha_0, \beta_0).
\]
\(\bar\theta\) には、このデータの単位では広い正規 prior を置く。\(\tau^2\) の逆ガンマは、正規分布の分散パラメータと共役 （掛けると同じ族に戻る）なので、ここでは Gibbs の仕組みを手計算で示すために使う。以下では
\[
p(\tau^2)\propto (\tau^2)^{-\alpha_0-1}\exp(-\beta_0/\tau^2)
\]
と定義する。\(\beta_0\) は逆ガンマ側の scale であり、同値に \(1/\tau^2\sim\mathrm{Gamma}(\alpha_0,\mathrm{rate}=\beta_0)\) である。とくに \(\mathrm{Inv\text{-}Gamma}(0.01,0.01)\) を「弱情報」とは呼ばない。0付近と長い裾に強い形を持ち、少数グループでは結果に影響しうるため、後で感度分析する。\(n_i,\bar y_i\) は店 \(i\) のデータ数と標本平均。
3本の完全条件付き分布はこうなる（導出は下の callout）。
\[
\theta_i \mid \cdots \ \sim\ N(m_i, v_i),\quad
v_i = \left(\frac{n_i}{\sigma^2} + \frac{1}{\tau^2}\right)^{-1},\quad
m_i = v_i\left(\frac{n_i \bar y_i}{\sigma^2} + \frac{\bar\theta}{\tau^2}\right).
\]
\[
\bar\theta \mid \cdots \ \sim\ N\!\left(\ \frac{\sum_i \theta_i / \tau^2 + \mu_0/\kappa^2}{N/\tau^2 + 1/\kappa^2},\ \ \left(\frac{N}{\tau^2} + \frac{1}{\kappa^2}\right)^{-1}\right).
\]
\[
\tau^2 \mid \cdots \ \sim\ \text{Inv-Gamma}\!\left(\alpha_0 + \frac{N}{2},\ \ \beta_0 + \frac{1}{2}\sum_i (\theta_i - \bar\theta)^2\right).
\]
\(\theta_i\) の条件付き平均 \(m_i\) をよく見ると、§3の shrinkage 公式と同じ形 （データと prior 平均の精度加重平均）だ。Gibbs は毎ステップ、この shrinkage を実行しているにすぎない。
\(\theta_i\) の条件付き\(\theta_i\) に関係する項は、尤度 \(\prod_j N(y_{ij}; \theta_i, \sigma^2)\) と prior \(N(\theta_i; \bar\theta, \tau^2)\) だけ。\(\sum_j (y_{ij}-\theta_i)^2 = n_i(\theta_i - \bar y_i)^2 + \text{const}\) を使うと、\(\theta_i\) について正規×正規。§3の平方完成と同じ計算で、精度 \(n_i/\sigma^2 + 1/\tau^2\) 、平均は精度加重平均 \(m_i\) になる。
\(\bar\theta\) の条件付き\(\bar\theta\) に関係するのは \(\prod_i N(\theta_i; \bar\theta, \tau^2)\) と prior \(N(\bar\theta; \mu_0, \kappa^2)\) 。\(\theta_i\) たちを「\(\bar\theta\) のまわりの \(N\) 個の観測」とみなせば、これも正規×正規。精度 \(N/\tau^2 + 1/\kappa^2\) 。\(\kappa^2 \to \infty\) （無情報）なら、単に \(\bar\theta \sim N(\bar\theta_{\text{post}}, \tau^2/N)\) 、\(\bar\theta_{\text{post}} = \frac1N\sum_i\theta_i\) に近づく。
\(\tau^2\) の条件付き\(\tau^2\) に関係するのは \(\prod_i N(\theta_i; \bar\theta, \tau^2) \propto (\tau^2)^{-N/2}\exp(-\frac{1}{2\tau^2}\sum_i(\theta_i-\bar\theta)^2)\) と prior Inv-Gamma\((\alpha_0,\beta_0) \propto (\tau^2)^{-\alpha_0-1}\exp(-\beta_0/\tau^2)\) 。掛けると再び逆ガンマの形になる（共役）。形状 \(\alpha_0 + N/2\) 、逆ガンマ側の尺度 \(\beta_0 + \frac12\sum_i(\theta_i-\bar\theta)^2\) 。
逆ガンマからのサンプリングは、ガンマの逆数で作る：\(\tau^2 = 1/G\) 、\(G \sim \text{Gamma}(\text{shape}=\alpha_0+N/2,\ \text{rate}=\beta_0+\frac12\sum(\theta_i-\bar\theta)^2)\) 。R の rgamma に渡すのは、精度 \(G=1/\tau^2\) の rate である。
6.2 コードは20行
完全条件付き分布さえ導ければ、あとは「順番に引く」だけ。ベクトル化して書く。この回の DGP と同じデータ（§1 の month1＝店ごとに \(n_i\) 観測ある想定の標本平均）で回す。
<- function (ybar, n_i, sigma2, n_draws = 2000 , burn = 500 ,mu0 = 0 , kappa2 = 1e6 ,alpha0 = 0.01 , beta0 = 0.01 ,seed = NULL , init_mu = NULL , init_tau2 = NULL ) {if (! is.null (seed)) set.seed (seed)<- length (ybar)<- ybar # 初期値 <- if (is.null (init_mu)) mean (ybar) else init_mu<- if (is.null (init_tau2)) max (var (ybar), 0.1 ) else init_tau2<- matrix (0 , n_draws, N)<- numeric (n_draws)<- numeric (n_draws)for (it in 1 : (n_draws + burn)) {# (1) theta_i を引く（店ごと、ベクトル化） <- 1 / (n_i / sigma2 + 1 / tau2)<- v_i * (n_i * ybar / sigma2 + mu / tau2)<- rnorm (N, m_i, sqrt (v_i))# (2) mu（全体平均）を引く <- N / tau2 + 1 / kappa2<- (sum (theta) / tau2 + mu0 / kappa2) / prec<- rnorm (1 , mean_, sqrt (1 / prec))# (3) tau2（店舗間分散）を引く：逆ガンマ = 1/ガンマ <- alpha0 + N / 2 <- beta0 + 0.5 * sum ((theta - mu)^ 2 )<- 1 / rgamma (1 , shape = shape, rate = precision_rate)# burn-in を過ぎたら保存 if (it > burn) {<- it - burn<- theta<- mu<- tau2list (theta = keep_theta, mu = keep_mu, tau2 = keep_tau2)<- lapply (1 : 4 , function (ch) {gibbs_hier (n_draws = 2000 , burn = 500 ,seed = 20269 + ch,init_mu = mean (month1) + c (- 1 , - 0.3 , 0.3 , 1 )[ch],init_tau2 = c (0.2 , 0.6 , 1.5 , 3 )[ch]# 事後要約には4本を結合するが、診断までは必ず鎖を分けて保持する <- list (theta = do.call (rbind, lapply (fits, ` [[ ` , "theta" )),mu = unlist (lapply (fits, ` [[ ` , "mu" )),tau2 = unlist (lapply (fits, ` [[ ` , "tau2" ))
1本の鎖の中身は rnorm と rgamma を交互に呼ぶ短いループだ。ただし、収束を1本だけで判断すると「同じ場所に閉じ込められた」ことを見落とすので、異なる初期値から4本の鎖 を走らせる。
6.3 収束と事後分布を診る
まずtraceplot （反復ごとのサンプルの軌跡）で、4本が同じ領域を行き来しているかを見る。ただし見た目だけでは足りないので、split-\(\widehat R\) 、有効サンプルサイズ（ESS）、事後平均の Monte Carlo 標準誤差（MCSE）も確認する。
<- bind_rows (lapply (seq_along (fits), function (ch) {data.frame (iter = seq_along (fits[[ch]]$ mu), chain = factor (ch),mu = fits[[ch]]$ mu, tau2 = fits[[ch]]$ tau2)|> pivot_longer (c (mu, tau2), names_to = "param" , values_to = "value" )ggplot (df_trace, aes (iter, value, color = chain)) + geom_line (linewidth = 0.25 , alpha = 0.8 ) + facet_wrap (~ param, scales = "free_y" ,labeller = as_labeller (c (mu = "全体平均 mu" , tau2 = "店舗間分散 tau2" ))) + labs (title = "4本の Gibbs chain の traceplot" , x = "反復" , y = NULL , color = "chain" )
# 教育用の基本診断。実務では posterior パッケージの rank-normalized Rhat と bulk/tail ESS を使う。 <- function (x) {<- as.matrix (x) # 行=反復、列=chain <- nrow (x); m <- ncol (x); n2 <- floor (n / 2 )<- do.call (cbind, lapply (1 : m, function (j) {cbind (x[1 : n2, j], x[(n - n2 + 1 ): n, j])<- mean (apply (split_x, 2 , var))<- n2 * var (colMeans (split_x))<- (n2 - 1 ) / n2 * W + B / n2<- sqrt (var_plus / W)<- min (500 , n - 1 )<- sapply (1 : m, function (j) {as.numeric (acf (x[, j], lag.max = lag_max, plot = FALSE )$ acf)[- 1 ]<- rowMeans (ac)<- 2 * floor (length (rho) / 2 )<- rho[seq (1 , last_even, by = 2 )] + rho[seq (2 , last_even, by = 2 )]<- which (pair_sum <= 0 )[1 ]<- if (is.na (first_nonpos)) seq_along (pair_sum) else seq_len (first_nonpos - 1 )<- max (1 , 1 + 2 * sum (pair_sum[keep]))<- min (n * m, n * m / tau_int)<- sd (as.vector (x)) / sqrt (ess)c (split_Rhat = rhat, ESS = ess, MCSE_mean = mcse)<- sapply (fits, ` [[ ` , "mu" )<- sapply (fits, ` [[ ` , "tau2" )<- rbind (mu = basic_mcmc_diag (mu_chains),tau2 = basic_mcmc_diag (tau2_chains))round (diag_table, 4 ) |> kable (caption = "MCMC診断（4 chains）" )
MCMC診断（4 chains）
mu
1.0010
4652.442
0.0023
tau2
1.0001
3156.384
0.0045
\(\widehat R\) は1に近いほど鎖間・鎖内のばらつきが整合し、ESS は自己相関を割り引いた実効的な標本数、MCSE は有限回の MCMC による事後平均の数値誤差だ。目安として split-\(\widehat R<1.01\) を確認し、ESS と MCSE が目的の精度に十分かを見る。上の関数は概念を見せる簡易版で、実務では rank-normalized split-\(\widehat R\) と bulk/tail ESS を使う。次に、\(\bar\theta\) と \(\tau^2\) の事後分布 を経験ベイズの点推定と比較する。
<- data.frame (mu = fit$ mu, tau2 = fit$ tau2) |> pivot_longer (everything (), names_to = "param" , values_to = "value" )<- data.frame (param = c ("mu" , "mu" , "tau2" , "tau2" ),value = c (mu_hat, mu_true, tau2_hat, tau2_true),type = c ("経験ベイズ" , "真値" , "経験ベイズ" , "真値" )ggplot (df_post, aes (value)) + geom_histogram (aes (y = after_stat (density)), bins = 40 , fill = "#1f77b4" , alpha = 0.5 ) + geom_vline (data = vlines, aes (xintercept = value, color = type, linetype = type), linewidth = 0.8 ) + facet_wrap (~ param, scales = "free" ,labeller = as_labeller (c (mu = "全体平均 mu" , tau2 = "店舗間分散 tau2" ))) + scale_color_manual (values = c ("経験ベイズ" = "#d62728" , "真値" = "grey30" )) + scale_linetype_manual (values = c ("経験ベイズ" = "dashed" , "真値" = "dotted" )) + labs (title = "超パラメータの事後分布" , x = NULL , y = "密度" , color = NULL , linetype = NULL ) + theme (legend.position = "top" )
図 6: 超パラメータの事後分布。破線は経験ベイズの点推定、点線は真値。この設定での近さを比較する。
cat ("mu : 事後平均" , round (mean (fit$ mu), 3 )," / 経験ベイズ" , round (mu_hat, 3 ), " / 真値" , mu_true, " \n " )
mu : 事後平均 2.145 / 経験ベイズ 2.17 / 真値 2
cat ("tau2 : 事後平均" , round (mean (fit$ tau2), 3 )," / 経験ベイズ" , round (tau2_hat, 3 ), " / 真値" , tau2_true, " \n " )
tau2 : 事後平均 0.912 / 経験ベイズ 0.881 / 真値 1
この DGP、この標本サイズ、この prior では、Gibbs の事後平均と経験ベイズの点推定が近くなる。これは一般定理ではない。少数グループ、境界に近い \(\tau^2\) 、異なる hyperprior、モデルのミススペシフィケーションでは差が大きくなりうる。フルベイズは hyperparameter の不確実性を事後分布へ伝播させる一方、その結果は hyperprior にも依存する。
6.4 prior 感度分析
\(\mathrm{Inv\text{-}Gamma}(0.01,0.01)\) は共役計算のデモには便利だが、「ほぼ無情報」ではない。\(\tau^2\) の結論が hyperprior の選択でどれだけ動くか、同じデータで比較する。
<- data.frame (prior = c ("IG(0.01, 0.01): 共役デモ" , "IG(3, 2)" , "IG(5, 4)" ),alpha0 = c (0.01 , 3 , 5 ),beta0 = c (0.01 , 2 , 4 )<- lapply (seq_len (nrow (prior_grid)), function (g) {<- lapply (1 : 4 , function (ch) gibbs_hier (n_draws = 1200 , burn = 300 ,alpha0 = prior_grid$ alpha0[g], beta0 = prior_grid$ beta0[g],seed = 30000 + 100 * g + ch,init_mu = mean (month1) + c (- 1 , - 0.3 , 0.3 , 1 )[ch],init_tau2 = c (0.2 , 0.6 , 1.5 , 3 )[ch]<- unlist (lapply (fs, ` [[ ` , "tau2" ))<- unlist (lapply (fs, ` [[ ` , "mu" ))data.frame (prior = prior_grid$ prior[g], mu_mean = mean (md),tau2_mean = mean (td), tau2_lo = unname (quantile (td, 0.025 )),tau2_hi = unname (quantile (td, 0.975 )))|> bind_rows ()|> mutate (across (where (is.numeric), ~ round (.x, 3 ))) |> kable (caption = "逆ガンマ hyperprior を変えた感度分析" )
逆ガンマ hyperprior を変えた感度分析
IG(0.01, 0.01): 共役デモ
2.151
0.907
0.513
1.476
IG(3, 2)
2.148
0.871
0.521
1.401
IG(5, 4)
2.149
0.885
0.547
1.378
この例で結論が似ていれば、それはデータに店舗が十分あり prior の影響が小さいという結果 であって、事前に保証された性質ではない。店舗数を減らすと差が広がりやすい。実務では \(\tau\) への half-normal なども候補にし、事前予測チェックと感度分析で尺度が妥当か確認する（その場合はこの共役 Gibbs 以外のサンプラーが必要になる）。
6.5 店舗別事後平均と経験ベイズ shrinkage を比較する
各店の \(\theta_i\) の事後平均（Gibbs サンプルの平均）を、§3の経験ベイズ plug-in 公式と比較する。この設定では近いはずだが、同じ推定量ではない。
<- colMeans (fit$ theta)ggplot (data.frame (eb = eb_est, gibbs = post_theta_mean), aes (eb, gibbs)) + geom_abline (slope = 1 , intercept = 0 , color = "grey60" ) + geom_point (color = "#1f77b4" , size = 2 , alpha = 0.8 ) + labs (title = "Gibbs 事後平均と経験ベイズ shrinkage の比較" ,x = "経験ベイズ closed form" , y = "Gibbs 事後平均" )cat ("Gibbs事後平均 と EB の相関:" , round (cor (post_theta_mean, eb_est), 5 ), " \n " )
Gibbs事後平均 と EB の相関: 0.99998
cat ("最大のズレ:" , round (max (abs (post_theta_mean - eb_est)), 4 ), " \n " )
この設定では点が対角線の近くに並ぶ。Gibbs の各更新で使う条件付き平均が shrinkage の形を持つことと、最終的なフルベイズ事後平均が経験ベイズ plug-in 推定と同一であること は別である。後者の近さは、この DGP・標本サイズ・hyperprior での経験的結果だ。MCMC は条件付き分布の束から同時事後分布を数値的に近似するアルゴリズムであり、診断で近似誤差を確かめる必要がある。
6.6 キャタピラープロット：店舗別の信用区間
フルベイズの最大の利点は、点推定だけでなく各店の事後分布まるごと が手に入ること。店舗別の95%信用区間を、キャタピラープロット（毛虫プロット）で並べる。事後分布の \(\theta_i\) サンプルから、2.5%・97.5%分位点を取るだけだ。
<- apply (fit$ theta, 2 , quantile, 0.025 )<- apply (fit$ theta, 2 , quantile, 0.975 )<- data.frame (store = 1 : N, mean = post_theta_mean, lo = ci_lo, hi = ci_hi,truth = theta_true, n_i = n_i|> arrange (mean) |> mutate (rank = row_number ())ggplot (df_cat, aes (rank, mean)) + geom_hline (yintercept = mu_hat, linetype = "dashed" , color = "grey60" ) + geom_errorbar (aes (ymin = lo, ymax = hi), color = "#1f77b4" , width = 0 , alpha = 0.7 ) + geom_point (size = 1.2 , color = "#1f77b4" ) + geom_point (aes (y = truth), color = "#d62728" , size = 1 , alpha = 0.8 ) + labs (title = "店舗別の事後平均と95%信用区間（毛虫プロット）" ,subtitle = "青=事後平均と信用区間、赤=真値、破線=全体平均" ,x = "店舗（事後平均の順にソート）" , y = "店舗効果 theta_i" )
赤点（真値）の大半が青い信用区間に収まっている。データが少ない店ほど区間が広い（不確実性が正直に表現される）——no pooling では暴れた点推定しか出なかったが、フルベイズは「この店はよく分からない」を区間の広さで教えてくれる。この不確実性の可視化 こそ、点推定ランキングにはできない芸当だ。
信用区間の頻度的な被覆も、定義を明示してシミュレーションで確認しよう。ここでは最初に生成した \(\theta_1,\ldots,\theta_N\) と \(n_i\) を固定 し、同じ真値のまわりで測定誤差だけを繰り返し発生させる。したがって調べるのは fixed-truth repeated-sampling coverage であり、\(\theta_i\) を引き直す random-effects の周辺被覆でも、hyperparameter まで引き直す prior-predictive coverage でもない。
<- function (seed) {set.seed (seed)<- rnorm (N, theta_true, sqrt (s2_i))<- lapply (1 : 4 , function (ch) gibbs_hier (n_draws = 800 , burn = 250 ,seed = 100000 + 10 * seed + ch,init_mu = mean (y) + c (- 1 , - 0.3 , 0.3 , 1 )[ch],init_tau2 = c (0.2 , 0.6 , 1.5 , 3 )[ch]<- do.call (rbind, lapply (fs, ` [[ ` , "theta" ))<- apply (draws, 2 , quantile, 0.025 )<- apply (draws, 2 , quantile, 0.975 )>= lo & theta_true <= hi<- sapply (1 : 30 , coverage_rep) # 行=店舗、列=反復データセット cat ("fixed-truth 反復標本での平均被覆率:" , round (mean (cover), 3 ), " \n " )
fixed-truth 反復標本での平均被覆率: 0.943
cat ("店舗別被覆率の範囲:" , paste (round (range (rowMeans (cover)), 3 ), collapse = " -- " ), " \n " )
全店舗・全反復を平均した被覆率は、モデルが整合的なら0.95付近になることが期待されるが、30反復なので Monte Carlo 誤差もある。また、全体平均が0.95に近くても各店舗を固定した条件付き被覆が一様に0.95とは限らない 。端の \(\theta_i\) 、小さい \(n_i\) 、hyperprior の選択では店舗別被覆がずれうる。各回で \(\theta_i\) だけを固定 hyperparameter の分布から引き直すなら random-effects の周辺被覆 、hyperparameter まで hyperprior から引き直すなら prior-predictive coverage であり、いずれも今の fixed-truth 条件付き評価とは別物である。
7. マーケの異質性推定の統一風景
今日の「異質性を表し、情報を個体間で共有する」という発想は、この講義の様々なモデルに現れる。ただし、すべてが同じ階層 prior を使うわけではない。一枚の表にまとめよう。
mixed logit / RC logit
消費者ごとの選好係数 \(\beta_i\)
\(\beta_i \sim N(\beta, \Sigma)\) 第4回
hierarchical Bayes conjoint
個人ごとの部分効用 \(\beta_i\)
\(\beta_i \sim N(\bar\beta, V)\) を Gibbs で推定第6回
正規階層モデル（今日）
店舗効果 \(\theta_i\)
\(\theta_i \sim N(\bar\theta, \tau^2)\) 第9回
BG/NBD の購買異質性
顧客ごとの購買率 \(\lambda\)
\(\lambda \sim \text{Gamma}(r, \alpha)\) 第10回（次回）
CATE（異質処置効果）
共変量ごとの処置効果 \(\tau(x)\)
causal forest・meta-learner・BART 等、方法ごとに異なる正則化
第11回
mixed logit、HB conjoint、正規階層モデル、BG/NBD は「個体パラメータを母集団分布から生じると置く」という階層構造を明示的に共有する。第4回の \(\beta_i\sim N(\beta,\Sigma)\) は今日のモデルの多変量版で、第6回の HB conjoint はその構造を数値的に推定する。一方、一般の CATE 推定は同じ prior を置かない 。causal forest や meta-learner は木・回帰・交差適合などで正則化し、Bayesian CATE（例：BART）の一部だけが prior を使う。共通するのは「異質性をそのまま細分化して過学習せず、何らかの partial pooling／regularization を使う」という設計原理である。
問い ：顧客の購買履歴 は、ターゲティング（誰にクーポンを送るか）にとってどれだけの価値があるのか。履歴を使わず全員に一律で送るのと、履歴から個人の反応を推定して送り分けるのとで、利益はどれだけ違うか。
データ ：Nielsen のツナ缶スキャンパネルから抽出した400世帯（平均約13回の購買）。ブランド選択、価格・販促環境、人口統計を使う。
モデル ：顧客ごとのランダム係数を持つ多項プロビット 選択モデルを使い、係数の母集団分布を人口統計変数にも結びつける。潜在効用をデータ拡張した Gibbs sampler で、個人係数と母集団パラメータを推定する。履歴が短い顧客は母集団分布へ強く縮小され、履歴が長い顧客では個人の選択履歴がより効く。この事後分布をクーポン配布の意思決定に使う。
主要な発見 ：短い購買履歴を使うターゲット・クーポンの純収益増分は一律配布の約2.5倍で、1回の購買機会だけを使う場合でも一律配布より約50%大きいと報告した。モデル依存の意思決定計算ではあるが、情報集合の価値を金銭で比較した点が重要である。
なぜこの回と繋がるか ：選択モデル自体は今日の正規観測モデルより複雑だが、「履歴の薄い顧客ほど母集団情報を借り、履歴が溜まるにつれて個別化する」という shrinkage の発想は同じである。これは §8 のコールドスタートにもつながる。
問い ：コンジョイント調査で、1人あたりの設問数を減らしたい（回答者の負担を軽くしたい）。だが設問が少ないと、個人ごとの部分効用は推定できないのではないか。少ない設問数で、どこまで個人の選好を復元できるか。
データ・設定 ：回答者が仮想的な商品・サービスの full profile を評価する metric conjoint （評定データ）であり、選択データではない。各回答者が答えた full-profile evaluation から calibration 用プロファイルの回答をランダムに削り、短い調査票を模した reduced design を作る。残しておいた validation responses を使い、部分効用の回復と予測を検証する。
識別戦略・モデル ：個人ごとの部分効用 \(\beta_i\) に階層 prior を置く hierarchical Bayes conjoint。Gibbs サンプラーで、全回答者の情報を借りながら各個人の \(\beta_i\) を推定する。設問が少ない個人は、他の回答者から推定された母集団分布（prior）に縮小される。
主要な発見 ：階層ベイズを使うと、個人あたりの calibration profiles 数が個人パラメータ数を下回っていても （個人別最小二乗が存在しない状況でも）、個人の部分効用の異質性をかなり良く回復できることを示した。全回答者でプールした情報が各個人の推定を支え、削除せずに残した validation responses の予測にも使える。ここで予測しているのは評定反応であって、選択確率ではない。
なぜこの回と繋がるか ：第6回のコンジョイントで「個人異質性を推定したい、でも1人あたりのデータが少ない」と言った、まさにその問題の解決策がこれだ。この論文の metric conjoint は線形・正規の構造なので、今日の Gibbs を行列に拡張した形に近い。ただし分散共分散行列など追加のパラメータがあり、単なるコード置換ではない。
問い ：マーケティングのモデルは、消費者異質性をどう扱うべきか。異質性を無視した集計モデルは、何を、どのように取り違えるのか。
データ・設定 ：複数のスキャンパネルデータセットにわたる、ブランド選択・購買のモデリング。個人レベルの選好・反応の分布を推定する。
識別戦略・モデル ：階層ベイズを一貫した枠組みとして提示する。個人ごとの選好パラメータに母集団分布（prior）を置き、Gibbs で個人と母集団を同時推定する。連続的な異質性（正規分布など）と離散的な異質性（潜在クラス）を、同じベイズの枠組みで比較・整理した。
主要な論点 ：集計モデルは個人の異質性を「平均的な1人」に押し込めるため、非線形な選択モデルでは価格反応や市場シェア予測を取り違えうる。とくにターゲティングのように個人差を使う意思決定では、異質性分布と個人別不確実性の扱いが重要になる。
なぜこの回と繋がるか ：この論文（と Rossi-Allenby-McCulloch の教科書）は、今日の講義全体の設計図だ。「消費者は違う、その違いを階層ベイズで推定し、意思決定に使う」という思想を、quantitative marketing の中心に据えた。今日の shrinkage・Gibbs・ターゲティングは、すべてこの路線の上にある。
8. 意思決定への接続：ランキングでなく確率で
最後に、事後分布を「意思決定」に繋ぐ。§1のダッシュボードの罠を思い出そう。点推定でトップ店をランキングすると、上振れした店を拾ってしまう。フルベイズなら、もっと賢い選び方ができる。
「上位の店に追加投資したい」を、点推定のランキング でなく、事後確率 \(P(\theta_i > c \mid \text{data})\) で判断する。「効果が閾値 \(c\) を超える確率が十分高い店」だけを選ぶ。これは事後分布のサンプルから、\(\theta_i > c\) となった割合を数えるだけで出せる（§5で言った「確率で語れる」がベイズの強みだ）。
ここでは説明のため「効果が \(c=3.0\) を超える事後確率が80%以上の店」を選び、no pooling の点推定トップ10と比べる。\(c=3.0\) と0.80は例示的な基準 で、実務では追加投資の便益・費用、誤選択の損失、予算制約から決める。事後確率だけで効用関数は自動的に決まらない。
<- 3.0 <- colMeans (fit$ theta > c_thresh) # P(theta_i > c | data) <- which (prob_gt >= 0.80 ) # 確率ルールで選抜 <- order (month1, decreasing = TRUE )[1 : 10 ] # no pooling 点推定トップ10 cat ("確率ルール P(theta>3.0)>=0.8 で選ばれた店数:" , length (sel_prob), " \n " )
確率ルール P(theta>3.0)>=0.8 で選ばれた店数: 5
cat ("no pooling トップ10 のうち、確率ルールに落ちた店:" ,length (setdiff (top10_np, sel_prob)), "店 \n " )
no pooling トップ10 のうち、確率ルールに落ちた店: 5 店
# 点推定では上位なのに、確率ルールで落ちた店の正体を表にする <- setdiff (top10_np, sel_prob)data.frame (store = dropped,n_i = n_i[dropped],ybar_nopool = round (month1[dropped], 2 ),post_mean = round (post_theta_mean[dropped], 2 ),prob_gt_c = round (prob_gt[dropped], 2 ),truth = round (theta_true[dropped], 2 )|> arrange (n_i) |> kable (caption = "no pooling では上位、だが確率ルールで落ちた店（n_i の小さい順）" )
no pooling では上位、だが確率ルールで落ちた店（n_i の小さい順）
20
7
3.91
3.22
0.64
4.10
47
7
3.58
3.02
0.51
4.27
23
12
3.18
2.89
0.41
2.62
29
12
3.14
2.86
0.39
2.56
8
24
3.19
3.03
0.53
3.03
# 落ちた店の観測数と、全店の観測数の中央値を比べる（前者が小さいはず） cat ("落ちた店の n_i の中央値:" , median (n_i[dropped])," / 全店の n_i の中央値:" , median (n_i), " \n " )
落ちた店の n_i の中央値: 12 / 全店の n_i の中央値: 17
no pooling のトップ10のうち、いくつかの店が確率ルールで落ちる。その落ちた店を \(n_i\) の小さい順に並べた表を見てほしい。落ちた店の観測数の中央値は、全店の中央値より小さい はずだ——データが少ないのに \(\bar y_i\) がたまたま上振れして点推定ランキングの上位に来た店たちである。事後分布で見ると、これらの店は「効果が本当に \(c\) を超えている確信が持てない」（\(P(\theta_i > c)\) が0.8に届かない）。真値 truth の列を見ると、\(\bar y_i\) の派手さのわりに真の効果が控えめな店（\(c=3.0\) に届いていない店）も混じっている。確率ルールは、上振れした小さい店への無駄な投資リスクを下げうる 。
図でも見よう。点推定 \(\bar y_i\) と事後確率 \(P(\theta_i > c)\) の関係を、\(n_i\) で色分けする。
<- data.frame (ybar = month1, prob = prob_gt, n_i = n_i)ggplot (df_tgt, aes (ybar, prob, color = n_i)) + geom_hline (yintercept = 0.80 , linetype = "dashed" , color = "grey50" ) + geom_vline (xintercept = c_thresh, linetype = "dotted" , color = "grey50" ) + geom_point (size = 2.2 , alpha = 0.85 ) + scale_color_viridis_c (name = "観測数 n_i" ) + annotate ("text" , x = min (month1), y = 0.82 , label = "選抜ライン P=0.8" ,hjust = 0 , color = "grey40" , size = 3 ) + labs (title = "点推定でなく事後確率で選ぶ" ,x = "no pooling 点推定 (ybar_i)" , y = "P(theta_i > c | data)" )
同じ点推定 \(\bar y_i\) でも、\(n_i\) が小さい店（濃い色）ほど事後確率が低くかもしれない 」と「効果が高いと確信できる 」は違う。追加投資は後者にすべきだ。ベイズは、この区別を確率という共通言語で扱えるようにする。
レコメンドのコールドスタート ：新商品・新店舗には履歴がない。素朴に「クリック率0/0」と扱うと計算できないし、「たまたま最初の1人がクリックしたらCTR100%」と暴れる。解決策の一つは、新アイテムを全体平均から始め、データが溜まるにつれて縮小を緩める こと。多腕バンディットの Thompson sampling も事後分布から行動を無作為化するが、使う尤度・prior は問題ごとに異なる。
ベイジアンA/BテストツールのUI ：最近のA/Bテストツールは「Bが勝つ確率85%」といった表示を出す。あの数字の正体は、今日の \(P(\theta_B > \theta_A \mid \text{data})\) そのものだ。頻度論のp値（「帰無仮説のもとで、観測された検定統計量以上に極端な統計量が出る確率」）と違い、「Bが勝っている確率」を直接答える。第7回のA/Bテストを、ベイズで見ると事後分布の比較になる。UIの裏で走っているのは、今日書いたような事後分布計算である。
まとめ
「個別に見たいが個別のデータは少ない」というジレンマは、partial pooling で解く。complete pooling（異質性が死ぬ）と no pooling（ノイズまみれ）の中間。
shrinkage の公式 \(\hat\theta_i = B_i\bar\theta + (1-B_i)y_i\) 、\(B_i = \sigma_i^2/(\sigma_i^2+\tau^2)\) 。データが少ない店ほど全体平均に強く縮む 。この DGP では no pooling より平均 MSE が下がるが、どの実現データでも勝つとは限らない。\(\bar\theta,\tau^2\) をデータから推定するのが経験ベイズ 。古典的な正規平均問題では、固定ターゲットへの James–Stein 縮小は \(p\ge3\) 、標本総平均への Lindley 型縮小は非自明に \(p\ge4\) 。保証は期待合計二乗誤差についてである。ベイズ ：事後 ∝ 尤度 × prior。階層モデルでも交換可能性・分布形・hyperprior は分析者の仮定で、prior の判断は消えない。信用区間の確率解釈にも「指定したモデル・prior のもとで」が付く。Gibbs サンプラー ：完全条件付き分布から交互に乱数を引く。異なる初期値の複数 chain を走らせ、traceplot だけでなく \(\widehat R\) ・ESS・MCSE を診る。経験ベイズとの近さは、この DGP・標本サイズ・prior での結果である。意思決定は点推定ランキングでなく \(P(\theta_i>c\mid\text{data})\) で。多くの異質性モデルは partial pooling／regularization の発想を共有するが、一般の CATE が同じ階層 prior を持つわけではない。
宿題
ブラウザ実験室（playground9.html ）で、今回のシミュレーションの数値を自分でいじって遊んでみよう。インストール不要、ブラウザだけでRが動く（初回ロードのみ30秒〜1分かかる）。
本文 §3 のコードで、sigma2（1観測あたりのノイズ分散）を 4.0 から 1.0 に下げて実行し直せ。縮小係数 \(B_i\) はどう変わるか。ノイズが小さくなると、全体平均に縮む必要は増えるか減るか。直感と照らし合わせよ。
§6 の Gibbs で、n_draws を 2000 から 200 に減らして traceplot を見よ。事後分布のヒストグラムはどう変わるか。「サンプルが足りない」とはどういう見た目かを体感する。
ChatGPT に「経験ベイズとフルベイズの違いを、階層モデルの文脈で説明して」と聞き、本文 §5・§6 の説明と突き合わせよ。とくに「\(\tau^2\) の不確実性の扱い」がどう違うかに注目。
ChatGPT に「固定ターゲットへの James–Stein 推定量の \(p\ge3\) と、標本総平均への Lindley 型推定量の \(p\ge4\) はなぜ違うのか」と聞いてみよ。さらに「勝つ」が期待合計二乗誤差の比較であり、毎データセット・各成分の勝利ではないことを自分の言葉で説明せよ。
coding 課題は assignment9.qmd を参照。60店舗のクーポン施策 DGP を自作し、経験ベイズと Gibbs を自分の手で実装し、確率ベースのターゲティングまで通す。この回のすべてを自分の手で回す課題だ。
次回予告
第10回は、今日の階層的な発想を購買履歴 に適用する。顧客ごとに「どれくらいの頻度で買うか（購買率 \(\lambda\) ）」「いつ離脱するか」が違う。BG/NBD モデル では購買率の顧客間異質性をガンマ分布で表し、離脱過程も別にモデル化するため、今日の正規階層モデルと「分布を置き換えるだけ」ではない。ただし、異質性を母集団分布で表し、短い履歴について集団情報を借りるという発想は共通している。
なお次回、ガンマ分布のパラメータとして \(\alpha\) が登場するが、これは需要推定の回（第2〜6回）で使った価格係数 \(\alpha\) とは全くの別物 だ。BG/NBD の標準的な記法では \(\alpha\) は rate パラメータ（scale ではない） で、R の rgamma(shape = r, rate = alpha) に対応する。
参考文献
Rossi, P. E., Allenby, G. M., and McCulloch, R. (2005). Bayesian Statistics and Marketing. Wiley.
Allenby, G. M., and Rossi, P. E. (1999). “Marketing Models of Consumer Heterogeneity.” Journal of Econometrics , 89(1–2), 57–78.
Rossi, P. E., and Allenby, G. M. (2003). “Bayesian Statistics and Marketing.” Marketing Science , 22(3), 304–328.
Lenk, P. J., DeSarbo, W. S., Green, P. E., and Young, M. R. (1996). “Hierarchical Bayes Conjoint Analysis: Recovery of Partworth Heterogeneity from Reduced Experimental Designs.” Marketing Science , 15(2), 173–191. https://doi.org/10.1287/mksc.15.2.173
Rossi, P. E., McCulloch, R. E., and Allenby, G. M. (1996). “The Value of Purchase History Data in Target Marketing.” Marketing Science , 15(4), 321–340.
Efron, B., and Morris, C. (1975). “Data Analysis Using Stein’s Estimator and Its Generalizations.” Journal of the American Statistical Association , 70(350), 311–319.
Efron, B., and Morris, C. (1977). “Stein’s Paradox in Statistics.” Scientific American , 236(5), 119–127.
Stein, C. (1956). “Inadmissibility of the Usual Estimator for the Mean of a Multivariate Normal Distribution.” Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability , 1, 197–206.
James, W., and Stein, C. (1961). “Estimation with Quadratic Loss.” Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability , 1, 361–379.
Gelman, A., and Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
