共役事前分布によるベイズ推定

statistics

正規分布

平均

ベイズ推定

尤度関数を正規分布 $\mathcal{N}(x; \mu, \sigma^{2})$ と仮定して、$\mu$ の共役事前分布を正規分布 $\mathcal{N}(\mu; \mu_{0}, \sigma_{0}^{2})$ とする。あたらしくサンプル $x$ が得られたときの、$\mu_{0}$ と $\sigma_{0}^{2}$ の更新式を求めたい。 ただし $\sigma^{2}$ は既知とする。

事後分布は

\begin{aligned}
\mathcal{N}(\mu; \mu_{1}; \sigma_{1}^{2})
  &\propto \mathcal{N}(x; \mu, \sigma^{2}) \cdot
           \mathcal{N}(\mu; \mu_{0}, \sigma_{0}^{2}) \\

  &= \left( \frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{-\frac{(x - \mu)^2}{2 \sigma^{2}}} \right)
     \left( \frac{1}{\sqrt{2 \pi \sigma_{0}^{2}}} e^{-\frac{(\mu - \mu_{0})^2}{2 \sigma_{0}^{2}}} \right) \\

  &\propto \left( e^{-\frac{(x - \mu)^2}{2 \sigma^{2}}} \right)
           \left( e^{-\frac{(\mu - \mu_{0})^2}{2 \sigma_{0}^{2}}} \right) \\

  &= e^{-\frac{1}{2} \left(\frac{(x - \mu)^2}{\sigma^{2}} + \frac{(\mu - \mu_{0})^2}{\sigma_{0}^{2}} \right)} \\

  &= e^{-\frac{1}{2} \left(\frac{\sigma_{0}^{2} (x - \mu)^2 + \sigma^{2} (\mu - \mu_{0})^2}{\sigma_{0}^{2} \sigma^{2}} \right)} \\

  &= e^{-\frac{1}{2} \left(\frac{\sigma_{0}^{2} x^{2} - 2 \sigma_{0}^{2} x \mu + \sigma_{0}^{2} \mu^{2} + \sigma^{2} \mu^{2} - 2 \sigma^{2} \mu \mu_{0} + \sigma^{2} \mu_{0}^2}{\sigma_{0}^{2} \sigma^{2}} \right)} \\

  &= e^{-\frac{1}{2} \left(\frac{ (\sigma_{0}^{2} + \sigma^{2}) \mu^{2} - 2 (\sigma_{0}^{2} x + \sigma^{2} \mu_{0}) \mu + (\sigma_{0}^{2} x^{2} + \sigma^{2} \mu_{0}^2)}{\sigma_{0}^{2} \sigma^{2}} \right)} \\

  &= e^{-\frac{1}{2} \left(\frac{ \mu^{2} - 2 \frac{\sigma_{0}^{2} x + \sigma^{2} \mu_{0}}{\sigma_{0}^{2} + \sigma^{2}} \mu + \frac{\sigma_{0}^{2} x^{2} + \sigma^{2} \mu_{0}^2}{\sigma_{0}^{2} + \sigma^{2}}}{\frac{\sigma_{0}^{2} \sigma^{2}}{\sigma_{0}^{2} + \sigma^{2}}} \right)} \\

  &= e^{-\frac{1}{2} \left(\frac{ \mu^{2} - 2 \frac{\frac{x}{\sigma^{2}} + \frac{\mu_{0}}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}} \mu + \frac{\frac{x^{2}}{\sigma^{2}} + \frac{\mu_{0}^2}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}}}{\frac{1}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}}} \right)} \\

  &= e^{-\frac{1}{2} \left( \frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}} \right) \left( \mu^{2} - 2 \frac{\frac{x}{\sigma^{2}} + \frac{\mu_{0}}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}} \mu + \frac{\frac{x^{2}}{\sigma^{2}} + \frac{\mu_{0}^2}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}} \right)} \\

  &\propto e^{-\frac{1}{2} \left( \frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}} \right) \left( \mu^{2} - 2 \frac{\frac{x}{\sigma^{2}} + \frac{\mu_{0}}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}} \mu \right)} \\

  &\propto e^{-\frac{1}{2} \left( \frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}} \right) \left( \mu - \frac{\frac{x}{\sigma^{2}} + \frac{\mu_{0}}{\sigma_{0}^{2}}}{\frac{1}{\sigma^{2}} + \frac{1}{\sigma_{0}^{2}}} \right)^{2}}
\end{aligned}

よって

\begin{aligned}
\mu_{1} &= \frac{\frac{\mu_{0}}{\sigma_{0}^{2}} + \frac{x}{\sigma^{2}}}{\frac{1}{\sigma_{0}^{2}} + \frac{1}{\sigma^{2}}} \\
\frac{1}{\sigma_{1}^{2}} &= \frac{1}{\sigma_{0}^{2}} + \frac{1}{\sigma^{2}}
\end{aligned}

となる。

複数サンプルの場合

TODO

分散

平均 $\mu$ は既知とする。

ガンマ分布

\mathcal{G}(x; k, \theta) = x^{k-1} \frac{e^{-\frac{x}{\theta}}}{\Gamma(k) \theta^{k}}

$x > 0$ がガンマ分布に従う。$k > 0$ と $\theta > 0$ はパラメータ。

