Bayes推定

はじめに

先の投稿でBayes推定に触れた。今回は、回帰問題を例に取り、最尤推定・MAP推定・Bayes推定を統一的に解説したい。

尤度・事後分布・事前分布

最初に先の解説の一部を再掲載する。
いま、X=\{\vec{x}_1,\cdots,\vec{x}_N\},Y=\{y_1,\cdots,y_N\}が観測されているとする。\vec{x}_iが説明変数、y_iが目的変数である。このとき、潜在変数\thetaを導入し、同時確率分布p(X,Y,\theta)を考える。この分布にBayesの定理を適用すると次式を得る。

(1)   \begin{equation*} p(\theta|X,Y) = \frac{p(Y|X,\theta)p(\theta)}{p(Y|X)} \end{equation*}

ただし、式変形の途中で、p(X)=p(X|\theta)を用いた。式(1)の左辺にあるp(\theta|X,Y)を事後分布、右辺の分子にあるp(Y|X,\theta)を尤度、p(\theta)を事前分布、右辺分母にあるp(Y|X)をモデルエビデンスと呼ぶ。いま、事後分布の\theta依存性に注目すると、式(1)の右辺分母は定数と見なすことができる。したがって

(2)   \begin{equation*} p(\theta|X,Y) \propto p(Y|X,\theta)p(\theta) \end{equation*}

と書くことができる。尤度と事前分布に、\thetaをパラメータとして持つ適当な関数を仮定し回帰問題を解くことになる。

最尤推定

最尤推定では式(2)の右辺にある尤度p(Y|X,\theta)を潜在変数\thetaについて最大化することにより回帰問題を解く。事前分布は考慮しない。観測値が独立同分布から得られていると仮定すると、尤度は以下のように変形できる。

(3)   \begin{equation*} p(Y|X,\theta)=\prod_{n=1}^N\;p(y_n|\vec{x}_n,\theta) \end{equation*}

\vec{x}_nをM次元のベクトルとし、尤度p(y_n|\vec{x}_n,\theta)として次の正規分布を仮定する。

(4)   \begin{equation*} p(y_n|\vec{x}_n,\theta)=\mathcal{N}(y_n|\vec{w}^{\;T}\vec{x}_n,\lambda^{-1}) \end{equation*}

ここで、M次元ベクトル\vec{w}と\lambdaがパラメータ\thetaに相当する。\lambdaは分散の逆数であり、精度（precision）と呼ばれる量である。式(3)の対数をとり、式(4)を代入すると次式を得る。

(5)   \begin{eqnarray*} \ln{p(Y|X,\theta)}&=&\sum_{n=1}^N\ln{\left\{\left(\frac{1}{2\pi}\right)^{1/2}\lambda^{1/2}\exp{\left[-\frac{\lambda}{2}\left(y_n-\vec{w}^{\;T}\vec{x}_n\right)^2\right]}\right\}}\\ &=&-\frac{N}{2}\ln{2\pi}+\frac{N}{2}\ln{\lambda}-\frac{\lambda}{2}\sum_{n=1}^{N}(y_n-\vec{w}^{\;T}\vec{x}_n)^2 \end{eqnarray*}

式(5)を\vec{w}で微分した式を0とおけば、最尤推定による\vec{w}（\vec{w}_{ML}と書くことにする）が決定される。

(6)   \begin{equation*} \vec{w}_{ML}=\left(\Phi^T\Phi\right)^{-1}\Phi^T\vec{y} \end{equation*}

ここで、\PhiはN\times Mの行列であり、その成分は\Phi_{ni}=x_{ni}で定義される。x_{ni}はベクトル\vec{x}_nの第i成分である。また、\vec{y}=(y_1,\cdots,y_N)^Tとした。同様に、\lambdaで微分した式を0とおけば、\lambdaが求まる。

