ゴール
- ベイズ推測における予測/汎化の指標についての全体感を得る。
- looやArvizパッケージを通したstan等のPPL活用において、モデル評価で何をやっているのかの理解を深め、ドキュメントなどの読解力を高める。
やること
- Gelman先生のこのPDFの内容を上から自分なりに適宜整理していく（出版されているものとも若干内容が違うかも）。
- のちに、こちらのブログで整理されている渡辺ベイズの用語や定義との比較？を行う。
前提となる量の定義
Log Predictive Density
Log predictive densityは以下のように定義される。\thetaがパラメータでyが手元のデータなので、対数尤度と同義だと思う。
Log predictive accuracy (or log predictive density for a future data point)
モデルフィットの評価指標として、out-of-sample predictive performanceを知りたい。\tilde{y_i}を未来のデータのある一点の実現値(未来のデータの実現値とは？となってしまうが)として、以下の量をlog predictive accuracy for \tilde{y_i}としている。
\begin{equation}
\log \, \int p(\tilde{y_i} \mid \theta) \, p(\theta \mid y) \, \text{d}\theta
= \log \, \text{E}_{\text{post}}\lbrack p(\tilde{y_i} \mid \theta) \rbrack
= \log \, p_{\text{post}} (\tilde{y_i})
\end{equation}
そこまで厳密な言葉の使い分けがなされている風ではないが、predictive accuracyと言っているときは未知データに対するもので、predictive densityと言っているときは手元のデータに対するものと思っておくといい気がしている。また、上式は、あるデータポイント1つについてのみであることにも注意する。
Expected log predictive density (elpd)
未来のあるデータポイント\tilde{y_i}が分かっていれば、(1)式のlog predictive accuracyを計算できるが、\tilde{y_i}は実際には不明でありその現れ方の揺らぎを考慮をした量を考えたくなる。そこで、(1)式を\tilde{y_i}の分布で期待値を取った量を考える：
\begin{equation}
\text{elpd} = \int \, \, \lbrace \log p_{\text{post}}(\tilde{y_i}) \rbrace f(\tilde{y_i}) \, \text{d}\tilde{y_i} = \text{E}_{\text{f}} \, \lbrack \log \, p_{\text{post}} (\tilde{y_i}) \rbrack
\end{equation}
Expected log pointwise predictive density (elppd)
(2)式はある一つのデータポイントに対する量なので、データ全体のnで和をとった量を考える：
\text{elppd} = \sum_{i=1}^n \text{E}_{\text{f}} \, \lbrack \log \, p_{\text{post}}(\tilde{y_i}) \rbrack
なお、elpdをあるデータ\tilde{y_i}だけでなくデータ全体\, \tilde{y} = (\tilde{y_1}, \tilde{y_2}, \cdots , \tilde{y_n})に対して考えると、
\int_{\tilde{y_1}} \int_{\tilde{y_2}} \cdots \int_{\tilde{y_n}} \log p_{\text{post}}(\tilde{y_1}, \tilde{y_2}, \dots, \tilde{y_n} \mid \text{data}) \, \text{d}\tilde{y_1} \, \text{d}\tilde{y_2} \cdots \text{d}\tilde{y_n}
と書ける(当該PDFには記載なし)。このとき\tilde{y_i}が独立であれば同時分布が積で表せるので、\text{elpd}=\text{elppd}になる。
小まとめ
基本的には、log predictive accuracyを求めたいが、未知のデータは手元にないし現れ方に揺らぎがあるので、現れ方の期待値を取ったelpdが、out-of-sample predictive accuracyを評価するのに適切であるということになる（と思う）。で、これを推定するための指標が後半に色々と出てくる。
Log pointwise predictive density (lppd)
一旦elpdから離れる。
Log predictive accuracyにおいて \log \, p_{\text{post}} (\tilde{y})は分からないので、一旦手元のデータを全て使って得られる量、\log \, p_{\text{post}} (y)を考えてみる。ただし、この量はout-of-sample predictive performanceを見る観点ではバイアスのかかった指標となる。y = (y_1, y_2, \cdots , y_n)が独立のとき
\begin{align*}
\text{lppd} &= \log p_{\text{post}} (y_1, y_2, \cdots, y_n) \\
&= \log \, \prod_{i}^{n} p_{\text{post}} (y_i) \\
&= \sum_{i=1}^n \log \, p_{\text{post}} (y_i) \\
&= \sum_{i=1}^n \log \, \int p(y_i \mid \theta) p(\theta \mid y) \text{d} \theta
\end{align*}
という量を考えられる。これは得られた事後サンプルを用いて、次のように計算できる：
\widehat{\text{lppd}} = \sum_{i=1}^n \log \left ( \frac{1}{S} \sum_{s=1}^S p(y_i \mid \theta^s) \right)
ここでSはMCMCのiterations/drawsの数。それぞれの事後サンプル\theta^sのときのy_iの尤度を足し合わせてSで割ったものの対数を、i=1,\cdots,nまでの個体で足し合わせる。実際の計算はコードを見るのが分かりやすいので、Rがわかる人はこちらのブログなどを参照すると良い。
情報量基準
predictive accuracyを測る指標は情報量基準と呼ばれる。情報量基準は、within-sample biasを補正したものである。補正していない指標としては上述したlppdなどがある。それとは別に、cross-validationによるpredictive accuracyの測定という方法もある。
AIC
log predictive accuracyは事後分布で平均化されていたが、点推定値をプラグインして計算することもできる\text{elpd}_{\hat{\theta}} = \text{E}_{f} \lbrack \log \, p(\tilde{y} \mid \hat{\theta}) \rbrack。ただし、\tilde{y}の分布fが分からないので直接は計算できない。
そこで、lppdのように手元のデータを使って計算する代わりに、評価を手元のデータで行ったことによるバイアスを補正する項k(パラメタ数)を差し引いた量を考える：
\widehat{\text{elpd}}_{\text{AIC}} = \log \, p(y \mid \hat{\theta}_{\text{MLE}}) - k
これを用いてAICは、
\text{AIC} = -2 \cdot \widehat{\text{elpd}}_{\text{AIC}} = -2 \log \, p(y \mid \hat{\theta}_{\text{MLE}}) + 2k
と定義される。
しかし、informative priorsを使う場合や階層モデルの場合、overfittingが抑えられているためこのkによる補正が不適切になる。
DIC
DICもAICと同じように考えていく。
\widehat{\text{elpd}}_{\text{DIC}} = \log \, p(y \mid \hat{\theta}_{\text{Bayes}}) - p_{\text{DIC}}
ここで、\hat{\theta}_{\text{Bayes}}= \text{E} \lbrack \theta \mid y \rbrackで事後平均。
kの代わりに用いられるeffective number of parameters p_{\text{DIC}}は、
\begin{align*}
\text{p}_{\text{DIC}}
&= 2 \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}} ) - \text{E}_{\text{post}} \lbrack \log \, p (y \mid \theta) \rbrack \right) \\
&= 2 \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}} ) - \int \log \, p (y \mid \theta) p(\theta \mid y) \text{d}\theta \right)
\end{align*}
で定義される。これが事前分布からサンプリングした時と事後平均でのDevianceの差を事後分布での期待値として定義されていると考えることができるらしい(上の式では尤度で書かれているので符号関係が逆になっている)。モデルを複雑にするほど事後平均をプラグインした際にデータからの逸脱度は小さくなるため、そのDevianceの差(事前分散でのDevianceから事後平均でのDevianceを引いた量)が大きくなる＝ペナルティが大きくなるという気持ちだと思う。別の見方として、Deviance(尤度を-2倍をしたもの)の事後分布での平均(2倍される括弧内の二項目)から事後分布を代入したDeviance(2倍される括弧内の一項目)を差し引いたものとも言える(WinBUGSの説明などはこれだった)
また、サンプルサイズが大きければkに近づくそうだ。chatGPTに相談して自分で納得した厳密でない証明スケッチとしては次の通りである。 \log p (y \mid \theta)を事後平均など\hat{\theta}_{\text{Bayes}}の周りでテーラー展開して、サンプルサイズが大きく対数尤度のヘシアンが正則などの条件の下で事後分布が正規分布に近似できるとき(こちらのブログに条件の記載あり)は、テーラー展開された二次の項の部分の事後分布での期待値が、自由度=パラメタ数のカイ二乗分布に従うと考えられて、以下のように示せる:
\begin{align*}
\text{p}_{\text{DIC}}
&= 2 \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}} ) - \int \log \, p (y \mid \theta) p(\theta \mid y) \text{d}\theta \right) \\
&\approx 2 \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}} ) - \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}}) - \frac{k}{2}) \right) \right) \\
&= k
\end{align*}
さて話を元に戻す。\widehat{\text{p}_{\text{DIC}}}は手元のデータと事後サンプルでは、
\widehat{\text{p}_{\text{DIC}}} = 2 \left(\log \, p (y \mid \hat{\theta}_{\text{Bayes}} ) - \frac{1}{S} \sum_{s=1}^S \log \, p (y \mid \theta^s) \right)
と求める。
ただし、この\widehat{p_{\text{DIC}}}は負になることがあるため、
p_{\text{DIC'}} = 2 \text{Var}_{\text{post}} \lbrack \log \, p (y \mid \theta) \rbrack
という求め方もある。実際には、N=n \cdot Sとして、
\widehat{p_{\text{DIC'}}} = \frac{1}{N-1} \sum_{s=1}^N \lbrace \log \, p(y_i \mid \theta^s) - \overline{\log \, p(y \mid \theta)} \rbrace ^2
と事後サンプル全体を用いて計算できると思う(この記載はPDFにはない)。そしてDICはいずれかのp_{\text{DIC}}を用いて、
\text{DIC} = -2 \cdot \widehat{\text{elpd}}_{\text{DIC}}= -2 \log \, p(y \mid \hat{\theta}_{\text{Bayes}}) + 2p_{\text{DIC}}
となる。
WAIC
WAICはlppdに補正を加えた情報量基準。補正項の定義は2つあるらしい。一つ目はp_{\text{DIC}_1}と似ていて、
p_{\text{WAIC}_1} = 2 \sum_{i=1}^n \lbrace \log \, \text{E}_{\text{post}} \lbrack p(y_i \mid \theta) \rbrack - \text{E}_{\text{post}} \lbrack \log \, p(y_i \mid \theta) \rbrack \rbrace
である(一つ目の項が事後分布での期待値のlogになっている)。実際には事後サンプルを使って次のように求められる：
\widehat{p_{\text{WAIC}_1}} = 2 \sum_{i=1}^n \lbrace \log \, \frac{1}{S} \sum_{s=1}^S p(y_i \mid \theta^s) - \frac{1}{S} \sum_{s=1}^S \log \, p(y_i \mid \theta^s) \rbrace
もう一つの補正項は、p_{\text{DIC}_2}と同じように分散を用いて（ただしpointwise）、
p_{\text{WAIC}_2} = \sum_{i=1}^n \text{Var}_{\text{post}} \lbrack \log \, p(y_i \mid \theta) \rbrack
と定義される。実際には、手元データと事後を用いて次のように計算する：
\widehat{p_{\text{WAIC}_2}} = \sum_{i=1}^n \left (\frac{1}{S-1} \sum_{s=1}^S \lbrace \log \, p(y_i \mid \theta^s) - \overline{\log \, p(y_i \mid \theta)} \rbrace ^2 \right)
なお、参照しているGelman先生のPDFだと分散を用いて定義されたDICの補正項に関して事後サンプルで計算する式は書かれていない。また、p_{\text{WAIC}_2}の計算ではS-1で割っているが、上で紹介したブログではvarで計算されているのでSで割っていると思われる。ただ、MCMCのiterationsが大きければ無視できる違いなのであまりに気する必要はないはず。
p_{\text{WAIC}_1} または p_{\text{WAIC}_2}を用いて、lppdを補正して得られるelppdのWAIC推定値が得られ、
\widehat{\text{elppd}}_{\text{WAIC}} = \text{lppd} - p_{\text{WAIC}}
ここでは添え字がelpdではなくelppdになっていることに注意。elpdが実際に知りたい指標であるとするならば、データが独立であることを仮定しているということになるのかもしれない。そしてWAICは、
\text{WAIC} = -2 \cdot \widehat{\text{elppd}}_{\text{WAIC}}
と定義される。なお、資料によるとp_{\text{WAIC}_2}の方がおすすめらしい。
LOO-CV lppd
lppdのバイアス補正をLOO-CVで行うことを考えると、
\begin{align*}
\text{lppd}_{\text{LOO-CV}}
&= \sum_{i=1}^n \log \, p_{\text{post}_{(-i)}} (y_i) \\
&= \sum_{i=1}^n \log \, \int p(y_i \mid \theta) \, p(\theta \mid y_{-i}) \, \text{d}\theta
\end{align*}
という量を定義できる。log predictive accuracyのところでp_{\text{post}} (y_i)は\text{E}_{\text{post}}\lbrack p(y_i \mid \theta) \rbrack = \log \, \int p(y_i \mid \theta) \, p(\theta \mid y) \, \text{d}\thetaをまとめたものとして記載していた。LOO-CVはあるデータポイントiを除いたデータでデータポイントiを得られた未知のデータにfitを行うので、yをy_iを除いたy_{-i}=(y_1, \cdots, y_{i-1}, y_{i+1}, \cdots, y_{n})に基づいた事後分布になっている。
手元で計算する場合には、
\widehat{\text{lppd}_{\text{LOO-CV}}} = \sum_{i=1}^n \log \, \left( \frac{1}{S} \sum_{s=1}^S p(y_i \mid \theta_{-i}^s) \right)
となる。ここで、\theta_{-i}はp(\theta \mid y_{-i})による平均化であることを表すために表記を変更したものである。
バイアス補正項(effectie number of parameters)は、
p_{\text{LOO-CV}} = \text{lppd} - \text{lppd}_{\text{LOO-CV}}
で得られる。
渡辺ベイズとの用語比較
ここからは私の整理である。
予測分布
lppdにおけるp_{\text{post}}(y_i) = \small\int p(y_i \mid \theta) p(\theta \mid y) \text{d} \thetaの部分を「予測分布 p^*(y)」と呼んでいる。須山先生のベイズ本にある記載（あるいはご本人のこちらのブログ）でも、予測分布(predictive distribution)と説明されている。
個人的によく混乱するのが、\small\int p(y_i \mid \theta) p(\theta \mid y) \text{d} \thetaと記載があったときに、今手元にあるデータの尤度をパラメータの事後分布で平均化したものと考えるのか、事後分布からサンプルされたパラメータを用いて確率モデルから新たに生み出されるy_iを考えているのか、の区別である。専門家としては区別しなくても良いのかもしれないが、用語が色々あることも相まって素人的にはよく分からなくなる。
Gelman先生の記事中はyを観測されたデータと捉えているので、\small\int p(y_i \mid \theta) p(\theta \mid y) \text{d} \thetaは尤度を平均化したものであると分かる。\small\int p(\tilde{y_i} \mid \theta) p(\theta \mid y) \text{d} \thetaとすると尤度でないことは分かるのだが、\tilde{y_i}が未知のデータと捉えるのか、その確率モデルからこれからサンプルされるデータと捉えるのかによっても、predictive accuracyの話をしているのか、predictive distributionの話をしているのかが分かりづらいように思える。Stanユーザーガイドでは、y^\text{rep}としてyのreplicationsを取り出すことを明示的にしているようだ（この場合、事後分布を用いているので"posterior" predictive distributionであることに注意。)
改めて渡辺先生の本の記載を見て、本記事での記法に添いつつ予測分布の定義を書くと、
p^*(y) = \int p(y \mid \theta) p(\theta \mid Y^n) \text{d} \theta
となる。Y^n=(Y_1, Y_2, \cdots, Y_n)がサンプルの確率変数で、y^nはその実現値と説明されている（p.2)ので、尤度ぽい見方をしているのだろうか。
改めて須山先生のブログの記載を見て、本記事での記法に添いつつ予測分布の定義を書くと、
p(y \mid D) = \int p(y \mid \theta, D) p(\theta \mid D) \text{d}\theta
となる。観測されたデータをDとしてまとめてあり、yは新しいデータとされている。その下に「まだ入ってきていない新データに関する確率分布（平均とか分散とか）を計算しようというわけです。」と書かれているので、どちらかというとposterior predictive distributionのことを意識していると感じる。
何か解決したわけではないが、個人的なメモはここで終わり。
サンプルの現れ方の平均
改めてこの用語を見ると、なんじゃそれはというふうに思ってしまっていたのだが、elpdの考え方を整理しているときに出てきたことがまさにそうなのだと思う（以下、elpdの記載の一部を再掲）：
elpdをあるデータ\tilde{y_i}だけでなくデータ全体\, \tilde{y} = (\tilde{y_1}, \tilde{y_2}, \cdots , \tilde{y_n})に対して考えると、
\int_{\tilde{y_1}} \int_{\tilde{y_2}} \cdots \int_{\tilde{y_n}} \log p_{\text{post}}(\tilde{y_1}, \tilde{y_2}, \dots, \tilde{y_n} \mid \text{data}) \, \text{d}\tilde{y_1} \, \text{d}\tilde{y_2} \cdots \text{d}\tilde{y_n}
経験誤差(Training error)
経験誤差は次のように定義されている：
\text{T}_n = \frac{1}{n} \sum_{i=1}^n - \log \, p^*(y_i)
上で、p^*(y_i)がp_{\text{post}}(y_i)と同じ量であることを見たので、\text{lppd}で\text{T}_nを表すことを考える。\text{lppd} = \sum_{i=1}^n \log \, p_{\text{post}} (y_i)であるから、
\begin{align*}
\text{T}_n
&= - \frac{1}{n} \sum_{i=1}^n \log \, p^*(y_i) \\
&= - \frac{1}{n} \sum_{i=1}^n \log \, p_{\text{post}} (y_i) \\
&= - \frac{1}{n} \text{lppd}
\end{align*}
となるはずである。
WAIC（再び）
渡辺本では、WAICは逆温度\betaと汎関数分散\text{V}_nを用いて、
\text{W}_n = \text{T}_n + \frac{\beta \text{V}_n}{n}
と定義されている。
汎関数分散\text{V}_nは本記事の記法で書くと：
\text{V}_n = \sum_{i=1}^n \lbrace \text{E}_{\text{post}} \lbrack (log \, p(y_i \mid \theta))^2 \rbrack - (\text{E}_{\text{post}} \lbrack log \, p(y_i \mid \theta) \rbrack)^2 \rbrace
と定義されている。これは上で定義を確認したGelman先生のp_{\text{WAIC}_2}と同じ定義である。気持ちとしては、対数尤度の事後分布で平均化したもののばらつきをiごとに見ているので、この値が大きいほど、モデルp(y_i \mid \theta)のデータへの適合度がばらつく = バイアス・バリアンス・トレードオフ的には過適合しているということになる。より柔軟なモデルほどこの値は大きくなるので、柔軟=モデルのパラメータが多い、ということになるのだと思う。
逆温度\beta=1の下で、Gelman先生の方の定義の\text{WAIC}を渡辺先生の定義W_nに書き替えていくと、
\begin{align*}
\text{WAIC}
&= -2 \cdot \widehat{\text{elppd}}_{\text{WAIC}} \\
&= -2 (\text{lppd} - p_{\text{WAIC}}) \\
&= -2(-n\text{T}_n - p_{\text{WAIC}}) \\
&= 2(n\text{T}_n + p_{\text{WAIC}}) \\
&= 2 n(\text{T}_n + \frac{p_{\text{WAIC}}}{n}) \\
&= 2 n \cdot \text{W}_n
\end{align*}
最後の等式はp_{\text{WAIC}_2}=\text{V}_nならば成り立つことは上で見ていてる(p_{\text{WAIC}_1}でもnが大きくなれば成り立つのだろうと思われる)。参照しているPDFでは「渡辺先生のオリジナルのWAICの定義では、2がかけられておらずnで割られている」とあることが最低限整理できた。Gelman先生の方では、AICなどの他のdevianceベースの指標と揃えられるようにリスケールしたらしい。
後記
- 色々参照して、解像度は高まった気がする。
- elpdの推定量が色々出てくるところからハットをつけるべきなのかどうか（PDFではついている）が悩ましいなと思い始めた（手元のデータで得られたMLEの話をしているのか、MLEという推定量の話をしているのかで、ハットをつけなくて良い後者のような考え方でもいいのでは）。
- 渡辺ベイズでは汎化誤差や自由エネルギーを真の分布の近さの指標の基礎として、情報量基準がそれを近似するかどうかで整理しているので枠組みとして良いなと思った。証明とかも含めていつかやりたい...(真のモデルはなという前提の元でも成り立つ議論なのかはよく分からないけど)。
Discussion