逆ガンマ分布

$x$ がガンマ分布に従うとき、$\frac{1}{x}$ が従う分布が逆ガンマ分布。

\mathcal{G}^{-1}(x; k , \theta) = x^{-k-1} \frac{e^{-\frac{\theta}{x}}}{\Gamma(k) \theta^{-k}}

$x > 0$ が逆ガンマ分布に従う。$k > 0$ と $\theta > 0$ はパラメータで、ガンマ分布のものと意味は異なる ($k$ はそのまま、$\theta$ は逆数をとれば意味的には一致する)。

ベイズ推定において、正規分布の分散の事前分布として利用される。

ベイズ推定

尤度関数を正規分布 $\mathcal{N}(x; \mu, \sigma^{2})$ と仮定して、その分散 $\sigma^{2}$ の共役事前分布を逆ガンマ分布 $\mathcal{G}^{-1}(\sigma^{2}; \frac{k}{2}, \frac{\theta}{2})$ とする。あたらしくサンプル $x$ が得られたときの、$k$ と $\theta$ の更新式を求めたい。ただし $\mu$ は既知とする。

\begin{aligned}
\mathcal{G}^{-1}\left(\sigma^{2}; \frac{k'}{2}, \frac{\theta'}{2}\right)

  &\propto \mathcal{N}(x; \mu, \sigma^{2}) \cdot
           \mathcal{G}^{-1}\left( \sigma^{2}; \frac{k}{2}, \frac{\theta}{2} \right) \\

  &\propto \left( \frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{-\frac{(x - \mu)^{2}}{2 \sigma^{2}}} \right) \left( (\sigma^{2})^{-\frac{k}{2}-1} \frac{e^{-\frac{\theta}{2 \sigma^{2}}}}{\Gamma(k) \theta^{-k}} \right) \\

  &\propto \left( (\sigma^{2})^{-\frac{1}{2}} e^{-\frac{(x - \mu)^{2}}{2 \sigma^{2}}} \right) \left( (\sigma^{2})^{-\frac{k}{2}-1} e^{-\frac{\theta}{2 \sigma^{2}}} \right) \\

  &= (\sigma^{2})^{-\frac{k + 1}{2}-1} e^{-\frac{\theta + (x - \mu)^{2}}{2 \sigma^{2}}}
\end{aligned}

よって

\begin{aligned}
k' &= k + 1 \\
\theta' &= \theta + (x - \mu)^{2}
\end{aligned}

となる。

複数サンプルのとき

$N$ 個のサンプル $\{x_{i}\}_{i=1}^{N}$ が与えられるとする。尤度関数が $\prod_{i=1}^{N} \mathcal{N}(x_{i}; \mu, \sigma^{2})$ となる以外は $1$ 個のときと同じになる。

\begin{aligned}
\mathcal{G}^{-1}\left(\sigma^{2}; \frac{k'}{2}, \frac{\theta'}{2}\right)

  &\propto \left(\prod_{i=1}^{N} \mathcal{N}(x_{i}; \mu, \sigma^{2}) \right) \cdot \mathcal{G}^{-1}\left( \sigma^{2}; \frac{k}{2}, \frac{\theta}{2} \right) \\

  &\propto \left( \left(\frac{1}{\sqrt{2 \pi \sigma^{2}}}\right)^{N} e^{-\frac{\sum_{i=1}^{N} (x_{i} - \mu)^{2}}{2 \sigma^{2}}} \right) \left( (\sigma^{2})^{-\frac{k}{2}-1} \frac{e^{-\frac{\theta}{2 \sigma^{2}}}}{\Gamma(k) \theta^{-k}} \right) \\

  &\propto \left( (\sigma^{2})^{-\frac{N}{2}} e^{-\frac{\sum_{i=1}^{N} (x_{i} - \mu)^{2}}{2 \sigma^{2}}} \right) \left( (\sigma^{2})^{-\frac{k}{2}-1} e^{-\frac{\theta}{2 \sigma^{2}}} \right) \\

  &= (\sigma^{2})^{-\frac{k + N}{2}-1} e^{-\frac{\theta + \sum_{i=1}^{N} (x_{i} - \mu)^{2}}{2 \sigma^{2}}}
\end{aligned}

よって

\begin{aligned}
k' &= k + N \\
\theta' &= \theta + \sum_{i=1}^{N} (x_{i} - \mu)^{2}
\end{aligned}

となる。

平均と分散

正規逆ガンマ分布

正規分布の平均 $\mu$ と分散 $\sigma^{2}$ の共役事前分布として利用される。

p(\mu, \sigma^{2}) = p(\mu \mid \sigma^{2}) p(\sigma^{2})

より