(7)   \begin{equation*} \frac{1}{\lambda}_{ML}=\frac{1}{N}\sum_{n=1}^{N}(y_n-\vec{w}_{ML}^{\;T}\vec{x}_n)^2 \end{equation*}

未知変数\vec{x}_{*}が与えられた時、最尤推定が予測する目的変数の値y_{*}は次式で求められる。

(8)   \begin{equation*} y_*=\vec{w}_{ML}^{\;T}\vec{x}_* \end{equation*}

その分散は

(9)   \begin{equation*} \sigma_{ML}^2=\frac{1}{\lambda}_{ML} \end{equation*}

である。分散は定数であり、未知変数\vec{x}_*に依存しない。最尤推定ではパラメータ\thetaの取り得る範囲に制限を与えず最適な値を探索する。これが過学習をもたらす原因である。\thetaの取り得る範囲に制限を付加したものが次に示すMAP推定である。

MAP推定

MAP推定では事前分布p(\theta)を考慮することで、\thetaの取り得る範囲に制限を与える。

(10)   \begin{equation*} p(Y|X,\theta)p(\theta) = \prod_{n=1}^{N}p(y_n|\vec{x}_n,\theta)p(\theta) \end{equation*}

この式を最大化することで回帰問題を解く。p(y_n|\vec{x}_n,\theta)として先と同じ正規分布を仮定し、p(\theta)として次の正規分布を仮定する。

(11)   \begin{equation*} p(\theta)=\mathcal{N}(\vec{w}|\vec{m},\Lambda^{-1}) \end{equation*}

ここで、\vec{m}はM次元ベクトル、\LambdaはM\times M行列である。ともにあらかじめ決めておく量である。(\vec{w},\lambda)は最尤推定と同じく計算により求める量であることに注意する。式(10)に式(4)と式(11)を代入し、対数を取る。

(12)   \begin{eqnarray*} \ln{p(Y|X,\theta)p(\theta)} = &-& \frac{N}{2}\ln{2\pi}+\frac{N}{2}\ln{\lambda}-\frac{\lambda}{2}\sum_{n=1}^{N}(y_n-\vec{w}^{\;T}\vec{x}_n)^2 \\ &-& \frac{M}{2}\ln{2\pi}+\frac{1}{2}\ln{\det{\Lambda}}-\frac{1}{2}\left(\vec{w}-\vec{m}\right)^T\Lambda\left(\vec{w}-\vec{m}\right) \end{eqnarray*}

この式を\vec{w}で微分した式を0とおくと、MAP推定による\vec{w}（\vec{w}_{MAP}と書くことにする）を得る。

(13)   \begin{equation*} \vec{w}_{MAP}=\left(\lambda\Phi^T\Phi+\Lambda\right)^{-1}\left(\lambda\Phi^T\vec{y}+\Lambda\vec{m}\right) \end{equation*}

\lambdaで微分すると次式を得る。

(14)   \begin{equation*} \frac{1}{\lambda}_{MAP}=\frac{1}{N}\sum_{n=1}^{N}(y_n-\vec{w}_{MAP}^{\;T}\vec{x}_n)^2 \end{equation*}

いま簡単のため、\vec{m}=\vec{0}, \Lambda=\alpha Iとする。すなわち、\vec{w}は正規分布\mathcal{N}(\vec{w}|\vec{0},\alpha^{-1}I)から生成されると仮定する。このとき

(15)   \begin{eqnarray*} \vec{w}_{MAP}&=&\lambda S\Phi^T\vec{y} \\ S&\equiv&\left(\lambda\Phi^T\Phi+\alpha I\right)^{-1} \end{eqnarray*}

となる。未知変数\vec{x}_{*}が与えられた時、MAP推定が予測する目的変数の値y_{*}は次式で求められる。

(16)   \begin{equation*} y_*=\vec{w}_{MAP}^{\;T}\vec{x}_* \end{equation*}

その分散は

