課題9：階層モデルと shrinkage — 経験ベイズと Gibbs サンプラー

計量経済学II

重要この課題のねらい

第9回でやった partial pooling / shrinkage を、自分の手で最初から最後まで通す。60店舗のクーポン施策を題材に、(1) データを自作し、(2) 経験ベイズのモーメント法で縮小推定を作り、(3) MSE を no pooling と比較し、(4) Gibbs サンプラーを複数 chain で実装・診断してフルベイズの事後分布を出し、(5) 点推定ランキングでなく事後確率でターゲティングする、という流れだ。

• Part 1：60店舗のクーポン施策 DGP を自作（40点）。
• Part 2：経験ベイズ・MSE比較・Gibbs実装・確率ベースのターゲティング（60点）。

配点は合計100点。手を動かせば2〜4時間で終わる分量。完全条件付き分布の式は下に与えるので、まずコードの穴を埋めることから始めよ。焦らず、まずデータを作り、縮小の様子を眺め、それから Gibbs に進もう。

この課題は講義ノート（lecture9.qmd）と同じ構造だが、設定（店舗数・ノイズの大きさ・閾値・文脈）は変えてある。ノートのコードをコピペするだけでは通らないので、自分で理解して埋めること。

---

Part 1：クーポン施策データを自作する（40点）

舞台設定

あるドラッグストアチェーンが、全 \(N=60\) 店舗でクーポン施策を打った。店舗 \(i\) の真のクーポン効果（1人あたり購入額の増分、単位：百円）を \(\theta_i\) とする。効果は店舗ごとに違い、共通の分布から生まれる。各店舗では \(n_i\) 人の来店客について効果が観測されるが、店舗の規模はバラバラ（大型店は客が多く、小型店は少ない）。

モデル（条件付き分布の束）は次の通り。店舗 \(i\)、客 \(j\)：

\[ \theta_i \sim N(\bar\theta,\ \tau^2), \qquad y_{ij} \mid \theta_i \sim N(\theta_i,\ \sigma^2). \]

店舗 \(i\) の標本平均 \(\bar y_i = \frac{1}{n_i}\sum_j y_{ij}\) が、その店の「生のクーポン効果推定値（no pooling 推定）」だ。標本平均の分布は \(\bar y_i \mid \theta_i \sim N(\theta_i,\ \sigma^2/n_i)\) なので、測定分散は \(s_i^2 = \sigma^2/n_i\)（客が少ない店ほど大きい）。

真のパラメータ（あなたはこれを後で推定で回収する）：

記号	意味	値
\(\bar\theta\)	全店平均のクーポン効果	2.0
\(\tau^2\)	店舗間の効果のばらつき（分散）	1.0
\(\sigma^2\)	客1人あたりの観測ノイズ（既知とする）	9.0
\(N\)	店舗数	60

\(\sigma^2\) は既知として扱う（実務では店内のばらつきから別途推定できるが、この課題では簡単のため与える）。\(n_i\) は \(5\)〜\(100\) の範囲でバラバラにする。

問1（15点）：DGP のパラメータと店舗別データを作る

以下のコードを完成させよ。店舗ごとの真の効果 \(\theta_i\)、観測数 \(n_i\)、そして店舗別の標本平均 \(\bar y_i\)（no pooling 推定）を生成する。\(\bar y_i\) は \(N(\theta_i,\ \sigma^2/n_i)\) から直接引いてよい（客ごとの生データを作ってから平均しても同じ）。

# 次元・真のパラメータ
N          <- 60
mu_true    <- 2.0     # 全店平均効果 theta_bar
tau2_true  <- 1.0     # 店舗間分散
sigma2     <- 9.0     # 客1人あたりのノイズ分散（既知）

set.seed(23)          # この課題は seed=23 で固定（採点の再現性のため）