\begin{aligned}
\mathcal{N}\textrm{-}\mathcal{G}^{-1}\left(\mu, \sigma^{2}; \mu_{0}, \gamma, \frac{k}{2}, \frac{\theta}{2}\right)

 &= \mathcal{N}\left(\mu; \mu_{0}, \frac{\sigma^{2}}{\gamma} \right)
     \mathcal{G}^{-1}\left(\sigma^{2}; \frac{k}{2}, \frac{\theta}{2} \right)
 \\

 &= \left(
      \frac{1}{\sqrt{2 \pi \frac{\sigma^{2}}{\gamma}}}
        e^{-\frac{\gamma (\mu - \mu_{0})^{2}}{2 \sigma^{2}}}
    \right)
    \left(
      (\sigma^{2})^{-\frac{k}{2}-1}
      \frac{e^{-\frac{\theta}{2 \sigma^{2}}}}{\Gamma(k) \theta^{-k}}
    \right)
 \\

 &\propto (\sigma^{2})^{-\frac{k+1}{2} - 1}
   e^{-\frac{\theta + \gamma (\mu - \mu_{0})^2}{2 \sigma^{2}}}

\end{aligned}

と定義される。

ベイズ推定

尤度関数を正規分布 $\mathcal{N}(x; \mu, \sigma^{2})$ と仮定して、$\mu$ と $\sigma^{2}$ の共役事前分布を正規逆ガンマ分布 $\mathcal{N}\textrm{-}\mathcal{G}^{-1}(\mu, \sigma^{2}; \mu_{0}, \gamma_{0}, \frac{k_{0}}{2}, \frac{\theta_{0}}{2})$ とする。あたらしくサンプル $x$ が得られたときの、$\mu_{0}$、$\gamma_{0}$、$k_{0}$、$\theta_{0}$ の更新式を求めたい。

事後分布は

\begin{aligned}
\mathcal{N}\textrm{-}\mathcal{G}^{-1}\left(\mu, \sigma^{2};
  \mu_{1}, \gamma_{1}, \frac{k_{1}}{2}, \frac{\theta_{1}}{2}\right)

 &\propto \mathcal{N}(x; \mu, \sigma^{2})
    \mathcal{N}\textrm{-}\mathcal{G}^{-1}\left(\mu, \sigma^{2};
      \mu_{0}, \gamma_{0}, \frac{k_{0}}{2}, \frac{\theta_{0}}{2} \right) 
\\

 &\propto
    \left(
      \frac{1}{\sqrt{2 \pi \sigma^{2}}}
        e^{-\frac{(x - \mu)^2}{2 \sigma^{2}}}
    \right)
    \left(
      (\sigma^{2})^{-\frac{k_{0}+1}{2} - 1}
        e^{-\frac{\theta_{0} + \gamma_{0} (\mu - \mu_{0})^2}{2 \sigma^{2}}}
    \right)
 \\

 &\propto
    (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{\theta_{0} + \gamma_{0} (\mu - \mu_{0})^2 + (x - \mu)^2}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k + 2}{2}-1}
      e^{-\frac{\theta_{0} + \gamma_{0} \mu^{2} - 2 \gamma_{0} \mu \mu_{0} + \gamma_{0} \mu_{0}^{2} + x^{2} - 2 x \mu + \mu^2}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + \gamma_{0} \mu^{2} - 2 \gamma_{0} \mu \mu_{0} + \gamma_{0} \mu_{0}^{2}
               + x^{2} - 2 x \mu + \mu^2}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \mu^{2} - 2 (\gamma_{0} \mu_{0} + x) \mu
               + (\gamma_{0} \mu_{0}^{2} + x^{2})}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu^{2} - 2 \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \mu \right)
               + (\gamma_{0} \mu_{0}^{2} + x^{2})}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               - \frac{(\gamma_{0} \mu_{0} + x)^{2}}{\gamma_{0} + 1}
               + (\gamma_{0} \mu_{0}^{2} + x^{2})}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               - \frac{(\gamma_{0} \mu_{0} + x)^{2}}{\gamma_{0} + 1}
               + \frac{(\gamma_{0} + 1) (\gamma_{0} \mu_{0}^{2} + x^{2})}{\gamma_{0} + 1}}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               + \frac{(\gamma_{0} + 1) (\gamma_{0} \mu_{0}^{2} + x^{2}) - (\gamma_{0} \mu_{0} + x)^{2}}{\gamma_{0} + 1}}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               + \frac{\gamma_{0}^{2} \mu_{0}^{2} + \gamma_{0} x^{2} + \gamma_{0} \mu_{0}^{2} + x^{2} - \gamma_{0}^{2} \mu_{0}^{2} - 2 \gamma_{0} \mu_{0} x - x^{2}}{\gamma_{0} + 1}}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               + \frac{\gamma_{0} x^{2} + \gamma_{0} \mu_{0}^{2} - 2 \gamma_{0} \mu_{0} x}{\gamma_{0} + 1}}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               + \frac{\gamma_{0}}{\gamma_{0} + 1} (x^{2} + \mu_{0}^{2} - 2 \mu_{0} x)}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
        \theta_{0} + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}
               + \frac{\gamma_{0}}{\gamma_{0} + 1} (x - \mu_{0})^{2}}{2\sigma^{2}}}
 \\

 &= (\sigma^{2})^{-\frac{k_{0} + 2}{2}-1}
      e^{-\frac{
         \left(\theta_{0} + \frac{\gamma_{0}}{\gamma_{0} + 1} (x - \mu_{0})^{2} \right) + (\gamma_{0} + 1) \left(\mu - \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \right)^{2}}{2\sigma^{2}}}
 \\

\end{aligned}

より