(17)   \begin{equation*} \sigma_{MAP}^2=\frac{1}{\lambda}_{MAP} \end{equation*}

である。分散は定数であり、未知変数\vec{x}_*に依存しない。

最尤推定・MAP推定の実験

最尤推定とMAP推定を用いて、具体的に回帰問題を解いてみる。考える正解曲線f_{G}(x)を次式で定義する。

(18)   \begin{equation*} f_G(x)&=&x+\sin{3x} \end{equation*}

適当に10(=N)個の点をサンプルしこれを観測点とする。観測点を生成する際に、以下のようにノイズを付加した。

(19)   \begin{eqnarray*} y&=&f_G(x)+\epsilon\\ \epsilon&\sim&\mathcal{N}(\epsilon|0,\sigma^2) \end{eqnarray*}

\sigma=0.015とした。また、\vec{x}=(1,x,x^2,\cdots,x^{M-1})とおく。すなわち、M-1次式で曲線を再現することを考える。最初に最尤推定の結果を示す。

次に、MAP推定の結果を示す。\alpha=0.1とした。

Mを大きくしていくと、最尤推定の方は、左側にある観測点で過学習となり始め、右側の観測点を再現できなくなる。一方、MAP推定の方はどの次元においても観測点近傍を通る曲線を再現できている。各図のタイトルの横に標準偏差を示した。最尤推定・MAP推定から求まる分散は\vec{x}に依存しない定数なので、各予測点の確からしさを表す量ではない。最尤推定とMAP推定は、「点」を予測する手法であり、点推定（Point Estimation）と総称される。

Bayes推定

再度、式(1)に戻って議論を行う。

(20)   \begin{equation*} p(\theta|X,Y) = \frac{p(Y|X,\theta)p(\theta)}{p(Y|X)} \end{equation*}

分母にあるモデルエビデンスは、分子を\thetaで積分したものである。

(21)   \begin{equation*} p(Y|X)=\int d\theta\;p(Y|X,\theta)p(\theta) \end{equation*}

事後分布p(\theta|X,Y) を求めることができれば、次式により未知変数\vec{x}_*に対する予測値y_*を計算することができる。

(22)   \begin{equation*} p(y_*|\vec{x}_*,X,Y)=\int d\theta\;p(y_*|\vec{x}_*,\theta)p(\theta|X,Y) \end{equation*}

Bayes推定の目的は、事後分布p(\theta|X,Y) を求めることである。
尤度と事前分布として、最尤推定・MAP推定のときと同じ正規分布を仮定する。尤度に対しては式(4)を、事前分布に対しては式(11)を用いる。ここで、注意しなければならないことがある。MAP推定では天下り的に事前分布として\vec{w}に対する確率分布だけを与えたが、Bayes推定では、厳密に以下の場合を考える必要がある。

  1. \lambdaを固定し、\vec{w}を推定したい場合
  2. \vec{w}を固定し、\lambdaを推定したい場合
  3. 両方とも推定したい場合

