はじめに
この記事は、 RICORA Advent Calendar 2022 18日目です。
こんにちは、さぶろくです。
データに正規モデルを当てはめたときのベイズ推測についての記事です。
正規モデル
確率変数$ Y$が平均$ \theta$, 分散$ \sigma ^2$の正規分布に従うとは、$ Y$の密度関数が以下のように与えられることを言います。
p(y|\theta, \sigma^2) = \frac{1}{\sqrt {2\pi \sigma ^2}}\exp\{-\frac{1}{2}(\frac{y-\theta}{\sigma})^2\}, ~-\infty < y < \infty
このとき、$ Y \sim N(\theta, \sigma^2)$のように表記することにします。
また、正規分布の性質のうち重要なものを以下に示しておきます。
- $ \theta$に関して対称で、平均、中央値、最頻値はすべて$ \theta$に等しい。
- 分布の95%は平均から標準偏差の2倍（正確には、1.96倍）の距離の間にある。
- $ X \sim N(\mu, \tau^2), ~Y\sim N(\theta, \sigma^2), XとY$が独立なとき、
aX+bY \sim N(a\mu+b\theta, a^2\tau^2+b^2\sigma^2)
分散が与えられた下での平均の推測
モデルとして、$ {Y_1 ,\cdots, Y_n}\sim i.i.d. N(\theta, \sigma^2)$を考えよう。
このとき、同時密度関数は以下のようになります。
\begin{align}
p(y_1, \cdots, y_n|\theta, \sigma^2) &= \prod_{i=1}^n p(y_i|\theta, \sigma^2)\\
&= \prod_{i=1}^n \frac{1}{\sqrt {2\pi \sigma ^2}}\exp\{-\frac{1}{2}(\frac{y_i-\theta}{\sigma})^2\}\\
&= (2\pi \sigma^2)^{-\frac{n}{2}}\exp\{-\frac{1}{2}\sum(\frac{y_i-\theta}{\sigma})^2\}
\end{align}
正規モデルはパラメータを2つ持つ。
2つのパラメータを推測する場合、1つのパラメータを推測する操作を2回行うということを考えます。
そこで、$ \sigma^2$が既知の場合の$ \theta$の推測から考え、$ \theta$の共役事前分布を使って事後分布を導出する。
ただし、共役事前分布とは事後分布と同じクラスに属するような事前分布のことを言います。
任意の事前分布$ p(\theta|\sigma^2)$について、事後分布は以下の条件を満たします。
\begin{align}
p(\theta|y_1, \cdots, y_n, \sigma^2) &\propto p(\theta|\sigma^2)\times \exp\{-\frac{1}{2\sigma^2}\sum(y_i-\theta)^2\}\\
&\propto p(\theta | \sigma^2) \times \exp\{c_1(\theta - c_2)^2\}
\end{align}
上の計算より、$ p(\theta|\sigma^2)$が共役になるためには、同じクラスに属する必要があるので、少なくとも$e^{c_1(\theta - c_2)^2}$のような項を含んでいなければなりません。
その際に、採用される分布が正規分布です。
それでは、$p(\theta|\sigma^2)$が正規分布かつ、$ y_1, \cdots, y_n|\theta, \sigma^2 \sim i.i.d.N(\theta, \sigma^2)$である場合の事後分布$ p(\theta | y_1, \cdots, y_n, \sigma^2)$もまた正規分布になることを確かめてみましょう。
$ \theta \sim N(\mu_0, \tau_0^2)$であるとき、
\begin{align}
p(\theta|y_1, \cdots , y_n, \sigma^2) &= \frac{p(\theta|\sigma^2)p(y_1, \cdots, y_n|\theta, \sigma^2)}{p(y_1, \cdots, y_n)}(\because Bayes' rule)\\
&\propto p(\theta|\sigma^2)p(y_1, \cdots, y_n|\theta, \sigma^2)\\
&\propto \exp\{-\frac{1}{2\tau_0^2}(\theta-\mu_0)^2\}\exp\{-\frac{1}{2\sigma^2}\sum(y_i-\theta)^2\}\\
&=\exp[-\frac{1}{2}\{\frac{1}{\tau_0^2}(\theta-\mu_0)^2+\frac{1}{\sigma^2}\sum(y_i-\theta)^2\}]\\
\end{align}
となる。
ここで、
\frac{1}{\tau_0^2}(\theta-\mu_0)^2+\frac{1}{\sigma^2}\sum(y_i-\theta)^2=a\theta^2-2b\theta+c,\\
a = \frac{1}{\tau_0^2}+\frac{n}{\sigma^2}, b = \frac{\mu_0}{\tau_0^2}+\frac{\sum y_i}{\sigma ^2}, c=\frac{\mu_0^2}{\tau_0^2}+\frac{\sum y_i^2}{\sigma^2}
とおく。
このとき、
\begin{align}
p(\theta|\sigma^2, y_1, \cdots, y_n) &\propto \exp\{-\frac{1}{2}(a\theta^2-2b\theta)\}\\
&= \exp \{ -\frac{1}{2}a(\theta-\frac{b}{a})^2+\frac{1}{2}\frac{b^2}{a}\}\\
&\propto \exp \{ -\frac{1}{2}a(\theta-\frac{b}{a})^2\}\\
&= \exp\{-\frac{1}{2}(\frac{\theta-\frac{b}{a}}{\frac{1}{\sqrt a}})^2\}
\end{align}
を得る。
この形は正規分布の密度関数なので、$ p(\theta|\sigma^2, y_1, \cdots, y_n)$が正規分布の密度関数になることが確かめられました。
この分布の平均と分散をそれぞれ$ \mu_n$と$ \tau_n^2$として、
\tau_n^2 = \frac{1}{a}=\frac{1}{\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}}, ~~\mu_n = \frac{b}{a} =
\frac{\frac{1}{\tau_0^2}\mu_0+\frac{n}{\sigma^2}\bar{y}}{\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}}
のようにおけます。
ここで、事後分散の逆数$ \frac{1}{\tau_n^2}$は以下の式で表せます。
\frac{1}{\tau_n^2} = \frac{1}{\tau_0^2} + \frac{n}{\sigma^2}
分散の逆数のことを精度と言います。
この式より、事後分布の精度は事前分布の精度と標本平均の標本分布の精度からなる式で表せるということがわかります。
事後平均$ \mu_n$は、
~~\mu_n=
\frac{\frac{1}{\tau_0^2}\mu_0+\frac{n}{\sigma^2}\bar{y}}{\frac{1}{\tau_0^2}+\frac{n}{\sigma^2}}
となり、これは事前分布の精度と標本平均の標本分布の精度で重みづけられた、事前平均と標本平均の和となっています。
事前平均を$ Y_1, \cdots, Y_n$と同じ母集団から得られた$ \kappa _0$個の観測値の平均とみなすと、事前の観測値の平均の分散は$ \frac{\sigma^2}{\kappa _0}$のように表せます。
そこで、事前分散を$\tau_0^2 = \frac{\sigma^2}{\kappa _0}$としてみると、事後平均の式は
\mu_n = \frac{\kappa_0}{\kappa_0 +n}\mu_0+\frac{n}{\kappa_0 +n}\bar{y}
となり、事前のサンプルサイズと標本のサンプルサイズによって重みづけされた、事前平均と標本平均の和で表せます。
平均と分散の同時推定
$(\theta, \sigma^2)$の同時事後密度関数$ p(\theta, \sigma^2|y_1, \cdots, y_n)$を求めよう。
まず、共役な事前分布のクラスを求めるために、以下の確率の公理を用いる。
p(\theta, \sigma^2)=p(\theta|\sigma^2)p(\sigma^2)
$ p(\theta|\sigma^2)$に関しては$ N(\mu_0, \tau_0^2)$の密度関数となる。
ここで事前平均は事前の$ \kappa_0$個の観測値の平均とみなすことにし、$ \tau_0^2 = \frac{\sigma^2}{\kappa _0}$となる。
このとき、任意の事前分布$p(\theta, \sigma^2)$は以下のように計算できる。
p(\theta, \sigma^2) = p(\theta|\sigma^2)p(\sigma^2)=dnorm(\theta, \mu_0, \tau_0=\frac{\sigma}{\sqrt\kappa_0})\times p(\sigma^2),\\
dnorm(\theta, \mu_0, \tau_0=\frac{\sigma}{\sqrt\kappa_0})=\exp\{-\frac{1}{2\tau_0^2}(\theta-\mu_0)^2\}
今後、正規分布の密度関数を$ dnorm(\cdot, \cdot, \cdot)$であらわすことにする。
$ \sigma^2$について、$\sigma^2 > 0$より台が$(0, \infty)$となる事前分布族が必要になります。
その一例はガンマ分布族です。
しかし、この分布族は正規分布の分散について共役になりません。
ただ精度$(=\frac{1}{\sigma^2})$については共役であることが結果としてわかっています。
このとき、$ \sigma^2$は逆ガンマ分布に従うといいます。
ここまでの議論により、事前分布と標本モデルは以下のようになります。
\begin{align}
\frac{1}{\sigma^2}&\sim Ga(\frac{\nu_0^2}{2}, \frac{\nu_0^2}{2}\sigma_0^2),\\
\theta | \sigma^2 &\sim N(\mu_0, \frac{\sigma^2}{\kappa _0}),\\
Y_1, \cdots, Y_n &\sim i.i.d.N(\theta, \sigma^2)
\end{align}
事後分布においても、事前分布と同様な分解ができて、
p(\theta, \sigma^2|y_1, \cdots, y_n) = p(\theta|\sigma^2, y_1, \cdots, y_n)p(\sigma^2|, y_1, \cdots, y_n)
ここで、$p(\theta|\sigma^2, y_1, \cdots, y_n)$に関しては、
p(\theta|\sigma^2, y_1, \cdots, y_n) = dnorm(\theta, \mu_n, \frac{\sigma}{\sqrt\kappa_n}),\\
\kappa_n = \kappa_0 + n, \mu_n = \frac{\kappa_0\mu_0+ n\bar{y}}{\kappa _n}
となる。
$\sigma^2$の事後分布は
\begin{align}
p(\sigma^2|y_1, \cdots, y_n)&\propto p(\sigma^2)p(y_1, \cdots, y_n|\sigma^2)
\\&= p(\sigma^2)\int p(y_1, \cdots, y_n|\theta, \sigma^2)p(\theta|\sigma^2)d \theta\\
\end{align}
ここで、$p(y_1, \cdots, y_n|\theta, \sigma^2)p(\theta|\sigma^2)$を計算してみます。
\begin{align}
p(y_1, \cdots, y_n|\theta, \sigma^2)p(\theta|\sigma^2) &\propto (2\pi \sigma^2)^{-\frac{n}{2}}\exp\{-\frac{1}{2\sigma^2}\sum(y_i-\theta)^2\}\times(2\pi\tau_0^2)^\frac{1}{2}\exp\{-\frac{1}{2\tau_0^2}(\theta-\mu_0)^2\}\\
&\propto(2\pi \sigma^2)^{-\frac{n+1}{2}}\exp\{-\frac{\kappa_0}{2\sigma^2}(\theta-\mu_0)^2\}\exp\{-\frac{1}{2\sigma^2}\sum(y_i-\theta)^2\}\\
&=(2\pi \sigma^2)^{-\frac{n+1}{2}}\exp[-\frac{1}{2}\{\frac{\kappa_0}{\sigma^2}(\theta-\mu_0)^2+\frac{1}{\sigma^2}\sum(y_i-\theta)^2\}]\\
&=(2\pi \sigma^2)^{-\frac{n+1}{2}}\exp[-\frac{1}{2\sigma^2}\{\kappa_n\theta^2 - 2(\kappa_0\mu_0+\sum y_i)\theta+(\kappa_0\mu_0^2+\sum y_i^2)\}]\\
&=(2\pi \sigma^2)^{-\frac{n+1}{2}}\exp[-\frac{1}{2\sigma^2}\{\kappa_n(\theta - \frac{\kappa_0 \mu_0 + n\bar{y}}{\kappa_n})^2-\frac{(\kappa_0\mu_0+n\bar{y})^2}{\kappa_n}+\frac{\kappa_n(\kappa_0\mu_0^2+\sum y_i^2)}{\kappa_n}\}]\\
&=(2\pi \sigma^2)^{-\frac{n+1}{2}}\exp[-\frac{1}{2\sigma^2}\{\kappa_n(\theta - \frac{\kappa_0 \mu_0 + n\bar{y}}{\kappa_n})^2+\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]\\
&=(2\pi\sigma^2)^{-\frac{1}{2}}\exp\{-\frac{\kappa_n}{2\sigma^2}(\theta - \frac{\kappa_0 \mu_0 + n\bar{y}}{\kappa_n})^2\}\\
&\times
(2\pi\sigma^2)^{-\frac{n}{2}}\exp[-\frac{1}{2\sigma^2}\{\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]\\
&\propto dnorm(\theta |\sigma^2,y_1, \cdots, y_n, \mu_n, \frac{\sigma}{\sqrt{\kappa_n}})\\ &\times
(2\pi\sigma^2)^{-\frac{n}{2}}\exp[-\frac{1}{2\sigma^2}\{\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]
\end{align}
したがって、$ \frac{1}{\sigma^2} \sim Ga(\frac{\nu_0^2}{2}, \frac{\nu_0^2}{2}\sigma_0^2)$の変数変換より$p(\sigma^2)$の分布も含めて、
\begin{align}
p(\sigma^2|y_1, \cdots, y_n)&\propto p(\sigma^2)\int p(y_1, \cdots, y_n|\theta, \sigma^2)p(\theta|\sigma^2)d \theta\\
&\propto p(\sigma^2)\times(\sigma^2)^{-\frac{n}{2}}\exp[-\frac{1}{2\sigma^2}\{\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]\\
&\times\int dnorm(\theta |\sigma^2,y_1, \cdots, y_n, \mu_n, \frac{\sigma}{\sqrt{\kappa_n}})d \theta\\
&\propto \frac{(\frac{\nu_0 \sigma_0^2}{2})^\frac{\nu_0}{2}}{\Gamma (\frac{\nu_0}{2})}(\sigma ^2)^{-\frac{\nu_0}{2}-1}\exp\{-\frac{\nu_0 \sigma_0 ^2}
{2\sigma^2}\}\\
&\times(\sigma^2)^{-\frac{n}{2}}\exp[-\frac{1}{2\sigma^2}\{\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]\\
&\propto (\sigma ^2)^{-\frac{\nu_n}{2}-1}\exp[-\frac{\nu_n}{2\nu_n \sigma^2}\{
\nu_0\sigma_0^2+\sum(y_i-\bar{y})^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}\}]\\
&=(\sigma ^2)^{-\frac{\nu_n}{2}-1}\exp\{-\frac{\nu_n \sigma_n^2}{2\sigma^2}\}
\end{align}
となる。
ただし、
\begin{align}
\nu_n &= \nu_0 + n,\\
s^2 &= \frac{1}{n-1}\sum(y_i-\bar{y})^2,\\
\sigma_n^2 &= \frac{1}{\nu_n}[
\nu_0\sigma_0^2+(n-1)s^2+\frac{n\kappa_0(\bar{y}-\mu_0)^2}{\kappa_n}]
\end{align}
である。
これより、${\frac{1}{\sigma ^2}|y_1, \cdots, y_n}\sim Ga(\frac{\nu_n}{2}, \frac{\nu_n \sigma_n^2}{2})$となることがわかりました。
まとめると、事後分布は
\begin{align}
\sigma^2|y_1, \cdots, y_n &\sim inverse-Ga(\frac{\nu_n}{2}, \frac{\nu_n \sigma_n^2}{2})\\
\theta|\sigma^2, y_1, \cdots, y_n &\sim N(\mu_n, \frac{\sigma^2}{\kappa_n})\\
\end{align}
となります。
標本平均と事後平均推定量の比較
未知パラメータ$\theta$の点推定量とは、データをパラメータ空間$\Theta$の一点に集約できる関数です。
正規分布による標本モデルと共役事前分布において、$ \theta$の事後平均推定量は次のようになります。
\hat\theta_b(y_1, \cdots, y_n) = \mathbb E[\theta |y_1, \cdots, y_n] = \mu_n = \frac{n}{\kappa_0 +n}\bar{y}+\frac{\kappa_0}{\kappa_0 +n}\mu_0=w\bar{y}+(1-w)\mu_0
調査や実験が何度も繰り返された場合の推定量の振る舞いを標本の性質と言います。
母集団平均の真の値を$\theta_0$として、$\hat \theta_b$と標本平均$\hat\theta_e(y_1, \cdots, y_n) = \bar{y}$の標本の性質を比べてみよう。
$
\mathbb E[\hat\theta_e|\theta = \theta_0] = \theta_0$より、$\hat\theta_e$は不偏推定量である。
$\mathbb E[\hat\theta_b|\theta = \theta_0] = w\theta_0 + (1-w)\mu_0, \mu_0 \neq \theta_0$ならば、$\hat \theta_b$は不偏推定量ではない。
推定量の標本分布の中心が真の値にどれだけ近いかを表す指標をバイアスと言います。
不偏推定量はバイアスが0になるので、$\hat \theta_e$のほうが優れていると結論付けたくなるかもしれません。
しかし、バイアスだけでは推定値と真の値がどれだけ離れているかわかりません。
そこで、推定量$\hat \theta$が真の値$\theta_0$の近くにあるかどうかを評価するために平均二乗誤差(MSE)を使います。
$m=\mathbb E[\hat \theta|\theta_0]$とおくと、MSEは
\begin{align}
MSE[\hat\theta|\theta] &= \mathbb E[(\hat \theta - \theta_0)^2|\theta_0]\\
&= \mathbb E[(\hat \theta -m + m - \theta_0)^2|\theta_0]\\
&= \mathbb E[(\hat \theta - m)^2|\theta_0] + 2\mathbb E[(\hat \theta - m)(m - \theta_0)|\theta_0] + \mathbb E[(m- \theta_0)^2|\theta_0]\\
&= \mathbb Var[\hat \theta|\theta_0] + Bias^2[\hat \theta | \theta_0]
\end{align}
となる。
MSEの形より、$\hat \theta$と真の値$\theta_0$の距離は、$\theta_0$が$\hat \theta$の分布の中心にどれだけ近いか（バイアス）と、分布の散らばり（分散）に依存していることがわかります。
それでは、$\theta_b$と$\theta_e$を再度比較しましょう。$\theta_e$のバイアスは0でしたが、
\begin{align}
\mathbb Var[\hat \theta_e | \theta = \theta_0, \sigma^2] &= \frac{\sigma^2}{n},\\
\mathbb Var[\hat \theta_b | \theta = \theta_0, \sigma^2] &= w^2 \times \frac{\sigma^2}{n} < \frac{\sigma^2}{n}
\end{align}
となり、$\hat\theta_b$のほうが分散は小さい。
MSEは、
\begin{align}
MSE[\hat\theta_e|\theta_0] = Var[\hat \theta_e | \theta = \theta_0, \sigma^2] &= \frac{\sigma^2}{n}\\
MSE[\hat\theta_b|\theta_0] = \mathbb E[(\hat \theta_b-\theta_0)^2|\theta_0] &= \mathbb E[\{w(\bar{y}-\theta_0) + (1-w)(\mu_0-\theta_0)\}^2|\theta_0]\\
&=w^2\times \frac{\sigma^2}{n}+(1-w)^2(\mu_0-\theta_0)^2
\end{align}
差分をとってみると、
\begin{align}
MSE[\hat\theta_e|\theta_0] - MSE[\hat\theta_b|\theta_0]&=\frac{\sigma^2}{n}-w^2\times \frac{\sigma^2}{n}-(1-w)^2(\mu_0-\theta_0)^2\\
&=(1-w)(1+w)\times \frac{\sigma^2}{n} - (1-w)^2(\mu_0- \theta_o)^2\\
\end{align}
ここで、$MSE[\hat\theta_e|\theta_0]> MSE[\hat\theta_b|\theta_0]$となるのは、
\begin{align}
(\mu_0- \theta_o)^2 &< \frac{\sigma^2}{n}\frac{1+w}{1-w}\\
&= \sigma^2(\frac{1}{n}+ \frac{2}{\kappa_0})
\end{align}
が成り立つときです。
この結果より、ベイズ推定をする場合は、母集団について事前に知っていることがあれば、上の不等式を満たすような$ \mu_0$と$ \kappa _0$の値を見つけるべきだということがわかりました。
最後に
これで当記事は終わりです。
ここでの内容は、標準ベイズ統計学の5章に書いてあるので、興味を持った方は読んでみると良いと思います。
また、標準ベイズの演習問題5.3が分散の事後分布の導出になっているのですが、解答は見つからず、計算も骨が折れるので躓いた方はこの記事を参考にするとよいと思います。
参考文献