\begin{aligned}
\mu_{1} &= \frac{\gamma_{0} \mu_{0} + x}{\gamma_{0} + 1} \\
\gamma_{1} &= \gamma_{0} + 1 \\
k_{1} &= k_{0} + 1 \\
\theta_{1} &= \frac{\gamma_{0}}{\gamma_{0} + 1} (x - \mu_{0})^{2}
\end{aligned}

となる。

複数データの場合

TODO

多変量正規分布

平均ベクトル

$p$ 次元正規分布は

\mathcal{N}_{p}(\vec{x}; \vec{\mu}, \Sigma) =
  \frac{1}{(2 \pi)^\frac{p}{2} \sqrt{| \Sigma |}}
  e^{-\frac{1}{2}
    \left(\vec{x} - \vec{\mu} \right)^{\top} \Sigma^{-1}
    \left(\vec{x} - \vec{\mu} \right)}

と定義される。ここで $\vec{\mu}$ は $p$ 元平均ベクトル、$\Sigma$ は $p$ 行 $p$ 列の共分散行列である。共分散行列は対称行列なので、逆行列も対称行列になることに注意すると、あとの式変形の見通しがよくなる。

ベイズ推定

尤度関数を多変量正規分布 $\mathcal{N}_{p}(\vec{x}; \vec{\mu}, \Sigma)$ と仮定して、$\vec{\mu}$ の共役事前分布を多変量正規分布 $\mathcal{N}_{p}(\vec{\mu}; \vec{\mu}_{0}, \Sigma_{0})$ とする。あたらしくサンプル $\vec{x}$ が得られたときの、$\vec{\mu}_{0}$ と $\Sigma_{0}$ の更新式を求めたい。 ただし $\Sigma$ は既知とする。

事後確率は

\begin{aligned}

\mathcal{N}_{p}(\vec{\mu}; \vec{\mu}_{1}, \Sigma_{1})

 &\propto
    \mathcal{N}_{p}(\vec{x}; \vec{\mu}, \Sigma) \cdot
    \mathcal{N}_{p}(\vec{\mu}; \vec{\mu}_{0}, \Sigma_{0}) \\

 &\propto
    \left(
      e^{-\frac{1}{2}
        \left(\vec{x} - \vec{\mu} \right)^{\top} \Sigma^{-1}
        \left(\vec{x} - \vec{\mu} \right)}
    \right)
    \left(
      e^{-\frac{1}{2}
        \left(\vec{\mu} - \vec{\mu}_{0} \right)^{\top} \Sigma_{0}^{-1}
        \left(\vec{\mu} - \vec{\mu}_{0} \right)}
    \right) \\

 &\propto
    e^{-\frac{1}{2} \left\{
      \left(\vec{x} - \vec{\mu} \right)^{\top} \Sigma^{-1}
      \left(\vec{x} - \vec{\mu} \right) +
      \left(\vec{\mu} - \vec{\mu}_{0} \right)^{\top} \Sigma_{0}^{-1}
      \left(\vec{\mu} - \vec{\mu}_{0} \right)
    \right\}} \\

 &= e^{-\frac{1}{2} \left(
          \vec{x}^{\top} \Sigma^{-1} \vec{x}
        - 2 \vec{x}^{\top} \Sigma^{-1} \vec{\mu}
        + \vec{\mu}^{\top} \Sigma^{-1} \vec{\mu}
        + \vec{\mu}^{\top} \Sigma_{0}^{-1} \vec{\mu}
        - 2 \vec{\mu}_{0}^{\top} \Sigma_{0}^{-1} \vec{\mu}
        + \vec{\mu}_{0}^{\top} \Sigma_{0}^{-1} \vec{\mu}_{0}
    \right)} \\

 &= e^{-\frac{1}{2} \left\{
      \vec{\mu}^{\top} \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right) \vec{\mu}
        - 2 \left(\vec{x}^{\top} \Sigma^{-1} +
             \vec{\mu}_{0}^{\top} \Sigma_{0}^{-1} \right) \vec{\mu}
        + \left(
           \vec{x}^{\top} \Sigma^{-1} \vec{x} +
           \vec{\mu}_{0}^{\top} \Sigma_{0}^{-1} \vec{\mu}_{0}
        \right)
    \right\}} \\

 &\propto
    e^{-\frac{1}{2} \left\{
      \vec{\mu}^{\top} \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right) \vec{\mu}
        - 2 \left(\vec{x}^{\top} \Sigma^{-1} +
             \vec{\mu}_{0}^{\top} \Sigma_{0}^{-1} \right) \vec{\mu}
    \right\}} \\

 &=
    e^{-\frac{1}{2} \left\{
      \vec{\mu}^{\top} \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right) \vec{\mu}
        - 2 \vec{\mu}^{\top}
          \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)
          \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)^{-1}
          \left(\Sigma^{-1} \vec{x} +
             \Sigma_{0}^{-1} \vec{\mu}_{0} \right)
    \right\}} \\

 &\propto
    e^{-\frac{1}{2}
       \left\{ \vec{\mu} -
          \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)^{-1}
          \left(\Sigma^{-1} \vec{x}^{\top} +
             \Sigma_{0}^{-1} \vec{\mu}_{0} \right)
       \right\}^{\top}
       \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)
       \left\{ \vec{\mu} -
          \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)^{-1}
          \left(\Sigma^{-1} \vec{x} +
             \Sigma_{0}^{-1} \vec{\mu}_{0} \right)
       \right\}}
