simulate_assignment <- function(seed=2026101,N=240,Ttrain=12,Thold=8) {
set.seed(seed)
beta_known <- c(battery12=.5,nc=.8,major_brand=.3)
mu_true <- log(.15)-.45^2/2
tau_true <- .45
theta_true <- rnorm(N,mu_true,tau_true)
make_tasks <- function(T) {
id <- rep(seq_len(N),each=T); nr <- length(id)
# price,battery12,nc,major_brandをnr行2列の行列として作る
# quality=x'beta、choiceを0:2の選択として作る
list(id=id,price=price,battery12=battery12,nc=nc,
major_brand=major_brand,quality=quality,choice=choice,T=T)
}
list(train=make_tasks(Ttrain),test=make_tasks(Thold),N=N,
beta_known=beta_known,theta_true=theta_true,mu_true=mu_true,tau_true=tau_true)
}課題10:個人別需要とクーポン配分
計量経済学II
この課題のねらい
階層logitから得た個人の事後分布を使い、クーポンの配布ルールを選ぶ。推定した需要と施策の利益が、コードのどこでつながるかを確認しよう。目安は2〜4時間、配点は合計100点である。
Part 1:データを作る(30点)
問1:個人差のある選択履歴(20点)
イヤホンA/B/どちらも買わないの選択を、240人について作る。推定用12問、holdout用8問を生成し、seedは2026101とする。
\[u_{ijt}=x_{ijt}'\beta-\exp(\theta_i)p_{ijt}+\varepsilon_{ijt},\quad \theta_i\sim N(\log(0.15)-0.45^2/2,0.45^2).\]
属性係数は \(\beta=(0.5,0.8,0.3)\) とし、学習時にも既知と置く。価格は8・12・16千円、3つの属性はそれぞれ確率1/2で0か1を取り、個人の好みとは独立に割り当てる。outsideの効用の定数項は0である。講義ノートの simulate_cbc() を参考に、下の骨格を埋めよう。関数が返すリストの形式を変えないこと。
問2:履歴と価格の確認(10点)
全体の選択比率と、個人別の「どちらも買わない」比率を示せ。真の価格感応度の分布も描くこと。価格が顧客の好みと独立に割り当てられていることが、後のクーポン比較になぜ必要か、2〜3文で説明せよ。
Part 2:推定から施策・評価へ(70点)
問3:階層logitの推定(20点)
下の推定関数を使い、問1の推定用履歴だけから事後drawを得よ。\(\mu,\tau\) の真値・事後平均・95%信用区間・split \(\widehat R\) を表にし、holdoutの負の対数予測確率も計算する。推定用の回答をholdout評価に混ぜてはいけない。
第9回と同じ関数である。MCMC自体を書き直す必要はない。
# === CBC_DGP ===
choice_prob <- function(quality, price, alpha) {
v <- quality - price * alpha
vmax <- pmax(0, v[, 1], v[, 2])
ev <- exp(v - vmax)
cbind(outside = exp(-vmax), ev) / (exp(-vmax) + rowSums(ev))
}
individual_loglik <- function(theta, tasks, N) {
v <- tasks$quality - tasks$price * exp(theta[tasks$id])
vmax <- pmax(0, v[, 1], v[, 2])
logden <- vmax + log(exp(-vmax) + rowSums(exp(v - vmax)))
chosen <- numeric(length(tasks$choice))
inside <- which(tasks$choice > 0)
chosen[inside] <- v[cbind(inside, tasks$choice[inside])]
as.vector(rowsum(chosen - logden, tasks$id, reorder = FALSE))
}
fit_hb <- function(dat, seed = 202609, n_iter = 16000, burn = 4000, thin = 4,
prior_mu = log(0.2), prior_mu_sd = 1.5,
prior_shape = 2, prior_rate = 0.5) {
keep <- seq(burn + thin, n_iter, by = thin)
one_chain <- function(chain) {
set.seed(seed + chain)
theta <- rep(prior_mu + c(-0.4, 0.4)[chain], dat$N)
mu <- mean(theta); tau2 <- 0.5
ll <- individual_loglik(theta, dat$train, dat$N)
theta_draw <- matrix(NA_real_, length(keep), dat$N)
mu_draw <- tau_draw <- numeric(length(keep))
accepted <- numeric(dat$N); saved <- 0
for (iteration in seq_len(n_iter)) {
proposal <- theta + rnorm(dat$N, 0, 0.7)
ll_new <- individual_loglik(proposal, dat$train, dat$N)
logratio <- ll_new - ll + dnorm(proposal, mu, sqrt(tau2), log = TRUE) -
dnorm(theta, mu, sqrt(tau2), log = TRUE)
accept <- log(runif(dat$N)) < logratio
theta[accept] <- proposal[accept]; ll[accept] <- ll_new[accept]
if (iteration > burn) accepted <- accepted + accept
variance_mu <- 1 / (dat$N / tau2 + 1 / prior_mu_sd^2)
mean_mu <- variance_mu * (sum(theta) / tau2 + prior_mu / prior_mu_sd^2)
mu <- rnorm(1, mean_mu, sqrt(variance_mu))
tau2 <- 1 / rgamma(1, shape = prior_shape + dat$N / 2,
rate = prior_rate + sum((theta - mu)^2) / 2)
if (iteration %in% keep) {
saved <- saved + 1
theta_draw[saved, ] <- theta
mu_draw[saved] <- mu; tau_draw[saved] <- sqrt(tau2)
}
}
list(theta = theta_draw, mu = mu_draw, tau = tau_draw,
acceptance = accepted / (n_iter - burn))
}
chains <- lapply(1:2, one_chain)
list(theta = rbind(chains[[1]]$theta, chains[[2]]$theta),
mu = c(chains[[1]]$mu, chains[[2]]$mu), tau = c(chains[[1]]$tau, chains[[2]]$tau),
chain = rep(1:2, each = length(keep)), chains = chains,
settings = c(n_iter = n_iter, burn = burn, thin = thin))
}
# === PREDICTION ===
predict_existing <- function(fit, tasks) {
answer <- matrix(0, length(tasks$id), 3)
for (b in seq_len(nrow(fit$theta))) {
answer <- answer + choice_prob(tasks$quality, tasks$price, exp(fit$theta[b, tasks$id]))
}
answer / nrow(fit$theta)
}
predict_new <- function(fit, quality, price, n_integration = 80) {
# 母集団パラメータの不確実性と新しい個人の異質性の両方を積分。
nodes <- qnorm((seq_len(n_integration) - 0.5) / n_integration)
answer <- matrix(0, nrow(price), 3)
for (b in seq_along(fit$mu)) {
for (z in nodes) {
answer <- answer + choice_prob(quality, price, exp(fit$mu[b] + fit$tau[b] * z))
}
}
answer / (length(fit$mu) * n_integration)
}
log_score <- function(prob, choice) -mean(log(pmax(prob[cbind(seq_along(choice), choice + 1)], 1e-12)))
split_rhat <- function(chain1, chain2) {
n <- floor(length(chain1) / 2)
z <- cbind(chain1[1:n], tail(chain1, n), chain2[1:n], tail(chain2, n))
W <- mean(apply(z, 2, var)); B <- n * var(colMeans(z))
sqrt(((n - 1) * W / n + B / n) / W)
}問4:クーポンを選ぶ(20点)
Aの属性による効用は1.6、Bは0.8、通常価格はどちらも12千円とする。Aだけを販売する企業の原価は3千円、候補の値引きは0・1・2・3千円である。配信費用はない。
- 各顧客について、施策ごとの購買確率を事後drawで平均せよ。
- 値引き後のマージンを掛け、最適な一律クーポンと個人別クーポンを選べ。
- 配布なし・全員3千円・最適な一律・個人別の期待利益を比較せよ。
- 真の個人パラメータで同じルールを評価し、事後期待利益による比較と区別せよ。
利益を最大にするルールの選択に真値を使わないこと。なお、仮想回答と実購買で同じモデルが成り立つという教育用の仮定を置く。
問5:固定したルールを新しい実験で評価する(15点)
問4の個人別ルールを固定する。seedを2026102にし、各顧客を確率1/2でルール適用群か配布なし群に割り当て、新しい1回の購買結果を生成せよ。実現利益の平均差と標準誤差、95%信頼区間を示せ。標本が小さいため区間が広くても、それだけで誤りとはしない。
問6:経営への含意(15点)
「通常価格で買う確率が最も高い顧客から順にクーポンを配ればよい」という提案に、3行で回答せよ。購入予測、施策による増分利益、仮想回答から実購買への移行という3点を含めること。
提出物と採点基準
実行できるQMDと、そこから作ったHTMLを提出する。問1はDGPとデータ形式、問2は可視化と価格割当の説明、問3は事後分布とholdout、問4は積分と利益計算、問5は無作為割当と推論、問6は意思決定の解釈を評価する。MCMCや実験の1回の結果が真値にぴったり一致することは求めない。