統計的決定理論から見るベイズ推定

# はじめに
千葉大学・株式会社Nospareの[川久保](https://sites.google.com/site/ykawakubostat/)です．今回は，ベイズ推定という手法を，統計的決定理論（statistical decision theory）という枠組みで説明したいと思います．


# パラメータの点推定
観測データ$X_1,\dots,X_n$は，独立同一にパラメータ$\theta$を持つ確率分布（密度）$f(\cdot \mid \theta)$に従うとします．このとき，$X = (X_1,\dots,X_n)$の同時密度関数は，
```math
f_n(x \mid \theta) = \prod_{i=1}^n f(x_i \mid \theta)
```
となります．

### 損失とリスク
統計的決定理論では，パラメータ$\theta$を推定量$\delta(X)$で点推定したときの損失$L(\theta,\delta(X))$を定義します．最も代表的な損失関数（loss funciton）は二乗損失で，
```math
L(\theta, d) = (d - \theta)^2
```
です．損失$L(\theta,\delta(X))$は標本$X$の実現値に依存するので，推定量の良さを測るには，$X$の従う確率分布について期待値をとることが考えられます．
```math
R(\theta,\delta) := E[L(\theta,\delta(X))] = \int L(\theta, \delta(x)) f_n(x \mid \theta) dx
```
これを，リスクと呼びます．特に，二乗損失を用いた場合のリスクは，
```math
E[(\delta(X) - \theta)^2]
```
となり，平均二乗誤差（mean squared error, MSE）と呼ばれています．

どのようなパラメータ$\theta$の値に対しても一様に，リスクが最小になるような推定量を構成することは不可能です．それは，$\delta_0(X) = 0$という極端な推定量を考えることで説明されます．この推定量は，観測データ$X$によらず常に定数0として$\theta$を推定するという，明らかに不適切な推定量ですが，もし真のパラメータ値が$\theta = 0$であった場合は，$\delta_0(X)$のリスクは$R(0,\delta_0) = 0$になります．つまり，$\delta_0$を優越する（すべてのパラメータ値$\theta$においてリスクが小さくなる）推定量は存在しないのです．この議論は，任意のパラメータ空間上の点について成り立つため，一様に最良となる推定量は存在しないのです．

そこで，推定量に不偏性を課すなど推定量のクラスを制約して，そのクラスの中で最良の推定量を構成したり（例: 一様最小分散不偏推定量），最もリスクの高いパラメータ値のリスクを最小にするというミニマックス性を議論したりします．


### ベイズリスクとベイズ推定量
また，リスク$R(\theta,\delta)$の，異なる$\theta$の値にわたる「平均」を最小化しようとする考え方もあります．パラメータ$\theta$の確率密度$\pi(\theta)$を仮定し，$\pi(\theta)$について$R(\theta,\delta)$の期待値をとった
```math
\tag{1}
r(\pi,\delta) = \int R(\theta,\delta) \pi(\theta) d\theta
```
のことを，**ベイズリスク**と呼びます．そして，パラメータ$\theta$に仮定した確率分布はまさに，事前分布です．ベイズの定理から事後密度$\pi(\theta \mid x)$は，
```math
\pi(\theta \mid x) = \frac{f_n(x \mid \theta) \pi(\theta)}{\int f_n(x \mid \theta) \pi(\theta) d\theta}
```
であり，右辺の分母$f_\pi(x) = \int f_n(x \mid \theta)\pi(\theta)d\theta$は周辺尤度と呼ばれています．$f_n(x \mid \theta)\pi(\theta) = \pi(\theta \mid x) f_\pi(x)$を用いると，(1)式は，
```math
\tag{2}
\begin{split}
r(\pi,\delta) &= \int \left\{ \int L(\theta,\delta(x)) f_n(x \mid \theta) dx \right\} \pi(\theta) d\theta \\
&= \int \left\{ \int L(\theta, \delta(x)) \pi(\theta \mid x) d\theta \right\} f_\pi(x) dx
\end{split}
```
と変形できます．(2)式の最右辺における
```math
\int L(\theta, \delta(x)) \pi(\theta \mid x) d\theta = E[L(\theta, \delta(X)) \mid X]
```
は，事後期待損失（posterior expected loss）または事後リスク（posterior risk）と呼ばれていますが，ベイズリスクを最小にする推定量は，$X$の実現値$x$（すなわち観測データ）ごとに，事後リスクを最小にすれば良いことが分かります．事後リスクを最小にする，すなわちベイズリスクを最小にする推定量のことを，**ベイズ推定量**（Bayes estimator）と呼びます．

具体的に，二乗損失におけるベイズ推定量を考えてみます．事後リスクは，
```math
\begin{split}
& E\left[ \{ \delta(X) - \theta \}^2 \mid X \right] \\
=& \{\delta(X)\}^2 - 2E[\theta \mid X] \cdot \delta(X) + E[\theta^2 \mid X] \\
=& \left\{ \delta(X) - E[\theta \mid X] \right\}^2 + E[\theta^2 \mid X] - \{ E[\theta \mid X] \}^2
\end{split}
```
となり，これを最小にするのは，$\delta(X) = E[\theta \mid X]$，すなわち事後平均であることが分かります．また，絶対誤差損失$L(\theta,d) = |d - \theta|$においては，ベイズ推定量は事後中央値になります．





# 確率分布の推測
次に，パラメータの点推定ではなく，確率分布の推測問題を考えます．前節と同様，観測データ$X_1,\dots,X_n$は，パラメータ$\theta$を所与として独立同一に密度$f(\cdot \mid \theta)$を持つ分布に従うとし，パラメータ$\theta$の事前密度を$\pi(\theta)$とします．また，$\theta$を所与とした$X = (X_1,\dots,X_n)$の密度関数（尤度）を，$f_n(x \mid \theta) = \prod_{i=1}^n f(x_i \mid \theta)$と書きます．

ここで，パラメータ$\theta$を所与として，$X_1,\dots,X_n$と独立に$f(\cdot \mid \theta)$に従う確率変数$Y$（仮想的な将来の値）の従う確率密度（すなわち$f(y \mid \theta)$）を，データ$X$を観測した後にどのような分布で推測すれば良いかを考えます．これを予測分布（予測密度）と言い，ベイズ統計学の枠組みでは，以下の**事後予測密度**（posterior predictive density）
```math
f^\ast(y \mid x) = \int f(y \mid \theta) \pi(\theta \mid x) d\theta
```
を用います．ただし，$\pi(\theta \mid x)$は事後密度です．事後予測密度は，統計的決定理論の枠組みではどのように正当化されるのでしょうか．

ある予測密度$p(y \mid x)$で$f(y \mid \theta)$を予測したときの損失として，以下のKullback–Leibler（KL）ダイバージェンスを考えます．
```math
L(\theta,p) := \int \log \left\{ \frac{f(y \mid \theta)}{p(y \mid x)} \right\} f(y \mid \theta) dy
```
これを標本$X$の分布$f_n(x \mid \theta)$で期待値をとると，予測密度のリスクになります．
```math
R(\theta,p) := \int \left[ \int \log \left\{ \frac{f(y \mid \theta)}{p(y \mid x)} \right\} f(y \mid \theta) dy \right] f_n(x \mid \theta) dx
```
このリスクは$\theta$に依存しますが，ベイズ法の文脈では，さらにこれを$\theta$の事前分布で期待値をとったベイズリスクを考えます．
```math
r(\pi,p) := \iint \left[ \int \log \left\{ \frac{f(y \mid \theta)}{p(y \mid x)} \right\} f(y \mid \theta) dy \right] f_n(x \mid \theta) \pi(\theta) dx d\theta
```
2つの予測密度$p(y \mid x)$と$q(y \mid x)$のベイズリスクの差は，
```math
\begin{equation}
\tag{3}
\begin{split}
&r(\pi,p) - r(\pi,q) \\
=& \iint \left[ \int \log \left\{ \frac{q(y \mid x)}{p(y \mid x)} \right\} f(y \mid \theta) dy \right] f_n(x \mid \theta) \pi(\theta) dx d\theta \\
=&: D(p,q)
\end{split}
\end{equation}
```
であり，この値が正であれば，$q(y \mid x)$は$p(y \mid x)$より良い予測密度であると言えます．点推定量のベイズリスクを変形したときと同様に，$f_n(x\mid \theta) \pi(\theta) = \pi(\theta \mid x) f_\pi(x)$が成り立つことに気をつけ，(3)式の積分の順序交換を行うと，
```math
\begin{equation}
\tag{4}
\begin{split}
D(p,q) &= \iint \log \left\{ \frac{q(y \mid x)}{p(y \mid x)} \right\} \int f(y \mid \theta) \pi(\theta \mid x) d\theta f_\pi(x) dydx \\
&= \iint \log \left\{ \frac{q(y \mid x)}{p(y \mid x)} \right\} f^\ast(y \mid x) f_\pi(x) dydx
\end{split}
\end{equation}
```
といった式変形ができます．ただし$f^\ast(y \mid x)$は事後予測密度であることに気をつけてください．ここで，(4)式の$q(y \mid x)$を$f^\ast(y \mid x)$におきかえた
```math
D(p,f^\ast) = \iint \log \left\{ \frac{f^\ast(y \mid x)}{p(y \mid x)} \right\} f^\ast(y \mid x) f_\pi(x) dydx
```
を考えます．このとき，
```math
\int \log \left\{ \frac{f^\ast(y \mid x)}{p(y \mid x)} \right\} f^\ast(y \mid x) dy
```
はKLダイバージェンスであることから，$f^\ast(y \mid x)$以外の任意の予測密度$p(y \mid x)$に対して，$D(p,f^\ast) > 0$が言えます．すなわち，事後予測密度$f^\ast(y \mid x)$は，KLダイバージェンスを損失としたベイズリスクの意味で，最良な予測密度であることが示されました．

事後予測密度は，モデル診断（model checking）の場面でも有用で，実用上よく用いられます．また，KLダイバージェンスを損失とした予測密度のベイズリスクは，情報量規準の構成においても用いられています．



# おわりに
株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては[株式会社Nospare](https://bit.ly/38WtGZH) までお問い合わせください．
<a href="https://bit.ly/38WtGZH">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1023506/d5c625c0-d6f4-680b-d702-3c40f11d6c39.jpeg" alt="株式会社Nospare" width=400>
</a>