\end{aligned}

より

\begin{aligned}
\mu_{1} &= \left(\Sigma_{0}^{-1} + \Sigma^{-1} \right)^{-1}
           \left( \Sigma_{0}^{-1} \vec{\mu}_{0} +
              \Sigma^{-1} \vec{x} \right) \\
\Sigma_{1}^{-1} &= \Sigma_{0}^{-1} + \Sigma^{-1}
\end{aligned}

となる。

複数サンプルの場合

TODO

共分散行列

ウィシャート分布

多変量変量正規分布 $N_{p}(\vec{x}; \vec{0}, \Sigma)$ にしたがう $n$ 個の $p$ 元ベクトル $\{\vec{x}_{i}\}_{i=1}^{n}$ で構成される $p$ 行 $p$ 列の行列

X = \sum_{i=1}^{n} \vec{x}_{i} \vec{x}_{i}^{\top}

が従う確率分布を自由度 $n$ のウィシャート分布という。だたし $n \geq p$ である。

定義は

\mathcal{W}(X; \Sigma, p, n) =
  \frac{|X|^\frac{n-p-1}{2} e^{-\frac{\textrm{tr}(\Sigma^{-1} X)}{2}}}
       {c_{p}(n) |\Sigma|^\frac{n}{2}}

となる。ここで、$c_{p}(n)$ はウィシャート定数であり

c_{p}(n) = 2^\frac{pn}{2} \Gamma_{p}\left(\frac{n}{2}\right)

と定義され、$\Gamma_{p}(\cdot)$ は多変量ガンマ関数であり

\Gamma_{p}(n) = \pi^\frac{p (p-1)}{4}
  \prod_{i=1}^{p} \Gamma\left(\frac{n+(1-i)}{2}\right)

と定義される。

逆ウィシャート分布

ウィシャート分布に従う行列の逆行列が従う分布を、逆ウィシャート分布という。

定義は

\mathcal{W}^{-1}(X; \Sigma, p, n)
  = \frac{|X|^{-\frac{n+p+1}{2}} e^{-\frac{\textrm{tr}(\Sigma X^{-1})}{2}}}{c_{p}(n) |\Sigma|^{-\frac{n}{2}}}

となる。

$X \sim \mathcal{W}(\Sigma, p, n)$ ならば $X^{-1} \sim \mathcal{W}^{-1}(\Sigma^{-1}, p, n)$ となる。

ベイズ推定

尤度関数を多変量正規分布 $\mathcal{N}_{p}(\vec{x}; \vec{\mu}, \Sigma)$ と仮定して、$\Sigma$ の共役事前分布を逆ウィシャート分布 $\mathcal{W}^{-1}(\Sigma, \Sigma_{0}, p, n)$ とする。あたらしくサンプル $\vec{x}$ が得られたときの、$\Sigma_{0}$ の更新式を求めたい。ただし $\vec{\mu}$ は既知とする。

\begin{aligned}
\mathcal{W}^{-1}(\Sigma; \Sigma_{1}, p, n)

  &\propto
    \mathcal{N}_{p}(\vec{x}; \vec{\mu}, \Sigma) \cdot
    \mathcal{W}^{-1}(\Sigma; \Sigma_{0}, p, n)
  \\

  &\propto
    \left(
      |\Sigma|^{-\frac{1}{2}}
      e^{-\frac{1}{2} (\vec{x} - \vec{\mu})^{\top} \Sigma^{-1}
                      (\vec{x} - \vec{\mu})}
    \right)
    \left(
      |\Sigma|^{-\frac{n+p+1}{2}}
        e^{-\frac{\textrm{tr}(\Sigma_{0} \Sigma^{-1})}{2}}
    \right)
  \\

  &= |\Sigma|^{-\frac{(n+1)+p+1}{2}}
    e^{-\frac{1}{2} \left\{
        \textrm{tr}(\Sigma_{0} \Sigma^{-1}) +
        (\vec{x} - \vec{\mu})^{\top} \Sigma^{-1}
        (\vec{x} - \vec{\mu})
      \right\}}
  \\

  &= |\Sigma|^{-\frac{(n+1)+p+1}{2}}
    e^{-\frac{1}{2} \left\{
        \textrm{tr}(\Sigma^{-1} \Sigma_{0}) +
        (\vec{x} - \vec{\mu})^{\top} \Sigma^{-1}
        (\vec{x} - \vec{\mu})
      \right\}}
  \\

  &= |\Sigma|^{-\frac{(n+1)+p+1}{2}}
    e^{-\frac{1}{2} \left[
        \textrm{tr}(\Sigma^{-1} \Sigma_{0}) +
        \textrm{tr}\{
        (\vec{x} - \vec{\mu})^{\top} \Sigma^{-1}
        (\vec{x} - \vec{\mu})\}
      \right]}
  \\

  &= |\Sigma|^{-\frac{(n+1)+p+1}{2}}
    e^{-\frac{1}{2} \left[
        \textrm{tr}(\Sigma^{-1} \Sigma_{0}) +
        \textrm{tr}\{
        \Sigma^{-1}
        (\vec{x} - \vec{\mu}) (\vec{x} - \vec{\mu})^{\top} \}
      \right]}
   \\

  &= |\Sigma|^{-\frac{(n+1)+p+1}{2}}
    e^{-\frac{1}{2} \left[
        \textrm{tr}\{\Sigma^{-1} \left( \Sigma_{0} +
          (\vec{x} - \vec{\mu}) (\vec{x} - \vec{\mu})^{\top}\right)\}
      \right]}
