logo logo

March 20, 2019 10:04

PRMLの表紙を再現しながら学ぶベイズ線形回帰

ベイズ線形回帰の理論フロー

統計的機械学習は、ある真の確率分布から発生するサンプルの集合を観測し、元の真の確率分布を推定することを目標にします。以下は回帰問題を想定した表記になります。(予測分布$p(y | \bf{X},\bf{Y},\bf{x})$と呼ばれる。ここで$y$は回帰変数、$\bf{X}$、$\bf{Y}$はデータサンプルの特徴量と対応する値を表しています.
$\bf{x}$は入力特徴量です)古典的な線形モデルから最新の複雑な深層学習モデルまで、応用のパターンも多岐にわたりますが、ベースにある目標は共通しています。

Step1:2つの仮定モデル

さて、ベイズ学習ではこの真の確率分布を再現するため、パラメータがデータ観測を元に更新(学習される)されるようなある生成モデル $p(y|\bf{x},\bf{w})$を1つ仮定します。$\bf{w}$はこの生成モデルのパラメータです。さらにベイズ学習ではパラメータについてもある分布$p(\bf{w})$を仮定します。これを事前分布というのでした。これはデータが何も得られていないときに、何らかの尤もらしい仮定を事前分布を通して導入することに相当します。

Step.2 : パラメータについての事後分布の計算

次にパラメータ$\bf{w}$として尤もらしい値を、得られたデータを元に更新していきます。こうして得られる分布は事後分布と呼ばれ、$p(\bf{w}|\bf{X},\bf{Y})$と表記します。いわゆるベイジアンアップデートと呼ばれる過程です。この過程でベイズの定理が使われます。事後分布と事前分布を結ぶ重要な役割をするのが尤度関数で、$p(\bf{X},\bf{Y}|\bf{w})$と表記します。多くの場合仮定される独立同分布の下では、生成モデルの積として書くことができます。尤度関数を最大化する最尤推定はよく使われますが、ベイズ線形回帰の場合は以下のベイズの定理から従う式
$$p(\bf{w} | \bf{X}, \bf{Y}) = \frac{p(\bf{Y}| \bf{X}, \bf{w})p(\bf{w})}{\int dw p(\bf{Y} | \bf{X}, w)p(\bf{w})},$$
の分母も計算してあげないといけません。ちなみに事後分布を最大化するパラメータを探るMAP推定では、$\bf{w}$についての最大化が目的なので、$\bf{w}$に依存しない分母を計算する必要はありません。一般的には解析的に計算するにも数値的に計算するにもヘビーな計算になることが多いです。エビデンス周辺尤度と呼ばれます。

Step.3 : 予測分布の計算

こうしてデータを元に更新された$\bf{w}$のよりまともであると考えられる分布が得られたわけですが、最後に予測分布を得るには生成モデルと事後分布を加法定理によって結び付けます。もう一度$\bf{w}$についての積分を実行する必要があります。
$$p(y | \bf{X},\bf{Y}, \bf{x}) = \int dw p(y | \bf{w},\bf{x})p(\bf{w}|\bf{X},\bf{Y})$$

こうして新規の入力に対して、予測される$y$の平均や分散を議論することができ、予測の信頼性の議論もしやすくなります。以下ではガウス線形回帰モデルの仮定の下、予測分布がサンプルを集めていく毎に収束していく様子をシミュレーションしながら、PRMLの表紙を再現します。

ベイズ線形回帰

線形回帰モデルとは、入力変数$x$を変形する基底関数$\phi(x)$の線形結合によって、目標値$y$を予測するようなモデルの総称です。数式で書くと
$$y = \bf{w}^{T}\bf{\phi}(\bf{x}),$$
となります。初学者向けのよくある注意事項として、$\bf{\phi}(\bf{x})$は一般的に非線形関数でよく、線形であるのはあくまでパラメータベクトルと基底関数との関係です。そのため非線形変換を含んだ複雑な分布への対応力を残しつつ、線形な箇所についての解析の容易性をバランスよく含んだモデルであり、応用的にも理論的にも重要なモデルとされています。

