Playground 9：shrinkageの実験室

計量経済学II（ブラウザ上でRが動きます）

作者

Kei Ikegami

重要このページの使い方
  • このページのRコードは、あなたのブラウザの中で動いている（WebAssembly + webR）。インストールは一切不要。
  • 初回だけ、ページ読み込みに30秒〜1分かかる（R本体とパッケージのダウンロード）。2回目以降はキャッシュで速くなる。
  • 各セルのコードは自由に書き換えて「Run Code」で実行できる。何をどう壊しても誰にも迷惑はかからない。おかしくなったらページを再読み込みすれば元に戻る。
  • 元ネタは講義ノート（lecture9.html）。ここでは「shrinkage（縮小推定）がいつ・どれくらい効くのか」を、数値をいじって体感することに集中する。

実験1：shrinkageの綱引き — 矢印はどこまで引かれるか

店舗効果の推定は「自分のデータ \(y_i\)」と「全体平均 \(\bar\theta\)」の綱引きだ。真の店舗効果は \(\theta_i \sim N(2, \tau^2)\) から生まれ、各店は \(n_i\) 個の観測を持つので、店舗平均 \(y_i\) のノイズ分散は \(\sigma^2 / n_i\) になる。経験ベイズはこの \(y_i\)（no pooling 推定）を、次の式で全体平均へ引き寄せる。

\[ \hat\theta_i = B_i \, \bar\theta + (1 - B_i) \, y_i, \qquad B_i = \frac{\sigma_i^2}{\sigma_i^2 + \tau^2} \]

縮小の強さ \(B_i\) は「店ごとの測定ノイズ」と「店舗間の本物のばらつき」の比で決まる。データが語るなら聞き、語らないなら全体に頼る——この標語が矢印の長さにそのまま現れる。\(\tau\) と \(n_i\) をいじって確かめよう。

ノートやってみよう
  • tau <- 3（店舗どうしが全然似ていない世界）にすると、矢印はほとんど動かなくなる。逆に tau <- 0.2 ではほぼ全店が全体平均へ潰される。表の mean shrinkage B と合わせて、縮小が弱まる・強まるを確認しよう。
  • n_range <- c(50, 100)（どの店もデータが厚い）にすると、\(\tau\) が同じでも矢印が短くなることを確認しよう。c(2, 8) にすれば逆のことが起きる。
  • 矢印が長い店と短い店の色（\(n_i\)）を見比べよう。同じ \(\tau\)・同じ \(\sigma^2\) のもとで、店ごとに縮小の強さを変えているものは何か。「データが語るなら聞き、語らないなら全体に頼る」は、\(B_i\) の式のどの部分に対応しているだろうか。

実験2：MSE対決 — 3つのpooling、どれが当たるか

綱引きの絵がきれいでも、「本当に当たるのか」は別の問題だ。真値 \(\theta_i\) を知っているシミュレーションだからこそ、真値との二乗誤差（MSE）で答え合わせができる。データ生成を200回繰り返し、complete pooling（共通効果を \(n_i\) で精度加重推定）・no pooling（各店の \(y_i\) をそのまま）・partial pooling（経験ベイズ）のMSEを平均する。win_rate は「そのデータセットで最小MSEだった回数の割合」だ。どんな世界でどれが勝つのか、\(\tau\) と店舗数・観測数をいじって探そう。

ノートやってみよう
  • tau <- 5（店舗が全く似ていない世界）にすると、partial pooling の優位は消えるか。no pooling とのMSEの差と win_rate を、tau <- 1 のときと見比べよう。
  • tau <- 0.1 かつ N <- 100（全店ほぼ同じ効果）にすると、complete pooling が勝ちに来ることを確認しよう。次に n_range <- c(100, 200) として、no pooling と partial pooling の差がほぼ消えることも見よう。
  • partial pooling が勝つ回と負ける回の両方を探そう。経験ベイズは \(\hat\tau^2\) を通じて縮小を調整するが、有限標本での hyperparameter 推定誤差やモデルのずれがあるので、毎データセットで勝つ保証はない。表の MSE（反復平均）と win_rate（実現ごとの勝率）を混同しないこと。

実験3：ミニGibbsサンプラー — 事後分布を自分の目で見る

経験ベイズは \(\tau^2\) を1点に固定したが、その値自体にも不確実性がある。それも伝播させるのがフルベイズで、道具は講義ノート§6の Gibbs サンプラーだ（ここでは軽量版：各 chain draws 800・burn-in 200）。異なる初期値から4本を走らせ、traceplotに加えて split-\(\widehat R\)・ESS・MCSE を見る。さらに逆ガンマ prior を変えて感度を確かめる。

ここで \(\mathrm{Inv\text{-}Gamma}(\alpha_0,\beta_0)\) は \(p(\tau^2)\propto(\tau^2)^{-\alpha_0-1}\exp(-\beta_0/\tau^2)\)。\(\beta_0\) は逆ガンマの scale で、\(1/\tau^2\) を引く rgamma では rate として渡す。

ノートやってみよう
  • N <- 5（店舗が5つしかない世界）にすると、\(\tau^2\) の事後分布はどれくらい横に広がるか。表の95%区間の幅を N <- 30 のときと見比べよう。店舗間のばらつきは「店舗の数」からしか学べないので、少数グループでは階層モデルでも \(\tau^2\) をほぼ特定できない——これが階層モデルの限界だ。
  • tau_true <- 0.3 にすると、事後分布は0の近くに寄る。traceplotで \(\tau^2\) がときどき0のすぐ上まで潜る様子も観察しよう。tau_true <- 2 なら分布ごと右へ動く。
  • \(\tau^2\) の推定を良くするには、各店の観測数 \(n_i\) を増やすのと店舗数 \(N\) を増やすのと、どちらが本質的に効くのか。コード中の c(5, 30) を c(100, 300) に変えた場合と、N を増やした場合とで、95%区間の縮み方を比べて考えよう。
  • traceplotの4色が重なるか、split-\(\widehat R<1.01\) か、ESS と MCSE が十分かをセットで確認しよう。上の診断関数は教育用の簡易版で、実務では rank-normalized split-\(\widehat R\) と bulk/tail ESS を使う。
  • \(\mathrm{IG}(0.01,0.01)\) は「弱情報」ではなく共役計算のデモ。N <- 5 にすると \(\mathrm{IG}(3,2)\) との差が広がるか確認しよう。prior を変えて結論が動くなら、その不確実性を報告する。

次に読むもの

  • 講義ノート本体：lecture9.html（shrinkage公式の導出、James-Stein、ベイズ入門、Gibbsと信用区間の完全版）
  • coding課題：assignment9.qmd（60店舗のクーポン施策DGPを自作し、経験ベイズとGibbsをゼロから実装して、確率ベースのターゲティングまで通す）