\end{aligned}

より

\Sigma_{1} = \Sigma_{0} + (\vec{x} - \vec{\mu}) (\vec{x} - \vec{\mu})^{\top}

となる (更新後は、自由度 $n+1$ の逆ウィシャート分布になっている?)。

平均ベクトルと共分散行列

TODO

線形回帰の重み

尤度関数に正規分布 $\mathcal{N}(y; \vec{w}^{\top} \vec{x}, \sigma^{2})$ を、重みベクトル $\vec{w}$ の事前分布に多変量正規分布 $\mathcal{N}_{p}(\vec{w}; \vec{\mu}_{0}, \Sigma_{0})$ を仮定して、$\sigma$ を既知とする。

あたらしいサンプル $(\vec{x}, y)$ が得られたときの事後分布は

\begin{aligned}
\mathcal{N}_{p}(\vec{w}; \vec{\mu}_{1}, \Sigma_{1})

  &\propto \mathcal{N}(y; \vec{w}^{\top} \vec{x}, \sigma) \cdot
          \mathcal{N}_{p}(\vec{w}; \vec{\mu}_{0}, \Sigma_{0})
  \\

  &\propto \left( e^{-\frac{1}{2}
                     (\frac{y - \vec{w}^{T} \vec{x}}{\sigma})^{2}} \right)
          \left( e^{-\frac{1}{2}
            (\vec{w} - \vec{\mu}_{0})^{\top} \Sigma_{0}^{-1}
            (\vec{w} - \vec{\mu}_{0}) } \right)
  \\

  &= e^{-\frac{1}{2} \left\{
            (\frac{y - \vec{w}^{T} \vec{x}}{\sigma})^{2} +
            (\vec{w} - \vec{\mu}_{0})^{\top} \Sigma_{0}^{-1}
            (\vec{w} - \vec{\mu}_{0})
       \right\} }
  \\

  &= e^{-\frac{1}{2} \left\{
            \frac{y^{2}}{\sigma^{2}} - 2 y \vec{w}^{\top}
            \frac{\vec{x}}{\sigma^{2}} +
            \frac{(\vec{w}^{\top} \vec{x})^{2}}{\sigma^{2}} +

            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{w} - 2
            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{\mu}_{0} +
            \vec{\mu}_{0}^{\top} \vec{\mu}_{0}
       \right\} }
  \\

  &\propto e^{-\frac{1}{2} \left\{
            - 2 y \vec{w}^{\top}
            \frac{\vec{x}}{\sigma^{2}} +
            \frac{(\vec{w}^{\top} \vec{x})^{2}}{\sigma^{2}} +

            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{w} - 2
            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{\mu}_{0}
       \right\} }
  \\

  &= e^{-\frac{1}{2} \left\{
            - 2 y \vec{w}^{\top} \frac{\vec{x}}{\sigma^{2}} +
            \frac{\vec{w}^{\top} \vec{x} \vec{x}^{\top} \vec{w}}
                 {\sigma^{2}} +

            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{w} - 2
            \vec{w}^{\top} \Sigma_{0}^{-1} \vec{\mu}_{0}
       \right\} }
  \\

  &= e^{-\frac{1}{2} \left\{
            \vec{w}^{\top} (\Sigma_{0}^{-1} +
              \frac{\vec{x} \vec{x}^{\top}}{\sigma^{2}}) \vec{w} - 2
            \vec{w}^{\top} (\Sigma_{0}^{-1} \vec{\mu}_{0}
            + y \frac{\vec{x}}{\sigma^{2}})
       \right\} }
  \\

  &= e^{-\frac{1}{2} \left\{
            \vec{w}^{\top} (\Sigma_{0}^{-1} +
              \frac{\vec{x} \vec{x}^{\top}}{\sigma^{2}}) \vec{w} - 2
            \vec{w}^{\top} (\Sigma_{0}^{-1} \vec{\mu}_{0}
            + y \frac{\vec{x}}{\sigma^{2}})
       \right\} }
  \\

  &= e^{-\frac{1}{2} \left\{
            \vec{w}^{\top}
            (\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}}{\sigma^{2}})
            \vec{w}

            - 2 \vec{w}^{\top}
            (\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}}
                                    {\sigma^{2}})
            (\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}}
                                    {\sigma^{2}})^{-1}
              (
                \Sigma_{0}^{-1} \vec{\mu}_{0} +
                y \frac{\vec{x}}{\sigma^{2}})
       \right\} }
  \\

  &\propto e^{-\frac{1}{2} \left\{
            \left(\vec{w} - (\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}} {\sigma^{2}})^{-1} (\Sigma_{0}^{-1} \vec{\mu}_{0} + y \frac{\vec{x}}{\sigma^{2}})\right)^{\top}

            \left(\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}}{\sigma^{2}}\right)

            \left(\vec{w} - (\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}} {\sigma^{2}})^{-1} (\Sigma_{0}^{-1} \vec{\mu}_{0} + y \frac{\vec{x}}{\sigma^{2}})\right)
       \right\} }
  \\