# ---- ここを埋める ----
# 店ごとの観測数 n_i：5..100 の整数を一様に（不均等な店舗規模）。sample() を使う
n_i        <- ______
# 店ごとの真の効果 theta_i ~ N(mu_true, tau2_true)。rnorm の第3引数は「標準偏差」なので sqrt に注意
theta_true <- ______
# 店舗平均 ybar_i の測定分散 s2_i = sigma2 / n_i
s2_i       <- ______
# no pooling 推定 = 店舗別標本平均 ybar_i ~ N(theta_i, s2_i)。rnorm(N, 平均, 標準偏差)
ybar       <- ______

ヒントヒント：なぜ小さい店ほど推定が暴れるのか

\(\bar y_i \sim N(\theta_i,\ \sigma^2/n_i)\) だから、\(n_i\) が小さい店ほど分散 \(\sigma^2/n_i\) が大きく、真の効果 \(\theta_i\) から大きくズレた値が出やすい。\(n_i=6\) の店なら \(s_i^2 = 9/6 = 1.5\)、\(n_i=96\) の店なら \(s_i^2 = 9/96 \approx 0.09\)——16倍以上の差がある。この非対称なノイズが、shrinkage の効きどころだ。

問2（15点）：no pooling 推定と真値の散布図、過分散の確認

生成したデータについて、(a) no pooling 推定 \(\bar y_i\) と真値 \(\theta_i\) の散布図を描き（対角線を引く）、(b) 「no pooling 推定の散らばりが真の効果の散らばりより大きい（過分散）」ことを標準偏差の比較で示せ。

# ---- ここを埋める ----
# (a) no pooling 推定 vs 真値の散布図（点を n_i で色分けすると小さい店の暴れが見える）
df_np <- data.frame(truth = theta_true, nopool = ybar, n_i = n_i)
ggplot(df_np, aes(truth, nopool, color = n_i)) +
  geom_abline(slope = 1, intercept = 0, color = "grey60") +
  geom_point(size = 2, alpha = 0.85) +
  scale_color_viridis_c(name = "観測数 n_i") +
  labs(title = "no pooling 推定 vs 真値", x = "真の効果 theta_i", y = "no pooling 推定 (ybar_i)")

# (b) 過分散の確認
cat("真の効果の標準偏差       :", round(sd(theta_true), 3), "\n")
cat("no pooling 推定の標準偏差:", round(sd(ybar), 3), "\n")

ノートこの散布図の読み方

対角線（推定=真値）から縦に大きく外れている点ほど、推定誤差が大きい。色（\(n_i\)）に注目すると、濃い色＝小さい店ほど対角線から外れやすいはずだ。この「小さい店が外れる」構造が、次の Part 2 で全体平均に縮められる対象になる。

問3（10点）：測定分散と観測数の関係を可視化

\(s_i^2 = \sigma^2/n_i\)（測定分散）を \(n_i\) に対してプロットし、「小さい店ほど測定分散が大きい」ことを図で確認せよ。この \(s_i^2\) が、次節の縮小係数 \(B_i\) を決める。

# ---- ここを埋める ----
ggplot(data.frame(n_i = n_i, s2_i = s2_i), aes(n_i, s2_i)) +
  geom_point(color = "#1f77b4", size = 2, alpha = 0.8) +
  labs(title = "観測数と測定分散", x = "観測数 n_i", y = expression("測定分散 "*s[i]^2*" = sigma2 / n_i"))

---

Part 2：経験ベイズ・Gibbs・ターゲティング（60点）

問4（12点）：経験ベイズのモーメント法で shrinkage

モーメント法で \(\hat{\bar\theta}\) と \(\hat\tau^2\) を推定し、縮小推定 \(\hat\theta_i = B_i\hat{\bar\theta} + (1-B_i)\bar y_i\)、\(B_i = s_i^2/(s_i^2+\hat\tau^2)\) を計算せよ。次に (a) \(\hat{\bar\theta}, \hat\tau^2\) が真値の近くに来ることを確認し、(b) 「no pooling → shrinkage 後」の矢印プロットを描け。

