尤度は「このパラメータなら、このデータがどれくらい出やすいか」を表す。事前分布と組み合わせると、データを見た後の事後分布になる。
事後分布を平均や区間で要約することと、次の観測を予測することを区別する。
個人ごとのデータが少ないとき、全員を同じとする方法と、一人ずつ別々に推定する方法には、それぞれ難しさがある。
階層モデルは、個人差を残しながら情報を共有する。経験ベイズとフルベイズでは、母集団パラメータの不確実性の扱いが違う。
前回とのつながり・今回のゴール
前回は、商品の属性や価格を変えた選択課題を作り、回答から需要を推定した。ただし、一人の回答者に何百問も答えてもらうのは難しい。8問の回答からその人の好みを知りたい、と言われると、全員をまとめて推定したモデルだけでは物足りない。一方で、8問だけを使って一人ずつ推定すると、たまたまの回答に強く左右される。
そこで今回は、個人ごとの推定に、他の人のデータをどう使うか を考える。そのためにベイズ推定を導入する。いきなり選択モデルと一緒に扱うと、何が新しいのか分かりにくいので、まずは正規分布の簡単なモデルで仕組みを確認しよう。次回、この仕組みを前回の選択データに戻す。
1. 尤度の復習：データを固定し、パラメータを動かす
計量経済学IのOLSでは、回帰係数を変えたときに残差平方和がどう変わるかを見た。第2回の最尤法では、観測された選択が出る確率を、係数の関数として比較した。この「データを固定して係数を動かす」という見方が、尤度だった。
今回は、顧客の評価を標準化したスコアで考える。顧客の平均的な評価を \(\theta\) 、観測した平均スコアを \(y\) として、
\[
y\mid\theta\sim N(\theta,s^2)
\]
とする。\(s\) は平均スコアの標準誤差で、ここでは既知とする。例えば \(y=1.5,s=0.8\) が観測されたら、尤度は
\[
L(\theta;y)=\frac{1}{\sqrt{2\pi s^2}}
\exp\left[-\frac{(y-\theta)^2}{2s^2}\right]
\]
である。\(\theta=1.5\) の近くで大きくなる。最尤推定値は \(y\) そのものだ。ただし、尤度はパラメータについての確率分布として定義したものではない。この区別をつけてから、次に進もう。
2. 事前分布と事後分布
データを見る前に、\(\theta\) についてどの範囲を考えていたかを、事前分布 で表す。例えば
\[
\theta\sim N(m_0,v_0)
\]
とする。\(m_0\) は中心、\(v_0\) は事前の不確実性である。過去の調査や別のデータを参考に決めることもあれば、意味のある広い範囲を覆うように決めることもある。どちらの場合も、観測単位を変えれば妥当な値は変わる。
観測した \(y\) を使って更新した分布が事後分布 である。ベイズの公式は
\[
p(\theta\mid y)=\frac{L(\theta;y)p(\theta)}
{\int L(a;y)p(a)\,da}
\quad\Longleftrightarrow\quad
p(\theta\mid y)\propto L(\theta;y)p(\theta).
\]
分子は尤度と事前分布の積、分母は全体の面積を1にするための定数である。ここまでで使っているのは、条件付き確率の公式である。
正規分布どうしを組み合わせると、事後分布も正規になる。
\[
\theta\mid y\sim N(m_1,v_1),\qquad
v_1=\left(\frac1{s^2}+\frac1{v_0}\right)^{-1},\qquad
m_1=v_1\left(\frac{y}{s^2}+\frac{m_0}{v_0}\right).
\]
分散の逆数を精度 と呼ぶ。事後平均は、データと事前の中心を精度で重みづけした平均になっている。
<- 1.5 ; s_one <- 0.8 <- 0 ; v0 <- 1 <- 1 / (1 / s_one^ 2 + 1 / v0)<- v1 * (y_one / s_one^ 2 + m0 / v0)<- seq (- 3 , 4 , length.out = 400 )<- data.frame (theta = rep (grid, 3 ),density = c (dnorm (grid, m0, sqrt (v0)), dnorm (grid, y_one, s_one),dnorm (grid, m1, sqrt (v1))),distribution = rep (c ("事前" , "尤度の形（面積を1に調整）" , "事後" ), each = length (grid)))ggplot (update_plot, aes (theta, density, color = distribution)) + geom_line (linewidth = 0.9 ) + labs (x = expression (theta), y = "密度" , color = NULL ) + theme (legend.position = "bottom" ):: kable (data.frame (posterior_mean = m1, posterior_sd = sqrt (v1)), digits = 3 ,col.names = c ("事後平均" , "事後標準偏差" ))
事前分布と、観測によって更新した事後分布。尤度は形を比較するため面積を1に調整している。
図では、尤度の形を比較しやすくするために面積を1に調整している。一般のモデルで、尤度をそのまま事後分布と呼んでよいわけではない。
事前分布を変えると、どれくらい変わるか
事前分布を入れた以上、その選び方が結果に影響する。影響があること自体を隠すより、妥当だと考えられる別の事前分布でも結果を比べよう。
<- c (0.3 , 1 , 3 )<- do.call (rbind, lapply (prior_sd, function (sd0) {<- 1 / (1 / s_one^ 2 + 1 / sd0^ 2 )<- vv * y_one / s_one^ 2 data.frame (theta = grid, density = dnorm (grid, mm, sqrt (vv)),prior = paste0 ("事前SD = " , sd0))ggplot (sensitivity, aes (theta, density, color = prior)) + geom_line (linewidth = 0.9 ) + labs (x = expression (theta), y = "事後密度" , color = NULL )
事前分布を0の近くに強く集中させると、事後平均も0に寄る。データが少ない状況では、極端に狭い事前分布の影響をデータだけで打ち消すことは難しい。
登録したばかりの顧客が、最初の商品を高く評価したとする。その1回だけで「この人は何でも高く評価する」と決めてしまうのは早い。過去の顧客では評価がどれくらいばらついていたかと、今の顧客について何件観測したかを、一緒に考えたい。この後の階層モデルでは、過去の顧客についての情報もデータから学ぶ。
3. パラメータの区間と、次の観測の区間
事後分布の中央95%を取った区間を、95%信用区間 と呼ぶ。置いたモデルと事前分布のもとで、データを条件にした事後確率が95%になる区間である。計量Iで学んだ信頼区間の「同じ手続きを繰り返したとき、真値を含む割合」という定義とは異なる。
次に知りたいのが、同じ顧客の次のスコア \(y^{\mathrm{new}}\) なら、\(\theta\) の不確実性に、新しい観測のノイズも加える。
\[
p(y^{\mathrm{new}}\mid y)
=\int p(y^{\mathrm{new}}\mid\theta)p(\theta\mid y)\,d\theta.
\]
同じ測定分散 \(s^2\) なら、この例では \(y^{\mathrm{new}}\mid y\sim N(m_1,v_1+s^2)\) となる。事後平均だけを代入して \(N(m_1,s^2)\) とするのでは、パラメータの不確実性を落としてしまう。
<- data.frame (value = rep (grid, 2 ),density = c (dnorm (grid, m1, sqrt (v1)), dnorm (grid, m1, sqrt (v1 + s_one^ 2 ))),target = rep (c ("平均的評価 theta" , "次の観測 y_new" ), each = length (grid)))ggplot (prediction_plot, aes (value, density, color = target)) + geom_line (linewidth = 0.9 ) + labs (x = "評価スコア" , y = "密度" , color = NULL )
4. 全員同じ、一人ずつ、部分的に共有する
一人の更新が分かったので、顧客が何人もいる場合を考えよう。顧客 \(i\) の平均評価を \(\theta_i\) 、その人の観測平均を \(y_i\) とする。
\[
y_i\mid\theta_i\sim N(\theta_i,s_i^2),\qquad
\theta_i\mid\mu,\tau^2\sim N(\mu,\tau^2).
\]
\(\tau^2\) は顧客間の本当の違い、\(s_i^2\) はその人の平均を測るときのノイズである。1回の観測の分散が \(\sigma^2\) 、観測回数が \(n_i\) なら、\(s_i^2=\sigma^2/n_i\) となる。
complete pooling
全員が同じパラメータ
全員の観測をまとめる
no pooling
全員を別々に推定
他の人の情報を使わない
partial pooling
共通の分布から個人差が生じる
個人の情報と母集団の情報を合わせる
\(\mu,\tau^2\) を固定したとき、先ほどの正規分布の更新をそのまま使える。
\[
E[\theta_i\mid y_i,\mu,\tau^2]
=B_i\mu+(1-B_i)y_i,\qquad
B_i=\frac{s_i^2}{s_i^2+\tau^2}.
\]
\(B_i\) が大きいほど母集団平均へ寄せる。この動きを縮小（shrinkage） と呼ぶ。観測が少なく \(s_i^2\) が大きい人ほど強く寄せる一方、本人のデータが豊富なら \(y_i\) を重く使う。
<- expand.grid (n = 1 : 60 , tau = c (0.3 , 0.7 , 1.5 ))$ B <- (4 / weights$ n) / (4 / weights$ n + weights$ tau^ 2 )ggplot (weights, aes (n, B, color = factor (tau))) + geom_line (linewidth = 0.9 ) + labs (x = "個人の観測回数" , y = "母集団平均にかける重み B" , color = "個人差のSD" )
個人差が大きいと分かっている場合は、同じデータ量でもあまり寄せない。全員を平均へ寄せるほどよい、という方法ではない。
EfronとMorrisの野球の例では、18人の選手について、シーズン初めの45打数の成績から残りの期間の打率を予測する。早い時期の打率には本当の能力差と偶然の両方が含まれる。Efronによる紹介では、Clementeの最初の打率は0.400、残りの期間では0.346だった。全員を別々に扱う推定と、選手間で情報を共有する推定を比べる例である。ただし、縮小で全選手の誤差が必ず改善する、という主張ではない。著者による解説とデータ
5. 経験ベイズ：母集団の分布をデータから推定する
現実には \(\mu,\tau^2\) も分からない。そこで、まず全員のデータからこの2つを推定し、その推定値を固定して個人の分布を計算する方法が経験ベイズ である。
このモデルでは \(y_i\sim N(\mu,\tau^2+s_i^2)\) なので、この周辺尤度を最大にする \(\mu,\tau\) を求めればよい。ここで「周辺」と呼んでいるのは、観測できない \(\theta_i\) を積分して取り除いたためである。
60人のデータを作る。評価スコアを正規分布で近似する練習であり、購入の有無をこの正規モデルで推定しているわけではない。測定分散 \(\sigma^2=4\) は既知とする。
set.seed (202608 )<- 60 <- sample (c (3 , 8 , 20 , 50 ), N, replace = TRUE )<- 1 ; tau_true <- 0.7 ; sigma <- 2 <- sigma^ 2 / n_i<- rnorm (N, mu_true, tau_true)<- rnorm (N, truth, sqrt (s2))<- rnorm (N, truth, sqrt (s2))# 経験ベイズでは、周辺尤度で母集団パラメータを推定する。 <- function (par) {- sum (dnorm (y, par[1 ], sqrt (exp (2 * par[2 ]) + s2), log = TRUE ))<- optim (c (mean (y), log (0.5 )), neg_marginal, method = "L-BFGS-B" ,lower = c (- 5 , - 5 ), upper = c (5 , 2 ))<- eb_fit$ par[1 ]; tau_eb <- exp (eb_fit$ par[2 ])<- s2 / (s2 + tau_eb^ 2 )<- B * mu_eb + (1 - B) * y<- rep (weighted.mean (y, 1 / s2), N):: kable (data.frame (parameter = c ("母集団平均" , "個人差のSD" ),truth = c (mu_true, tau_true), empirical_bayes = c (mu_eb, tau_eb)), digits = 3 ,col.names = c ("パラメータ" , "真値" , "経験ベイズ" ))
母集団平均
1.0
0.977
個人差のSD
0.7
0.667
<- data.frame (truth = rep (truth, 3 ), estimate = c (complete, y, eb_mean),method = rep (c ("complete pooling" , "no pooling" , "partial pooling：EB" ), each = N))ggplot (pooling, aes (truth, estimate)) + geom_abline (slope = 1 , intercept = 0 , color = "grey60" ) + geom_point (alpha = 0.7 , color = "#1f77b4" ) + facet_wrap (~ method) + labs (x = "個人の真値" , y = "推定値" )<- data.frame (method = c ("complete" , "no pooling" , "EB partial" ),parameter_mse = c (mean ((complete - truth)^ 2 ), mean ((y - truth)^ 2 ), mean ((eb_mean - truth)^ 2 )),future_mse = c (mean ((complete - y_test)^ 2 ), mean ((y - y_test)^ 2 ), mean ((eb_mean - y_test)^ 2 ))):: kable (comparison, digits = 3 ,col.names = c ("方法" , "真の個人パラメータに対するMSE" , "別の観測平均に対するMSE" ))
complete
0.438
1.075
no pooling
0.490
1.026
EB partial
0.191
0.813
真値との比較。partial poolingは個人差を残しながら、ノイズによる散らばりを抑える。
真値に対する誤差はシミュレーションだから計算できる。実データでは、推定に使わなかった後の観測に対する予測誤差を確認する。観測の予測には測定ノイズも入るため、2つのMSEを同じものと考えないようにしよう。
6. フルベイズ：母集団パラメータにも事前分布を置く
経験ベイズでは \(\hat\mu,\hat\tau^2\) を固定した。この2つも有限のデータから推定しているので、本当は不確実である。フルベイズ では、その不確実性も事後分布に含める。
ここでは
\[
\mu\sim N(0,2^2),\qquad \tau^2\sim IG(2,1)
\]
とする。\(IG(a,b)\) は密度が \((\tau^2)^{-a-1}\exp(-b/\tau^2)\) に比例する逆ガンマ分布である。以下では、形状 \(a\) 、rate \(b\) のガンマ乱数の逆数として生成する。この事前分布を「何も仮定していない」と解釈してはいけない。
全員のパラメータと母集団パラメータの事後分布は
\[
p(\theta_1,\ldots,\theta_N,\mu,\tau^2\mid y)
\propto \left\{\prod_i p(y_i\mid\theta_i)p(\theta_i\mid\mu,\tau^2)\right\}
p(\mu)p(\tau^2).
\]
この分布から値を多数取り出せれば、個人ごとの平均や区間も、\(\mu,\tau\) の区間も計算できる。そのための方法として、今回はGibbsサンプリング を使う。他のパラメータを仮に固定し、\(\theta\) 、\(\mu\) 、\(\tau^2\) を順番に条件付き分布から引き直す方法である。点推定値へ収束させる計算ではなく、事後分布の中を動き回る計算だ。
\(\theta_i\) の更新には、すでに導いた正規分布を使う。\(\mu\) も正規分布、\(\tau^2\) も逆ガンマ分布になる。例えば \(\tau^2\) の条件付き分布は、形状が \(2+N/2\) 、rateが \(1+\sum_i(\theta_i-\mu)^2/2\) である。
<- function (y, s2, seed = 202608 , mu0 = 0 , mu_sd = 2 ,shape0 = 2 , rate0 = 1 , iterations = 5000 , burn = 1000 ) {set.seed (seed)<- length (y); theta <- y; mu <- mean (y); tau2 <- 1 <- matrix (0 , iterations - burn, N + 2 )for (iteration in seq_len (iterations)) {<- 1 / (1 / s2 + 1 / tau2)<- variance * (y / s2 + mu / tau2)<- rnorm (N, mean_theta, sqrt (variance))<- 1 / (N / tau2 + 1 / mu_sd^ 2 )<- rnorm (1 , variance_mu * (sum (theta) / tau2 + mu0 / mu_sd^ 2 ),sqrt (variance_mu))<- 1 / rgamma (1 , shape0 + N / 2 , rate0 + sum ((theta - mu)^ 2 ) / 2 )if (iteration > burn) draws[iteration - burn, ] <- c (mu, sqrt (tau2), theta)colnames (draws) <- c ("mu" , "tau" , paste0 ("theta" , seq_len (N)))<- normal_hb (y, s2, seed = 202608 )<- normal_hb (y, s2, seed = 202609 )<- rbind (chain1, chain2)<- posterior[, - (1 : 2 )]<- colMeans (post_theta)<- t (apply (post_theta, 2 , quantile, c (0.025 , 0.975 )))
<- data.frame (parameter = c ("mu" , "tau" ), truth = c (mu_true, tau_true),mean = colMeans (posterior[, 1 : 2 ]),lower = apply (posterior[, 1 : 2 ], 2 , quantile, 0.025 ),upper = apply (posterior[, 1 : 2 ], 2 , quantile, 0.975 )):: kable (posterior_summary, digits = 3 ,col.names = c ("パラメータ" , "真値" , "事後平均" , "2.5%点" , "97.5%点" ))
mu
mu
1.0
0.978
0.750
1.209
tau
tau
0.7
0.689
0.517
0.901
観測が少ない15人の事後平均と95%信用区間。赤い印はシミュレーションの真値。
<- order (n_i)[1 : 15 ]<- data.frame (id = factor (show_i, levels = show_i),mean = full_mean[show_i], lower = intervals[show_i, 1 ], upper = intervals[show_i, 2 ],truth = truth[show_i])ggplot (interval_plot, aes (id, mean)) + 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%信用区間" )
赤い印は真値で、実データでは見えない。点推定では順位がついていても、区間は大きく重なることがある。「推定値が一番高い人」と「確実に一番高い人」は違う。
計算が動いたことと、十分な標本が得られたこと
乱数を出せばよいわけではない。異なる乱数列で動かした2本の計算が、似た範囲を探索しているか確認する。最初の1,000回を捨て、残りを使っているが、この回数で常に十分になるわけではない。
<- data.frame (iteration = rep (seq_len (nrow (chain1)), 2 ),mu = c (chain1[, "mu" ], chain2[, "mu" ]), chain = rep (c ("1" , "2" ), each = nrow (chain1)))ggplot (trace_plot, aes (iteration, mu, color = chain)) + geom_line (alpha = 0.45 , linewidth = 0.3 ) + labs (x = "保存した反復" , y = expression (mu), color = "chain" )<- data.frame (chain = c ("1" , "2" ),mu_mean = c (mean (chain1[, "mu" ]), mean (chain2[, "mu" ])),tau_mean = c (mean (chain1[, "tau" ]), mean (chain2[, "tau" ]))):: kable (chain_comparison, digits = 3 )
1
0.979
0.690
2
0.977
0.688
異なる乱数seedから得た2列の母集団平均の抽出値。
連続した抽出値は相関する。保存した数が、そのまま独立な標本数になるわけではない。次回は、複数の計算列の比較を数値で調べる診断も使う。
Rubin (1981)の例では、8校で行われたSAT対策指導の無作為化実験をまとめる。学校ごとの効果推定値には、その学校の真の効果と推定誤差の両方が含まれる。階層モデルを使うと、各校を完全に別々に扱う方法と、すべて同じ効果とする方法の間で情報を共有できる。Rubin (1984)は、この例で同じ8校の実験を再び行ったときのデータを事後予測する考え方も説明している。因果効果を識別する実験の役割と、推定誤差を扱う階層モデルの役割は別である。原著者の解説 §5.5
7. 事後分布を使って予測する
既存の顧客なら、その人の \(\theta_i\) の事後抽出値を使う。新しい顧客なら、各反復の \(\mu,\tau\) から新しい \(\theta_{\mathrm{new}}\) を発生させる。どちらも、最後に観測のノイズを加える。
set.seed (2026081 )<- which.min (n_i)<- rnorm (nrow (posterior), post_theta[, i], sigma)<- rnorm (nrow (posterior), posterior[, "mu" ], posterior[, "tau" ])<- rnorm (nrow (posterior), theta_new, sigma)<- data.frame (value = c (future_existing, future_new),customer = rep (c ("既存顧客" , "新規顧客" ), each = nrow (posterior)))ggplot (predictive, aes (value, color = customer)) + geom_density (linewidth = 0.9 ) + labs (x = "次の1回の評価スコア" , y = "事後予測密度" , color = NULL )
新規顧客については、個人データがない分、母集団の情報に頼る。既存顧客の予測分布には、その人自身の履歴が反映される。ただし、次の観測に固有のノイズは、履歴を増やしてもなくならない。
同じ母集団から来た人だと扱うことは、モデルの仮定である。法人顧客と個人顧客、初回購入と長年の利用者などで体系的な違いがあるなら、母集団平均を属性で説明したり、分布を分けたりする必要がある。何でも一つの分布にまとめれば、意味のある差まで縮めることになる。
Lenk, DeSarbo, Green and Young (1996)は、パソコンの13属性を使ったコンジョイント調査を分析した。分析対象は179人で、各人に16個の推定用プロファイルと4個の検証用プロファイルがある。回答の一部を取り除き、階層ベイズが短い質問票でどこまで個人差を推定できるかを調べた。元の分析は評価値に対するモデルであり、次回のlogitと同じモデルではない。少ない個人データと母集団の情報を組み合わせ、別の回答で検証するという問題は共通している。原論文 pp.179–180
尤度と事前分布を組み合わせて、事後分布を求める。
事後分布の幅と、新しい観測の予測分布の幅を区別する。
partial poolingは、個人差と測定ノイズを区別しながら情報を共有する。
経験ベイズは母集団パラメータを推定して固定し、フルベイズはその不確実性も積分する。
推定に使っていないデータで予測を評価し、事前分布や共有する母集団の仮定も確認する。
宿題
観測回数を全員50回にしたら、complete/no/partial poolingの差はどうなるか。
本当の個人差を \(\tau=0.1\) と \(\tau=1.5\) に変え、どの推定が有利になるか比べよう。
「事後平均の予測」と「次の観測の予測」の違いを、式を使わずに説明してみよう。
coding課題は assignment8.qmd を参照。ブラウザで試す場合は Playground 8 を使ってほしい。
次回予告
次回は、第7回のイヤホンの選択課題に戻る。今回は連続的な評価スコアを観測したが、次回はどの商品を選んだか、あるいは何も買わなかったかを観測し、選択確率に入る個人の係数を推定する。個人ごとの選択履歴から事後分布を求め、次に何を選ぶかを予測する。
参考文献
Efron, B. (2010). “The Future of Indirect Evidence.” Statistical Science , 25(2), 145–157.
Rubin, D. B. (1981). “Estimation in Parallel Randomized Experiments.” Journal of Educational Statistics , 6(4), 377–401.
Rubin, D. B. (1984). “Bayesianly Justifiable and Relevant Frequency Calculations for the Applied Statistician.” Annals of Statistics , 12(4), 1151–1172.
Lenk, P. J., DeSarbo, W. S., Green, P. E., and Young, M. R. (1996). “Hierarchical Bayes Conjoint Analysis.” Marketing Science , 15(2), 173–191.
Train, K. E. (2009). Discrete Choice Methods with Simulation , 2nd ed., Chapter 12.
