第7回の選択課題に戻り、個人の価格係数を階層ベイズで推定する。
母集団の平均とばらつき、個人ごとの係数には、それぞれ事後の不確実性がある。
個人の回答が少ない場合は、母集団の情報を使っても、その人の好みを正確に言い当てられるとは限らない。
予測は、係数の事後平均をlogitに代入するだけでなく、事後分布にわたって選択確率を平均する。
既存顧客の次の選択と、履歴がない新規顧客の選択を区別して予測する。
前回とのつながり・今回のゴール
前回は、正規分布のモデルでベイズ推定を学んだ。データから個人の情報を取り出しつつ、母集団の情報も使うと、少ない観測に振り回されにくくなる。ただし、そこで観測していたのは連続的な評価スコアだった。
今回は第7回のイヤホンの選択課題に戻る。「Aを選んだ」「どちらも買わなかった」という回答から、価格への反応を人ごとに推定しよう。同じ人の次の選択を、全員同じ係数のモデルよりよく予測できるか が今回のゴールである。推定結果は次回、誰にどのクーポンを配るかを考えるときにも使う。
1. 同じ選択モデルに、個人の分布を加える
第7回と同じく、2つのイヤホンと「どちらも買わない」という外部選択肢を考える。価格の単位は千円で、調査では8、12、16から無作為に選ぶ。属性は、長時間バッテリー、ノイズキャンセリング、大手ブランドの3つである。
\[
u_{ijt}=x_{ijt}'\beta-\exp(\theta_i)p_{ijt}+\varepsilon_{ijt},\qquad j=A,B,
\qquad u_{i0t}=\varepsilon_{i0t}.
\]
\(\alpha_i=\exp(\theta_i)>0\) が価格係数で、値上げは効用を下げる。個人の分布は
\[
\theta_i\mid\mu,\tau^2\sim N(\mu,\tau^2)
\]
とする。したがって価格係数は対数正規分布になる。
今回のコードでは、仕組みを見通せるように、属性の係数 \(\beta=(0.5,0.8,0.3)\) は既知に固定し、価格係数の個人差だけを推定する 。第7回の推定値を無条件で真値とみなしたわけではない。実際の分析では、属性係数や定数項についても推定し、その不確実性を含める必要がある。ここで出す区間は、固定した属性係数が正しいという条件付きの区間である。
選択確率と、繰り返された選択の尤度
第2回のlogitと同じ計算により、\(\theta_i\) が与えられたときの選択確率は
\[
P_{ijt}(\theta_i)=
\frac{\exp\{x_{ijt}'\beta-\exp(\theta_i)p_{ijt}\}}
{1+\sum_{k=A,B}\exp\{x_{ikt}'\beta-\exp(\theta_i)p_{ikt}\}}.
\]
外部選択肢の確率は、分子を1にした式である。同じ人の8問に共通するのは \(\theta_i\) で、各問の誤差は条件付きで独立とする。観測された選択を \(y_{it}\) と書けば、個人の尤度は
\[
L_i(\theta_i)=\prod_{t=1}^{T_i}P_{i,y_{it},t}(\theta_i).
\]
積の中で、同じ \(\theta_i\) を使っているところが重要である。回答ごとに別の顧客を引き直したモデルでは、その人の一貫した好みを学習できない。
価格が高い商品を1度選んだだけでは、価格をあまり気にしない人なのか、たまたまその商品の属性がよかったのか分からない。属性と価格を変えた複数の課題を使うと、候補となる係数を比較できる。それでも8問で残る曖昧さは、点推定で消したことにせず、事後分布として持っておきたい。
2. データを作り、推定用と検証用を分ける
300人それぞれについて、推定用の8問と検証用の8問を作る。第7回と同じ生成手順・乱数seedを使う。同じ人の検証用の選択にも同じ真の \(\theta_i\) を使うが、商品属性と選択の誤差は新しく発生させる。
真の母集団パラメータは \(\tau=0.6\) 、\(\mu=\log(0.12)-0.6^2/2\) とする。このため、母集団における平均価格係数は \(E[\alpha_i]=0.12\) である。真の個人係数を推定へ渡すことはせず、後の答え合わせだけに使う。
第7回と共通のCBCデータ生成
<- function (quality, price, alpha) {<- quality - price * alpha<- pmax (0 , v[, 1 ], v[, 2 ])<- exp (v - vmax)cbind (outside = exp (- vmax), ev) / (exp (- vmax) + rowSums (ev))<- function (seed = 202607 , N = 300 , Ttrain = 8 , Thold = 8 ) {set.seed (seed)<- c (battery12 = 0.5 , nc = 0.8 , major_brand = 0.3 )<- log (0.12 ) - 0.6 ^ 2 / 2 <- 0.6 <- rnorm (N, mu_true, tau_true)<- function (T) {<- rep (seq_len (N), each = T)<- length (id)<- matrix (sample (c (8 , 12 , 16 ), nr * 2 , replace = TRUE ), nr, 2 )<- matrix (rbinom (nr * 2 , 1 , 0.5 ), nr, 2 )<- matrix (rbinom (nr * 2 , 1 , 0.5 ), nr, 2 )<- matrix (rbinom (nr * 2 , 1 , 0.5 ), nr, 2 )<- beta_known[1 ] * battery12 + beta_known[2 ] * nc + beta_known[3 ] * major_brand<- choice_prob (quality, price, exp (theta_true[id]))<- vapply (seq_len (nr), function (r) sample (0 : 2 , 1 , prob = prob[r, ]), integer (1 ))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)<- simulate_cbc ():: kable (data.frame (sample = c ("推定" , "検証" ),customers = dat$ N, tasks_per_customer = c (dat$ train$ T, dat$ test$ T),outside_share = c (mean (dat$ train$ choice == 0 ), mean (dat$ test$ choice == 0 ))),digits = 3 , col.names = c ("用途" , "人数" , "1人の課題数" , "非購入割合" ))
推定
300
8
0.442
検証
300
8
0.437
<- data.frame (alpha = exp (dat$ theta_true))ggplot (truth_plot, aes (alpha)) + geom_histogram (bins = 30 , fill = "#1f77b4" , color = "white" ) + labs (x = "価格係数の真値（千円あたり）" , y = "人数" , title = "生成した価格反応の個人差" )
これは調査の中で価格と属性を無作為化したデータである。実際の購買履歴で企業が価格を選んでいる場合には、第3回の価格内生性を別に考える必要がある。階層ベイズに変えるだけで、その問題まで解決するわけではない。
3. 母集団と個人の事後分布
前回と同じように、母集団パラメータにも事前分布を置く。
\[
\mu\sim N\{\log(0.2),1.5^2\},\qquad
\tau^2\sim IG(2,0.5).
\]
この数値は価格の単位が千円であることを踏まえた設定で、真の \(\mu,\tau\) を入れたものではない。\(IG(a,b)\) は前回と同じく、shape \(a\) 、rate \(b\) のガンマ変数の逆数である。例えばこの事前分布では \(E[\tau^2]=0.5\) なので、どの程度の個人差を想定しているかも明示できる。
全員の回答を \(D\) と書くと、事後分布は
\[
p(\theta_1,\ldots,\theta_N,\mu,\tau^2\mid D)
\propto
\left\{\prod_i L_i(\theta_i)\,\phi(\theta_i;\mu,\tau^2)\right\}
p(\mu)p(\tau^2).
\]
\(L_i\) は本人の回答から来る情報、\(\phi\) は母集団から来る情報である。母集団の分布自体も全員のデータから更新されるので、他の人の回答が間接的にその人の推定に使われる。
前回の正規モデルでは、個人の事後分布を正規分布の式で書けた。今回は尤度がlogitなので、その式は使えない。ただし、尤度と事前分布の積に比例するという考え方は同じである。
Metropolisの更新をGibbsの中に入れる
\(\mu,\tau^2\) を固定したとき、個人 \(i\) の更新では次の手順を使う。
現在の \(\theta_i\) の近くに、新しい候補 \(\theta_i^*\) を対称な正規分布から発生させる。
候補と現在の「尤度×母集団密度」を比べる。
比率が1より大きければ採用し、1未満ならその比率の確率で採用する。
この更新がMetropolis法である。悪い候補も時々採用するため、最大値だけに向かう最適化とは違う。続いて、今の全員の \(\theta_i\) を使って \(\mu,\tau^2\) を条件付き分布から更新する。
対称な提案分布では、採用確率は
\[
\min\left\{1,
\frac{L_i(\theta_i^*)\phi(\theta_i^*;\mu,\tau^2)}
{L_i(\theta_i)\phi(\theta_i;\mu,\tau^2)}\right\}
\]
になる。計算では小さな確率の積を直接掛けず、対数尤度を足して比較する。\(\mu\) の条件付き分布は正規、\(\tau^2\) は逆ガンマで、前回と同じ形である。
<- function (theta, tasks, N) {<- tasks$ quality - tasks$ price * exp (theta[tasks$ id])<- pmax (0 , v[, 1 ], v[, 2 ])<- vmax + log (exp (- vmax) + rowSums (exp (v - vmax)))<- numeric (length (tasks$ choice))<- which (tasks$ choice > 0 )<- v[cbind (inside, tasks$ choice[inside])]as.vector (rowsum (chosen - logden, tasks$ id, reorder = FALSE ))<- 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 ) {<- seq (burn + thin, n_iter, by = thin)<- function (chain) {set.seed (seed + chain)<- rep (prior_mu + c (- 0.4 , 0.4 )[chain], dat$ N)<- mean (theta); tau2 <- 0.5 <- individual_loglik (theta, dat$ train, dat$ N)<- matrix (NA_real_ , length (keep), dat$ N)<- tau_draw <- numeric (length (keep))<- numeric (dat$ N); saved <- 0 for (iteration in seq_len (n_iter)) {<- theta + rnorm (dat$ N, 0 , 0.7 )<- individual_loglik (proposal, dat$ train, dat$ N)<- ll_new - ll + dnorm (proposal, mu, sqrt (tau2), log = TRUE ) - dnorm (theta, mu, sqrt (tau2), log = TRUE )<- log (runif (dat$ N)) < logratio<- proposal[accept]; ll[accept] <- ll_new[accept]if (iteration > burn) accepted <- accepted + accept<- 1 / (dat$ N / tau2 + 1 / prior_mu_sd^ 2 )<- variance_mu * (sum (theta) / tau2 + prior_mu / prior_mu_sd^ 2 )<- rnorm (1 , mean_mu, sqrt (variance_mu))<- 1 / rgamma (1 , shape = prior_shape + dat$ N / 2 ,rate = prior_rate + sum ((theta - mu)^ 2 ) / 2 )if (iteration %in% keep) {<- saved + 1 <- theta<- mu; tau_draw[saved] <- sqrt (tau2)list (theta = theta_draw, mu = mu_draw, tau = tau_draw,acceptance = accepted / (n_iter - burn))<- 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))
2本の計算列を、異なる初期値・乱数seedから動かしている。各列16,000回のうち最初の4,000回を捨て、その後4回ごとに保存する。間引きは保存量を減らすためであり、間引いた値が独立になったと保証するものではない。
4. 母集団と個人を、分けて確認する
まず母集団の分布を見よう。\(\mu\) は対数価格係数の平均であり、価格係数の平均そのものではない。後者を知りたいときは、各抽出値から \(\exp(\mu+\tau^2/2)\) を計算する。
<- cbind (mu = fit$ mu, tau = fit$ tau,mean_alpha = exp (fit$ mu + fit$ tau^ 2 / 2 ))<- data.frame (parameter = colnames (hyper_draws),truth = c (dat$ mu_true, dat$ tau_true, 0.12 ),posterior_mean = colMeans (hyper_draws),lower = apply (hyper_draws, 2 , quantile, 0.025 ),upper = apply (hyper_draws, 2 , quantile, 0.975 )):: kable (hyper_table, digits = 3 ,col.names = c ("対象" , "真値" , "事後平均" , "2.5%点" , "97.5%点" ))
mu
mu
-2.30
-2.333
-2.448
-2.228
tau
tau
0.60
0.611
0.506
0.728
mean_alpha
mean_alpha
0.12
0.117
0.106
0.129
母集団の対数平均muと、個人差の標準偏差tauの事後分布。
<- data.frame (mu = fit$ mu, tau = fit$ tau) |> pivot_longer (everything (), names_to = "parameter" , values_to = "value" )ggplot (hyper_long, aes (value)) + geom_density (fill = "#1f77b4" , alpha = 0.4 ) + facet_wrap (~ parameter, scales = "free" ) + labs (x = "母集団パラメータ" , y = "事後密度" )
次に、一人ずつの \(\alpha_i\) を見る。指数関数は非線形なので、\(E[\alpha_i\mid D]\) は \(\exp\{E[\theta_i\mid D]\}\) と一般には一致しない。ここでも各抽出値を変換してから平均する。
<- exp (fit$ theta)<- colMeans (alpha_draws)<- t (apply (alpha_draws, 2 , quantile, c (0.025 , 0.975 )))<- data.frame (id = seq_len (dat$ N), truth = exp (dat$ theta_true),estimate = alpha_mean, lower = alpha_ci[, 1 ], upper = alpha_ci[, 2 ])ggplot (individual_table, aes (truth, estimate)) + geom_abline (slope = 1 , intercept = 0 , color = "grey60" ) + geom_point (alpha = 0.5 ) + labs (x = "個人価格係数の真値" , y = "個人価格係数の事後平均" )ggplot (individual_table[1 : 15 , ], aes (factor (id), estimate)) + geom_errorbar (aes (ymin = lower, ymax = upper), width = 0.2 ) + geom_point (color = "#1f77b4" ) + geom_point (aes (y = truth), shape = 4 , color = "#d62728" ) + labs (x = "最初の15人" , y = "価格係数と95%信用区間" )
母集団の平均を比較的よく推定できても、一人ずつの区間は広い。人数を増やすことと、その人の回答を増やすことは、違う情報を加えている。
Train (2009, Chapter 11)は、個人の係数を条件付き分布から求める方法を整理し、300人のシミュレーションで1人あたりの観測数を1、10、20、50と変えている。個人の観測が増えると、その人の係数についての情報が増える。一方で一人あたりの観測数が固定されたまま人数だけ増えても、各人の係数が必ず真値に収束するわけではない。ここでの「条件付き分布」は頻度論的に推定した母集団分布からも計算でき、すべてがフルベイズの議論ではない。今回の区間が何の不確実性を表すかを考えるときにも役立つ。著者公開版 pp.266–268
5. 計算の診断と、事前分布への感度
事後分布から十分に標本を取れていなければ、平均や区間も安定しない。計算列の時系列を描き、2本が似た範囲を探索しているかを確認しよう。
<- data.frame (draw = rep (seq_len (nrow (fit$ chains[[1 ]]$ theta)), 2 ),chain = factor (fit$ chain), mu = fit$ mu, tau = fit$ tau) |> pivot_longer (c (mu, tau), names_to = "parameter" , values_to = "value" )ggplot (trace, aes (draw, value, color = chain)) + geom_line (alpha = 0.45 , linewidth = 0.3 ) + facet_wrap (~ parameter, scales = "free_y" ) + labs (x = "保存した反復" , y = NULL , color = "chain" )
下では各計算列を前半・後半に分け、列内のばらつきと列間のばらつきを比較する簡単なsplit-\(\hat R\) も計算する。1に近いことは望ましいが、これだけで収束を証明できるわけではない。実務ではrank正規化した診断、実効標本数、Monte Carlo標準誤差なども合わせて確認する。
<- function (fit, tasks) {<- matrix (0 , length (tasks$ id), 3 )for (b in seq_len (nrow (fit$ theta))) {<- answer + choice_prob (tasks$ quality, tasks$ price, exp (fit$ theta[b, tasks$ id]))/ nrow (fit$ theta)<- function (fit, quality, price, n_integration = 80 ) {# 母集団パラメータの不確実性と新しい個人の異質性の両方を積分。 <- qnorm ((seq_len (n_integration) - 0.5 ) / n_integration)<- matrix (0 , nrow (price), 3 )for (b in seq_along (fit$ mu)) {for (z in nodes) {<- answer + choice_prob (quality, price, exp (fit$ mu[b] + fit$ tau[b] * z))/ (length (fit$ mu) * n_integration)<- function (prob, choice) - mean (log (pmax (prob[cbind (seq_along (choice), choice + 1 )], 1e-12 )))<- function (chain1, chain2) {<- floor (length (chain1) / 2 )<- cbind (chain1[1 : n], tail (chain1, n), chain2[1 : n], tail (chain2, n))<- mean (apply (z, 2 , var)); B <- n * var (colMeans (z))sqrt (((n - 1 ) * W / n + B / n) / W)<- data.frame (parameter = c ("mu" , "tau" , "theta1" , "theta150" , "theta300" ),split_Rhat = c (split_rhat (fit$ chains[[1 ]]$ mu, fit$ chains[[2 ]]$ mu),split_rhat (fit$ chains[[1 ]]$ tau, fit$ chains[[2 ]]$ tau),sapply (c (1 , 150 , 300 ), function (i)split_rhat (fit$ chains[[1 ]]$ theta[, i], fit$ chains[[2 ]]$ theta[, i])))):: kable (diagnostics, digits = 3 )
mu
1.003
tau
1.007
theta1
1.000
theta150
1.000
theta300
1.001
:: kable (data.frame (min_acceptance = min (fit$ chains[[1 ]]$ acceptance),max_acceptance = max (fit$ chains[[1 ]]$ acceptance)), digits = 3 )
採用率が高いほどよいとも限らない。候補が現在の値に近すぎると、ほとんど採用されても探索が遅い場合がある。今回の実装は、1次元の個人差を使って更新を理解するためのものである。
事前分布については、\(\tau^2\) のrateを0.5から1へ変えて再推定してみる。この変更では、個人差が大きい範囲に事前の重みが移る。他の設定を固定し、どこまで結論が変わるかを見る。
コード
<- fit_hb (dat, seed = 202610 , prior_rate = 1 )<- data.frame (prior_rate = c (0.5 , 1 ),mu_mean = c (mean (fit$ mu), mean (fit_alternative$ mu)),tau_mean = c (mean (fit$ tau), mean (fit_alternative$ tau)),tau_Rhat = c (split_rhat (fit$ chains[[1 ]]$ tau, fit$ chains[[2 ]]$ tau),split_rhat (fit_alternative$ chains[[1 ]]$ tau, fit_alternative$ chains[[2 ]]$ tau))):: kable (prior_comparison, digits = 3 )
0.5
-2.333
0.611
1.007
1.0
-2.339
0.624
1.005
結果が似ている場合も、試した事前分布の範囲での比較である。別の母集団分布、属性係数の固定、調査回答と実購買の違いまで確認できたことにはならない。
6. 次の選択は、事後分布で平均して予測する
既存顧客 \(i\) の新しい課題での予測は
\[
\Pr(y_i^{\mathrm{new}}=j\mid D,x^{\mathrm{new}},p^{\mathrm{new}})
=\int P_{ij}^{\mathrm{new}}(\theta_i)p(\theta_i\mid D)\,d\theta_i
\approx\frac1B\sum_{b=1}^B P_{ij}^{\mathrm{new}}(\theta_i^{(b)}).
\]
事後から係数を1組取り、選択確率を計算し、それを平均する。 logitは非線形なので、係数の事後平均を1度代入するだけの計算とは違う。
<- predict_existing (fit, dat$ test)<- choice_prob (dat$ test$ quality, dat$ test$ price,exp (colMeans (fit$ theta)[dat$ test$ id])):: kable (data.frame (mean_abs_difference = mean (abs (p_hb - p_plugin)),max_abs_difference = max (abs (p_hb - p_plugin))), digits = 4 ,col.names = c ("事後積分と平均代入の平均絶対差" , "最大絶対差" ))
complete pooling、no poolingとの予測比較
比較対象として、全員に同じ価格係数を当てはめるモデルと、一人ずつ独立に最尤推定するモデルも作る。比較の条件をそろえるため、この2つでも属性係数は同じ既知値に固定する。第7回の、属性係数も推定したモデルと完全に同じ推定ではない。
8問しかない人の最尤推定は、極端な値を取りうる。以下では計算のため価格係数を \([0.01,1]\) に制限する。その範囲に結果が依存することと、端に張り付いた人数も確認する。
<- optimize (function (theta)- sum (individual_loglik (rep (theta, dat$ N), dat$ train, dat$ N)),interval = c (log (0.01 ), log (1 )))$ minimum<- vapply (seq_len (dat$ N), function (i) {<- lapply (dat$ train, function (z)if (is.matrix (z)) z[dat$ train$ id == i, , drop = FALSE ] else z)$ id <- rep (1 , dat$ train$ T)$ choice <- dat$ train$ choice[dat$ train$ id == i]optimize (function (theta) - individual_loglik (theta, tasks_i, 1 ),interval = c (log (0.01 ), log (1 )))$ minimumnumeric (1 ))<- choice_prob (dat$ test$ quality, dat$ test$ price, rep (exp (pooled_theta), length (dat$ test$ id)))<- choice_prob (dat$ test$ quality, dat$ test$ price, exp (individual_theta[dat$ test$ id]))<- data.frame (method = c ("complete pooling" , "no pooling（範囲制限あり）" , "階層ベイズ" ),log_loss = c (log_score (p_pool, dat$ test$ choice), log_score (p_no, dat$ test$ choice),log_score (p_hb, dat$ test$ choice))):: kable (score_table, digits = 4 , col.names = c ("方法" , "検証データのlog loss" ))
complete pooling
1.0277
no pooling（範囲制限あり）
1.0338
階層ベイズ
0.9860
cat ("個別MLEが探索範囲の端付近となった人数:" ,sum (exp (individual_theta) < 0.0101 | exp (individual_theta) > 0.999 ), " \n " )
log lossは、実際に選ばれた選択肢へ割り当てた予測確率の対数を取り、符号を変えて平均した値である。小さいほどよい。推定用の回答をもう一度当てるのではなく、残しておいた8問だけで比較している。
ggplot (score_table, aes (method, log_loss)) + geom_col (fill = "#1f77b4" , width = 0.6 ) + labs (x = NULL , y = "検証log loss（小さい方がよい）" )<- data.frame (predicted = p_hb[, 1 ], observed = as.integer (dat$ test$ choice == 0 )) |> mutate (bin = cut (predicted, breaks = seq (0 , 1 , by = 0.1 ), include.lowest = TRUE )) |> group_by (bin) |> summarise (predicted = mean (predicted), observed = mean (observed), n = n (), .groups = "drop" )ggplot (calibration, aes (predicted, observed, size = n)) + geom_abline (slope = 1 , intercept = 0 , color = "grey60" ) + geom_point (color = "#1f77b4" ) + coord_equal (xlim = c (0 , 1 ), ylim = c (0 , 1 )) + labs (x = "予測した非購入確率" , y = "検証課題での非購入割合" , size = "課題数" )
この標本では階層ベイズのlog lossが最も小さい。ただし、これは今回の生成過程とモデルによる1回の比較である。予測がよかったことと、各人の価格係数を正確に回収できたことも同じではない。
Lenk, DeSarbo, Green and Young (1996)は、179人のパソコン評価データで、推定用16プロファイルと検証用4プロファイルを用いた。推定用の回答を減らし、階層モデルが個人差をどれくらい回収できるかを調べている。元の分析は評価値に対する階層モデルで、今回のCBC logitそのものではない。回答が少なくても情報を共有できることと、個人の係数が完全に分かることを区別する必要がある。別のプロファイルで予測を検証する考え方も、今回と共通している。原論文
7. 履歴がない新規顧客をどう予測するか
新規顧客については、その人の \(\theta\) の事後分布がまだない。そこで、母集団パラメータの各抽出値から新しい人の係数を発生させ、その人の選択確率を計算する。
\[
\Pr(y_{\mathrm{new}}=j\mid D)
=\iint P_j(\theta)\,\phi(\theta;\mu,\tau^2)\,
p(\mu,\tau^2\mid D)\,d\theta\,d\mu\,d\tau^2.
\]
個人の異質性と、母集団分布そのものの不確実性を、両方平均している。母集団平均の係数を持つ「代表的な一人」に置き換える計算ではない。
# Aは3属性すべてあり、Bは長時間バッテリーのみ。価格は千円。 <- matrix (c (sum (dat$ beta_known), dat$ beta_known[1 ]), 1 , 2 )<- matrix (c (12 , 8 ), 1 , 2 )<- predict_new (fit, quality_new, price_new, n_integration = 40 ):: kable (data.frame (choice = c ("買わない" , "A" , "B" ), probability = as.vector (p_new)),digits = 3 , col.names = c ("選択肢" , "新規顧客の予測確率" ))
買わない
0.354
A
0.432
B
0.214
stopifnot (abs (sum (p_new) - 1 ) < 1e-8 , max (abs (rowSums (p_hb) - 1 )) < 1e-8 )
履歴がない人を「価格に敏感な層」と決めつける情報は、今回のモデルにはない。同じ属性情報しかなければ、新規顧客の予測は同じになる。年齢や利用目的から母集団平均を説明するモデルへ拡張することはできるが、その関係もデータと仮定に基づいて学ぶ必要がある。
Rossi, McCulloch and Allenby (1996)は、購買履歴を使った個別クーポンを分析した。論文では、短い履歴を用いた個別化による純収益の増分が一律配布の2.5倍となり、1回の購買履歴だけでも50%大きいというモデルに基づく比較を報告している。これは実験で確認されたクーポン効果ではなく、推定した需要モデルを使った施策比較である。また元のモデルは階層probitで、今回のlogitとは異なる。個人の好みを推定した先で、何を目的に行動を選ぶかが重要になる例であり、次回はこの問いを扱う。原論文の要旨
階層ベイズは推定方法であり、logitとは別の選択モデルの名前ではない。
個人の尤度と母集団分布を組み合わせ、母集団パラメータも更新する。
母集団と個人の不確実性は分けて確認する。短い個人履歴から分かることには限界がある。
予測では、係数の事後分布にわたって選択確率を平均する。
計算の診断、事前分布への感度、未使用データでの予測評価を合わせて確認する。
宿題
推定用の課題を8問から16問に変え、母集団と個人の信用区間がどう変わるか見よう。
事後平均を代入した予測と、事後積分した予測の違いが大きい人を探そう。
「買いやすい顧客にクーポンを配ればよい」が正しいか、値引きで失う利益も含めて考えてみよう。
coding課題は assignment9.qmd を参照。ブラウザで小さな例を試す場合は Playground 9 を使ってほしい。
次回予告
次回は、今回の個人需要の事後分布を使ってクーポンを配る。購入確率を最大にするのではなく、値引きの費用も含めた期待増分利益を比べる。さらに、モデルで選んだ配布方針をRCTで検証するところまで進む。
参考文献
Train, K. E. (2009). Discrete Choice Methods with Simulation , 2nd ed., Chapters 11–12.
Lenk, P. J., DeSarbo, W. S., Green, P. E., and Young, M. R. (1996). “Hierarchical Bayes Conjoint Analysis.” Marketing Science , 15(2), 173–191.
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.