モーメント法の式（講義 §3.1）：\(\hat{\bar\theta} = \frac1N\sum_i \bar y_i\)、\(\hat\tau^2 = \widehat{\mathrm{Var}}(\bar y_i) - \overline{s_i^2}\)。

# ---- 経験ベイズ（モーメント法）----
mu_hat   <- ______     # ybar の平均
tau2_hat <- ______     # Var(ybar) から平均測定分散 mean(s2_i) を引く。負にならないよう max(..., 1e-6)
B_i      <- ______     # 縮小係数 s2_i / (s2_i + tau2_hat)
eb_est   <- ______     # 縮小公式 B_i * mu_hat + (1 - B_i) * ybar

# (a) 真値の回収を確認
cat("mu_hat   =", round(mu_hat, 3),   " （真値", mu_true, "）\n")
cat("tau2_hat =", round(tau2_hat, 3), " （真値", tau2_true, "）\n")

# (b) 矢印プロット（no pooling -> shrinkage 後）。色 = n_i。
df_arrow <- data.frame(nopool = ybar, eb = eb_est, n_i = n_i)
ggplot(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") +
  scale_x_continuous(breaks = c(1, 2), labels = c("no pooling", "shrinkage 後"), limits = c(0.9, 2.3)) +
  scale_color_viridis_c(name = "観測数 n_i") +
  labs(title = "経験ベイズによる縮小", x = NULL, y = "推定された店舗効果")

警告よくあるバグ：tau2_hat が負になる

var(ybar) - mean(s2_i) は、たまたま観測の散らばりが小さいと負になることがある（分散が負なのは無意味）。だから max(..., 1e-6) で下限を置く。もし tau2_hat が0近くに張り付いたら、それは「店舗差がほとんど無い（＝ complete pooling でよい）」という診断でもある。今回の DGP（真の \(\tau^2=1\)）ではそうはならないはずだ。

問5（10点）：MSE 比較表（complete / no pooling / EB）

3つの推定量の MSE \(= \frac1N\sum_i(\hat\theta_i - \theta_i)^2\) を計算し、表にせよ。1つの実現データでは経験ベイズが最小とは限らない点にも注意する。さらに、データ生成を1000回繰り返し、この DGP での期待 MSEをモンテカルロ近似して比較せよ。「平均的に勝つ」と「毎データセットで勝つ」の違いも1文で説明すること。

complete pooling は \(\theta_i=\theta\) を仮定するので、共通効果は測定精度 \(n_i/\sigma^2\) で重みづけた weighted.mean(ybar, n_i) で推定する。これは店舗を等重みにした経験ベイズの単純モーメント推定 mu_hat = mean(ybar) とは区別する。

# ---- ここを埋める ----
mse <- function(est) mean((est - theta_true)^2)
pooled_est <- weighted.mean(ybar, w = n_i) # complete pooling の共通効果MLE
data.frame(
  method = c("complete pooling", "no pooling", "empirical Bayes"),
  MSE    = round(c(mse(pooled_est), mse(ybar), mse(eb_est)), 4)
) |> kable(caption = "3つの推定量の MSE")

# モンテカルロ：1000 データセットで平均 MSE を比較
one_rep <- function() {
  n   <- sample(5:100, N, replace = TRUE)
  s2  <- sigma2 / n
  th  <- rnorm(N, mu_true, sqrt(tau2_true))
  y   <- rnorm(N, th, sqrt(s2))
  muh <- mean(y)
  pool <- weighted.mean(y, w = n)
  t2h <- max(var(y) - mean(s2), 1e-6)
  B   <- s2 / (s2 + t2h)
  eb  <- B * muh + (1 - B) * y
  c(complete = mean((pool - th)^2), nopool = mean((y - th)^2), eb = mean((eb - th)^2))
}
set.seed(20269)
sim <- replicate(1000, one_rep())
round(rowMeans(sim), 4)
cat("EB が no pooling に勝った割合:", round(mean(sim["eb", ] < sim["nopool", ]), 3), "\n")

問6（20点）：Gibbs サンプラーを手実装する

正規階層モデルの Gibbs サンプラーを完成させよ。完全条件付き分布は下に与える（講義 §6.1 の導出結果）。\(\bar\theta\) にはこのデータの尺度では広い正規 prior、\(\tau^2\) には共役計算を学ぶための逆ガンマ prior を置く。異なる初期値から4本の chain を走らせ、traceplot、split-\(\widehat R\)、ESS、MCSE を確認する。

ここで \(\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\) のガンマ分布では rate である。

完全条件付き分布（これを使え）：

\[ \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 + \frac12\sum_i(\theta_i - \bar\theta)^2\right) \quad\Longrightarrow\quad \tau^2 = \frac{1}{\texttt{rgamma}(1, \text{shape}, \text{rate})}. \]