問題設定

sin関数(青線)に正規ノイズがのった以下のようなサンプル(オレンジの点)に対する回帰を目標とします。





ガウス基底関数

まずは基底関数の形を定める。ここでは空間的に局所化したガウス基底関数
$$\phi _{j}(x) = \exp\left[-\frac{(x - \mu _{j})^2}{2s^{2}} \right]$$
として、$\mu$の値は0.1から0.1刻みで0.9まで、分散$s=0.1$とする。
また定数関数$\phi _{0}(x) = 1$をこの基底群に加えておきます。

基底関数をプロットした図が以下です。

理論との対応(Step1)

基底関数が決まったので、生成モデルの形を決めます。今回は精度$\beta$の正規ノイズが加わった形で定義します。なお、この$\beta$はモデル化するものがえいやで決めるハイパーパラメータです。
$$p(y | \bf{x},\bf{w}) = \mathcal{N}(y | \bf{w}^{T}\bf{\phi(\bf{x})},\beta^{-1})$$
さらにデータが何もないときの事前分布の形を仮定します。計算の便利さから、ここでもガウス分布を仮定しましょう(精度ハイパーパラメータを$\alpha$とします)。尤度関数との積もガウス分布となり、事後分布と事前分布の形が同じになるので、このような分布を共役事前分布と呼ぶのでした。
$$p(\bf{w}) = \mathcal{N}(\bf{w} | \bf{0}, \alpha^{-1}\bf{I})$$

理論との対応(Step2)

独立同分布の仮定により、生成モデルから尤度関数が次のように書けます。
$$p(\bf{Y} | \bf{X}, \bf{w}) = \prod _{n=1}^{N}p(y|\bf{x},\bf{w})$$
これで役者がそろったので、事後分布が計算できます。登場人物が全て正規分布であることから、この計算は解析的に行うことができ(とはいえ確かめるのは大変なのでPRML本書2章などを参照)、
$$p(\bf{w} | \bf{X}, \bf{Y}) = \mathcal{N}(\bf{w} | \bf{m _{N}}, \bf{S} _{N}),$$
$$\bf{m _{N}} = \bf{S} _{N}\left(\bf{S} _{0}^{-1}\bf{m} _{0} + \beta\bf{\Phi}^{T}\bf{t}\right),$$
$$\bf{S} _{N}^{-1} + \beta\Phi^{T}\Phi$$
となります。ただし事前分布の分散行列を$\bf{S} _{0} = \alpha^{-1}\bf{I}$としており、$\Phi$は計画行列とよばれるもので、その$ij$成分が$\phi _{j}(\bf{x} _{i})$で与えられる行列となります。

理論との対応(Step3)

最後に予測分布の計算を上にある周辺化によって行います。結果だけ記すことにすると
$$p(y | \bf{X}, \bf{Y}, \bf{x}) = \mathcal{N}(y | \bf{m} _{N}^{T}\bf{\phi}(\bf{x}),\sigma _{N}^{2}(\bf{x})),$$
$$\sigma _{N}^{2}(\bf{x}) = \frac{1}{\beta} + \bf{\phi}(\bf{x})^{T}\bf{S} _{N}\bf{\phi}(\bf{x})$$

実験結果

この予測分布をサンプル点を増やす毎に数値計算し、グラフ上にプロットしてみました。
赤オレンジの点がノイズの乗ったサンプル点で、青線が復元したいsin関数です。オレンジ線は予測分布の平均(上のしきにおける$\bf{m} _{N}^{T}\bf{\phi}(\bf{x})$)を表しており、サンプル点が増える毎に目標値に近づいている様子がわかります。また平均±分散のレンジを示したのが水色の影になります。サンプル点が少ない時に、これは予測の自信度として解釈できます。データが少ない時にある予測は今まで見てきたパターンに近いから自信度は高いが、見たことないものについては自信がないといった解釈が容易にできます。

n = 1

n = 2

n = 4

n = 16

n = 50

参照