それぞれに対し、任意の確率分布を仮定することができるが、次のようにおくことにより事後分布を解析的に扱いやすい形にすることができる。

  • 上記1の場合、事前分布を\vec{w}に対する正規分布とおくと、事後分布も\vec{w}に対する正規分布となる。
  • 上記2の場合、事前分布を\lambdaに対するガンマ分布とおくと、事後分布も\lambdaに対するガンマ分布となる。
  • 上記3の場合、事前分布を\vec{w},\lambdaに対するガウス・ガンマ分布とおくと、事後分布も\vec{w},\lambdaに対するガウス・ガンマ分布となる。

    このように事後分布を同じ関数形にできる事前分布を共役事前分布と呼ぶ。ここでは、1を採用したことになる。すなわち、\lambdaについてはあらかじめ決めておく定数とする。
    さて、これらを式(20)の右辺分子に代入する。分母の計算には式(21)を用いる。ガウス関数なので積分は容易に実行できる。最終的に事後分布も正規分布となる。

    (23)   \begin{equation*} p(\theta|X,Y)=\mathcal{N}(\vec{w}|\vec{m}_B,\Lambda_B^{-1}) \end{equation*}

    ここで、

    (24)   \begin{eqnarray*} \Lambda_B=\lambda\sum_{n=1}^N\vec{x}_n\vec{x}_n^T+\Lambda \\ \vec{m}_B=\Lambda_B^{-1}\left(\lambda\sum_{n=1}^{N}y_n\vec{x}_n+\Lambda\vec{m}\right) \end{eqnarray*}

    とした。さらに、式(22)も解析的に計算することができる。

    (25)   \begin{equation*} p(y_*|\vec{x}_*,X,Y)=\mathcal{N}(y_*|\mu_*,\lambda_*^{-1}) \end{equation*}

    ただし、

    (26)   \begin{eqnarray*} \mu_*&=&\vec{m}^{T}_B\vec{x}_* \\ \lambda_*^{-1}&=&\lambda^{-1}+\vec{x}^T_*\Lambda^{-1}_B\vec{x}_* \end{eqnarray*}

    とおいた。MAP推定のときと同じ仮定

    (27)   \begin{eqnarray*} \vec{m}&=&\vec{0}\\ \Lambda&=&\alpha I \end{eqnarray*}

    をおくと

    (28)   \begin{eqnarray*} \mu_*&=&\vec{w}_{MAP}^{\;T}\vec{x}_* \\ \vec{w}_{MAP}&=&\lambda S\Phi^T\vec{y} \\ \lambda_*^{-1}&=&\lambda^{-1}+\vec{x}^T_*S\vec{x}_*\\ S&=&\left(\lambda\Phi^T\Phi+\alpha I\right)^{-1} \end{eqnarray*}

    を得る。MAP推定のときに求めた\vec{w}_{MAP}(式(15))が現れることに注意する。Bayes推定では、式(25)の正規分布が求めるべき解となる。すなわち、回帰曲線を「点」ではなく確率として求める点が最尤推定やMAP推定との大きな違いである。また、\lambda,\alphaはあらかじめ決めておく量であることに注意する。
    先と同じ問題に適用した結果を以下に示す。今の場合、正解曲線からサンプルする際のノイズ分布は既知（平均0、標準偏差\sigma=0.015の正規分布）であるから、\lambda=1/\sigma^2とすることができるが、現実の問題に適応する際は、観測データについての何らかの事前知識が必要になる。\alphaについてはMAP推定のときと同じ値（0.1）を用いた。

    図の点線は正規分布（式(25)）の平均値をプロットしたものである。色ぬりした領域は標準偏差である。最尤推定やMAP推定との大きな違いは標準偏差が各点ごとに計算されることである。これがBayes推定から求まる予測値の不確かさ（uncertainty）に相当する。不確かさは、観測値が密にある領域では小さく、疎である領域では大きくなることが分かる。

    まとめ

    今回は、最尤推定、MAP推定、Bayes推定に対する統一的な解説を試みた。最尤推定は尤度を最大化するパラメータを無制限に探索する。MAP推定では事前分布を導入し、パラメータの探索範囲を制限することにより、過学習を抑制する。これら2つの手法は点推定と呼ばれ、予測値の不確かさを求めることはできない。一方、Bayes推定では、点ではなく確率（事後分布）を求めることが目的となる。確率を求めるのであるから、予測値の不確かさを自然に表現することができる。

    コードの場所

    今回のコードはここにある。pythonではなくてjuliaで実装した。

    • Kumada Seiya

    Kumada Seiya

    仕事であろうとなかろうと勉強し続ける、その結果”中身”を知ったエンジニアになれる

    最近の記事

    • 関連記事
    • おすすめ記事
    • 特集記事

    アーカイブ

    カテゴリー

    PAGE TOP