# ---- この関数を完成させる ----
gibbs_hier <- 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)
  N     <- length(ybar)
  theta <- ybar                            # 初期値
  mu    <- if (is.null(init_mu)) mean(ybar) else init_mu
  tau2  <- if (is.null(init_tau2)) max(var(ybar), 0.1) else init_tau2
  keep_theta <- matrix(0, n_draws, N)
  keep_mu    <- numeric(n_draws)
  keep_tau2  <- numeric(n_draws)
  for (it in 1:(n_draws + burn)) {
    # (1) theta_i を引く（ベクトル化）。上の theta_i の完全条件付き分布の v_i, m_i を実装
    v_i   <- ______     # 事後分散：精度 (n_i/sigma2 + 1/tau2) の逆数。分散でなく精度を足す
    m_i   <- ______     # 事後平均：v_i * (データ精度*ybar + prior精度*mu)
    theta <- rnorm(N, m_i, sqrt(v_i))       # rnorm の第3引数は標準偏差 = sqrt(v_i)
    # (2) mu を引く（mean_ と prec は用意済み）
    prec  <- N / tau2 + 1 / kappa2
    mean_ <- (sum(theta) / tau2 + mu0 / kappa2) / prec
    mu    <- ______     # rnorm(1, mean_, 標準偏差)。標準偏差 = sqrt(1/prec)
    # (3) tau2 を引く（逆ガンマ = 1/ガンマ）。shape, rate は用意済み
    shape          <- alpha0 + N / 2
    precision_rate <- beta0 + 0.5 * sum((theta - mu)^2)
    tau2  <- ______     # 1 / rgamma(1, shape = shape, rate = precision_rate)
    if (it > burn) {
      k <- it - burn
      keep_theta[k, ] <- theta
      keep_mu[k]      <- mu
      keep_tau2[k]    <- tau2
    }
  }
  list(theta = keep_theta, mu = keep_mu, tau2 = keep_tau2)
}

fits <- lapply(1:4, function(ch) gibbs_hier(
  ybar, n_i, sigma2, n_draws = 2000, burn = 500,
  seed = 2300 + ch,
  init_mu = mean(ybar) + c(-1, -0.3, 0.3, 1)[ch],
  init_tau2 = c(0.2, 0.6, 1.5, 3)[ch]
))

fit <- list(
  theta = do.call(rbind, lapply(fits, `[[`, "theta")),
  mu = unlist(lapply(fits, `[[`, "mu")),
  tau2 = unlist(lapply(fits, `[[`, "tau2"))
)

続いて、(a) 4本の \(\bar\theta,\tau^2\) の traceplot、(b) 下の簡易関数による split-\(\widehat R\)・ESS・MCSE、(c) 店舗別事後平均と経験ベイズ plug-in 推定の比較、(d) \(\mathrm{IG}(0.01,0.01)\) と \(\mathrm{IG}(3,2)\) の prior 感度分析、を示せ。\(\mathrm{IG}(0.01,0.01)\) は「弱情報」ではなく、あくまで共役計算のデモ設定である。