\end{aligned}

より

\begin{aligned}
\vec{\mu}_{1} &= \left(\Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}} {\sigma^{2}}\right)^{-1} \left(\Sigma_{0}^{-1} \vec{\mu}_{0} + y \frac{\vec{x}}{\sigma^{2}}\right) \\

\Sigma_{1}^{-1} &= \Sigma_{0}^{-1} + \frac{\vec{x} \vec{x}^{\top}}{\sigma^{2}}

\end{aligned}

となる。

実際に計算するには?

$N$ 個のサンプル ${(\vec{x}_{i}, y_{i})}_{i=1}^{N}$ が得られているとする。

\begin{aligned}
A &= \sigma^{2} \Sigma^{-1} \\
b &= \sigma^{2} \Sigma^{-1} \vec{\mu}
\end{aligned}

とおくと、更新式が

\begin{aligned}
A &\leftarrow A + \vec{x}_{i} \vec{x}_{i}^{\top} \\
\vec{b} &\leftarrow \vec{b} + y_{i} \vec{x}_{i}
\end{aligned}

となり、結果として

\begin{aligned}
A &= \sigma^{2} I_{p} + \sum_{i=1}^{N} \vec{x} \vec{x}^{\top} \\
b &= \sum_{i=1}^{N} y_{i} \vec{x}_{i}
\end{aligned}

が得られる。ここで、$\vec{\mu}$ と $\Sigma$ の初期値を、それぞれ零ベクトル $\vec{0}$ と単位行列 $I_{p}$ とした。

線形方程式 $A \vec{\mu} = \vec{b}$ の解として $\vec{\mu}$ の推定値が得られ、$\frac{A}{\sigma^{2}}$ の逆行列として $\Sigma$ が得られる。

指数分布

確率密度関数

\mathcal{E}(x; \lambda) = \lambda e^{-\lambda x}

ベイズ推定

尤度関数を指数分布 $\mathcal{E}(x; \lambda)$ と仮定して、$\lambda$ の共役事前分布をガンマ分布 $\mathcal{G}(\lambda; k, \theta)$ とする。あたらしくサンプル $x$ が得られたときの、$k$ と $\theta$ の更新式を求めたい。

事後分布は

\begin{aligned}
\mathcal{G}(\lambda; k', \theta')

  &= \mathcal{E}(x; \lambda) \cdot
     \mathcal{G}(\lambda; k, \theta) \\

  &\propto \left(\lambda e^{-\lambda x} \right)
     \left(\lambda^{k-1} e^{-\frac{\lambda}{\theta}} \right) \\

  &= \lambda^{(k+1)-1} e^{-\lambda \left(\frac{1}{\theta} + x \right)} \\

\end{aligned}

より

\begin{aligned}
k' &= k + 1 \\
\frac{1}{\theta'} &= \frac{1}{\theta} + x
\end{aligned}

となる。

複数サンプルの場合

$N$ 個のサンプル $\{x_{i}\}_{i=1}^{N}$ が与えられるとする。尤度関数が $\prod_{i=1}^{N} \mathcal{E}(x_{i}; \lambda)$ となる以外は、$1$ 個のサンプルのときと同じになる。

\begin{aligned}
\mathcal{G}(\lambda; k', \theta')

  &= \left( \prod_{i=1}^{N} \mathcal{E}(x_{i}; \lambda) \right)
     \cdot \mathcal{G}(\lambda; k, \theta) \\

  &\propto \left(\lambda^{N} e^{-\lambda \sum_{i=1}^{N} x_{i}} \right)
    \left(\lambda^{k-1} e^{-\frac{\lambda}{\theta}} \right) \\

  &= \lambda^{(k+N)-1} e^{-\lambda \left(\frac{1}{\theta} + \sum_{i=1}^{N} x_{i} \right)}
\end{aligned}

より

\begin{aligned}
k' &= k + N \\
\frac{1}{\theta'} &= \frac{1}{\theta} + \sum_{i=1}^{N} x_{i}
\end{aligned}

となる。

ポアソン分布

確率密度関数

\mathcal{P}(x; \lambda) = \frac{\lambda^{x} e^{-\lambda}}{\Gamma(x+1)}

ベイズ推定

尤度関数をポアソン分布 $\mathcal{P}(x; \lambda)$ と仮定して、$\lambda$ の共役事前分布をガンマ分布 $\mathcal{G}(\lambda; k, \theta)$ とする。あたらしくサンプル $x$ が得られたときの、$k$ と $\theta$ の更新式を求めたい。

事後分布は

\begin{aligned}
\mathcal{G}(\lambda; k', \theta')

  &= \mathcal{P}(x; \lambda) \cdot
     \mathcal{G}(\lambda; k, \theta) \\

  &\propto \left(\lambda^{x} e^{-\lambda} \right)
     \left(\lambda^{k-1} e^{-\frac{\lambda}{\theta}} \right) \\

  &= \lambda^{(k+x)-1} e^{-\lambda \left(\frac{1}{\theta} + 1 \right)} \\

\end{aligned}

より

\begin{aligned}
k' &= k + x \\
\frac{1}{\theta'} &= \frac{1}{\theta} + 1
\end{aligned}

となる。

複数サンプルの場合

$N$ 個のサンプル $\{x_{i}\}_{i=1}^{N}$ が与えられるとする。尤度関数が $\prod_{i=1}^{N} \mathcal{P}(x_{i}; \lambda)$ となる以外は、$1$ 個のサンプルのときと同じになる。

\begin{aligned}
\mathcal{G}(\lambda; k', \theta')

  &= \left( \prod_{i=1}^{N} \mathcal{P}(x_{i}; \lambda) \right)
     \cdot \mathcal{G}(\lambda; k, \theta) \\

  &\propto \left(\lambda^{\sum_{i=1}^{N} x_{i}} e^{-\lambda N} \right)
    \left(\lambda^{k-1} e^{-\frac{\lambda}{\theta}} \right) \\

  &= \lambda^{\left( k+\sum_{i=1}^{N} x_{i} \right)-1} e^{-\lambda \left(\frac{1}{\theta} + N \right)}
\end{aligned}

より

\begin{aligned}
k' &= k + \sum_{i=1}^{N} x_{i} \\
\frac{1}{\theta'} &= \frac{1}{\theta} + N
\end{aligned}

となる。

二項分布

確率密度関数

\mathcal{B}(k; n, p) =
  \left(\begin{array}{c} n \\ k \end{array}\right)
    p^{k} (1 - p)^{n-k}

$n = 1$ のときベルヌーイ分布になる。

ベータ分布

\mathcal{Beta}(x; \alpha, \beta) =
  \frac{x^{\alpha-1} (1 - x)^{\beta - 1}}{B(\alpha, \beta)}

ベイズ推定

尤度関数を二項分布 $\mathcal{B}(k; n, p)$ と仮定して、$p$ の共役事前分布をベータ分布 $\mathcal{Beta}(p; \alpha, \beta)$ とする。あたらしくサンプル $k$ が得られたときの、$p$ の更新式を求めたい。

事後分布は

\begin{aligned}
\mathcal{Beta}(x; \alpha', \beta')
  &\propto \mathcal{B}(k; n, p) \cdot
           \mathcal{Beta}(p; \alpha, \beta) \\

  &\propto \left( p^{k} (1 - p)^{n-k} \right)
    \left( p^{\alpha-1} (1 - p)^{\beta - 1} \right) \\

  &= p^{\alpha + k - 1} (1 - p)^{\beta + (n - k) - 1}
\end{aligned}

より

\begin{aligned}
\alpha' &= \alpha + k \\
\beta' &= \beta + n - k
\end{aligned}

となる。

多項分布

確率密度関数

\mathcal{M}\left(\{x_{k}\}_{k=1}^{K}; \{p_{k}\}_{k=1}^{K} \right) =
  \frac{\Gamma\left(\left(\sum_{k=1}^{K} x_{k} \right) + 1\right)}
       {\prod_{k=1}^{K} \Gamma(x_{k}+1)}
  \prod_{k=1}^{K} p_{k}^{x_{k}}

ディリクレ分布

\mathcal{D}\left(\{x_{k}\}_{k=1}^{K}; \{\alpha_{k}\}_{k=1}^{K} \right) =
  \frac{\prod_{k=1}^{K} x_{k}^{\alpha_{k} - 1}}
       {B\left(\{\alpha_{k}\}_{k=1}^{K} \right)}

ここで $B(\cdot)$ は多変量ベータ関数であり

B\left(\{\alpha_{k}\}_{k=1}^{K} \right) =
  \frac{\prod_{k=1}^{K} \Gamma(\alpha_{k})}
       {\Gamma\left(\sum_{k=1}^{K} \alpha_{k} \right)}

と定義される。

ベイズ推定

尤度関数を多項分布 $\mathcal{M}\left(\{x_{k}\}_{k=1}^{K}; \{p_{k}\}_{k=1}^{K} \right)$ と仮定して、$\{p_{k}\}_{k=1}^{K}$ の共役事前分布をディリクレ分布 $\mathcal{D}\left(\{p_{k}\}_{k=1}^{K}; \{\alpha_{k}\}_{k=1}^{K} \right)$ とする。あたらしくサンプル $\{x_{k}\}_{k=1}^{K}$ が得られたときの、$p$ の更新式を求めたい。

事後分布は

\begin{aligned}
\mathcal{D}\left(\{p_{k}\}_{k=1}^{K}; \{\alpha'_{k}\}_{k=1}^{K} \right)
  &\propto
    \mathcal{M}\left(\{x_{k}\}_{k=1}^{K}; \{p_{k}\}_{k=1}^{K} \right)
    \cdot \mathcal{D}\left(\{p_{k}\}_{k=1}^{K}; \{\alpha_{k}\}_{k=1}^{K} \right)
  \\

  &\propto \left( \prod_{k=1}^{K} p_{k}^{x_{k}} \right)
    \left( \prod_{k=1}^{K} p_{k}^{\alpha_{k} - 1} \right) \\

  &= \prod_{k=1}^{K} p_{k}^{\alpha_{k} + x_{k} - 1}
\end{aligned}

より

\alpha'_{k} = \alpha_{k} + x_{k}

となる。