# (a) 4 chains の traceplot
df_trace <- bind_rows(lapply(1:4, 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") +
  labs(title = "Gibbs の traceplot", x = "反復", y = NULL)

# (b) 簡易診断関数（実務では rank-normalized Rhat と bulk/tail ESS を使う）
basic_mcmc_diag <- function(x) {
  x <- as.matrix(x); n <- nrow(x); m <- ncol(x); n2 <- floor(n / 2)
  sx <- do.call(cbind, lapply(1:m, function(j)
    cbind(x[1:n2, j], x[(n - n2 + 1):n, j])))
  W <- mean(apply(sx, 2, var)); B <- n2 * var(colMeans(sx))
  rhat <- sqrt(((n2 - 1) / n2 * W + B / n2) / W)
  ac <- sapply(1:m, function(j)
    as.numeric(acf(x[, j], lag.max = min(500, n - 1), plot = FALSE)$acf)[-1])
  rho <- rowMeans(ac); le <- 2 * floor(length(rho) / 2)
  ps <- rho[seq(1, le, 2)] + rho[seq(2, le, 2)]
  stop_at <- which(ps <= 0)[1]
  keep <- if (is.na(stop_at)) seq_along(ps) else seq_len(stop_at - 1)
  ess <- min(n * m, n * m / max(1, 1 + 2 * sum(ps[keep])))
  c(split_Rhat = rhat, ESS = ess,
    MCSE_mean = sd(as.vector(x)) / sqrt(ess))
}
rbind(mu = basic_mcmc_diag(sapply(fits, `[[`, "mu")),
      tau2 = basic_mcmc_diag(sapply(fits, `[[`, "tau2"))) |> round(4)

# (c) Gibbs 事後平均 vs 経験ベイズ
post_theta_mean <- 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 事後平均と経験ベイズの比較", x = "経験ベイズ", y = "Gibbs 事後平均")
cat("Gibbs事後平均 と EB の相関:", round(cor(post_theta_mean, eb_est), 5), "\n")

# (d) prior 感度：IG(3,2) でも4 chainsを走らせ、tau2の平均・95%区間を比較する
fits_alt <- lapply(1:4, function(ch) gibbs_hier(
  ybar, n_i, sigma2, alpha0 = 3, beta0 = 2,
  seed = 2400 + ch, init_tau2 = c(0.2, 0.6, 1.5, 3)[ch]
))
tau_default <- unlist(lapply(fits, `[[`, "tau2"))
tau_alt <- unlist(lapply(fits_alt, `[[`, "tau2"))
rbind(
  `IG(0.01,0.01)` = c(mean = mean(tau_default), quantile(tau_default, c(.025, .975))),
  `IG(3,2)` = c(mean = mean(tau_alt), quantile(tau_alt, c(.025, .975)))
) |> round(3)

警告よくあるバグ：rgamma の rate と scale の取り違え

R の rgamma は rate と scale のどちらでも指定できるが、混同すると全く違う分布になる。ここで beta0 は逆ガンマ分布の scale だが、\(G=1/\tau^2\) を引く rgamma には rate = beta0 + 0.5*sum((theta-mu)^2) と渡す。逆数 1/rgamma(...) も忘れない。また v_i の計算で n_i / sigma2 を n_i * sigma2 と書き間違えるミスも多い。

問7（10点）：点推定でなく事後確率でターゲティング

追加予算を「効果の高い店」に投じたい。ここでは 「効果が \(c=3.0\)（百円）を超える確率が80%以上の店」 を選抜し、no pooling の点推定トップ10と比較せよ。どの店が「点推定では上位なのに確率ルールで落ちるか」、その店の \(n_i\) と真値に注目して1〜2文でコメントすること。

# ---- ここを埋める ----
c_thresh <- 3.0
# 各店について P(theta_i > c | data) を、事後サンプル fit$theta が c を超えた割合として計算
# ヒント：fit$theta は (draw 数 x 店数) の行列。列ごとに (fit$theta > c_thresh) の平均を取る
prob_gt  <- ______
sel_prob <- which(prob_gt >= 0.80)          # 確率ルールで選抜
top10_np <- order(ybar, decreasing = TRUE)[1:10]   # no pooling トップ10

cat("確率ルール P(theta>3.0)>=0.8 の店数:", length(sel_prob), "\n")
cat("no pooling トップ10 のうち確率ルールに落ちた店:", length(setdiff(top10_np, sel_prob)), "店\n")

# 点推定では上位、だが確率ルールで落ちた店の正体
dropped <- setdiff(top10_np, sel_prob)
if (length(dropped) > 0) {
  data.frame(
    store = dropped, n_i = n_i[dropped],
    ybar = round(ybar[dropped], 2), post_mean = round(post_theta_mean[dropped], 2),
    prob_gt = round(prob_gt[dropped], 2), truth = round(theta_true[dropped], 2)
  ) |> arrange(desc(ybar)) |> kable(caption = "点推定では上位、確率ルールで落ちた店")
}

ヒントヒント：何が起きているか

閾値 \(c=3.0\) は全店平均 \(\bar\theta=2.0\) より上にある。だから「効果が3.0を超えると確信できる」ためには、点推定が高いだけでなく、不確実性が小さい（＝ \(n_i\) が大きい）必要がある。\(n_i\) が小さい店は、\(\bar y_i\) がたまたま高くても事後分布が横に広く、\(P(\theta_i > 3.0)\) が0.8に届かない。真値 truth を見ると、落ちた店の中に「実は3.0を超えていない」店が混じることがある。確率ルールは、そうした上振れした小さい店への無駄打ちリスクを下げうる。

問8（8点）：経営含意を3行で

上の分析（矢印プロット、MSE、確率ベースのターゲティング）を踏まえ、このドラッグストアの本部がクーポン施策の店舗別評価をどう運用すべきかを、経営含意として3行で書け。とくに「ダッシュボードの店舗別ランキングをそのまま信じることの危険」と「小さい店の扱い」に触れること。

# 判断材料を1つの表にまとめておくと書きやすい
data.frame(
  approach = c("no pooling トップ10", "確率ルール P(>3.0)>=0.8"),
  n_selected = c(10, length(sel_prob))
) |> kable()

ノートこの課題の締め

第9回で学んだ通り、「個別に見たいが個別のデータは少ない」というジレンマは、あらゆるマーケの現場（店舗評価・顧客ターゲティング・営業評価・レコメンド）に現れる。今日あなたが手を動かした shrinkage と Gibbs は、次回（第10回）の顧客パネル分析（BG/NBD）で、購買履歴の異質性推定として再登場する。「異質性を分布でモデル化し、データの薄い個体を全体に縮める」という発想を、自分の道具にしてほしい。

---

提出方法

• この .qmd の穴埋め（______）をすべて埋め、各 chunk の #| eval: false を外して、最後まで render できる状態にする。
• 各問の出力（表・図）とコメントを含めること。
• render した HTML を提出。

重要配点（合計100点）

• 問1（15）：DGP のパラメータ・観測数・店舗別標本平均の生成
• 問2（15）：no pooling 推定 vs 真値の散布図、過分散の確認
• 問3（10）：測定分散と観測数の関係の可視化
• 問4（12）：経験ベイズのモーメント法・縮小推定・矢印プロット
• 問5（10）：MSE比較表とモンテカルロでの平均的優位
• 問6（20）：Gibbs の実装・4 chains・\(\widehat R\)/ESS/MCSE・EB比較・prior 感度
• 問7（10）：事後確率ベースのターゲティングと点推定との比較
• 問8（8）：経営含意（3行）