ベイズ統計を理解する（2）〜判別モデルと生成モデル〜

先週に引き続き、Computer Vision Modelsのテキストとベイズ推定とグラフィカルモデル：コンピュータビジョン基礎1 | Udemyを使用して、ベイズ統計について勉強中です。今週は、判別モデルと生成モデル、それぞれのモデリング方法での基礎的な学習方法を、回帰と識別それぞれのケースでどのように行うかを書いていきたいと思います。

ベイズ推論とモデリング

まずベイズを用いた学習では、観測変数xが与えられた時に、その時の状態wを求めるためのモデル $Pr(w|x)$ を学習により獲得することを目標とします。この事後分布 $Pr(w|x)$ を利用するプロセスは、

xとwを関連付けるモデルをパラメーターθを用いて表す
学習データからパラメーターθを学習する
与えられた観測データxから $Pr(w|x,\theta)$ を計算する。その計算が難しい場合は、サンプル結果やMAP推定値を返す

といったようになります。

1.で使用するモデリングの方法には、 $Pr(w|x)$ を直接モデル化しパラメーターを学習する方法（判別モデル：Discriminative Model）と、 $Pr(x|w)$ をモデル化しパラメーターを学習する方法（生成モデル：Generative Model）の2つがあります。

判別モデルでは直接推論を行うことができるので簡易でありまた計算コストも抑えることができます。一方、生成モデルはデータセットをサンプリングでき、また推論時にwの事前分布を使用できるため、ノイズの多いデータでは判別モデルより頑健性の高い結果を得ることができます。

それぞれのモデルで、回帰と分類で行われるモデリングの例を見ていきましょう。

まず、回帰の例です。判別モデルでは

${ Pr(w|x,\theta) = Norm_w(\phi_0 + \phi_1 x, {\sigma}^2) }$

と表現することができます。この場合、 $\sigma^2$ が0であれば、 $w = \phi_0 + \phi_1 x$ となることからわかるように、線形回帰モデルとしてxとwの関係が表現されています。もちろん、 $Pr(w|x)$ が正規分布であると仮定する必要もありませんし、wとxの関係が線形なものである必要も、分散に影響を与えないようなモデルにする必要もありません。モデル化後は、お好きな推定手法を用いてパラメーターθを推定すればOKです。

一方、生成モデルでは

${ Pr(x|w,\theta) = Norm_x(\phi_0 + \phi_1 w, {\sigma}^2) }$

と表現することができます。見ていただいたように、xとwが反転しているだけで、これも線形回帰モデルとなります。この場合も基本的にモデリング自体はもっと自由に行うことができます。生成モデルが判別モデルと異なるところは、事前分布を仮定してやる必要がある点です。この場合、事前分布 $Pr(w)$ を正規分布として仮定しており、必要に応じてパラメーターを学習する必要があります。

次に、分類の例です。判別モデルでは、

${ Pr(w|x,\theta) = Bern_w(sig(\phi_0 + \phi_1 x)) }$

と表現することができます。この場合、Pr(w)をベルヌーイ分布で仮定し、wとxの関係をシグモイド関数で表現していることから、ロジスティック回帰となります。もちろん、分布や関係式は自由に設定することができます。

一方、生成モデルでは大きく異なり、

${ Pr(x|w,\theta) = Norm_x(\mu_w, {\sigma}^2_w) }$

と表現されます。これは、wの値ごとにxの分布が変化することを確率密度関数として表しています。この時、 $Pr(w)$ はベルヌーイ分布で表現し、このパラメーターも必要に応じて学習する必要があります。

生成モデルによる分類

生成モデルによる分類は、確率密度関数の推定という形で行われます。

例えば、 $Pr(x|w) = Norm_w(\mu_w, \Sigma_w)$ とwの値ごとの確率変数xの確率密度関数を定義します。この時最尤推定を行うことにより、推定値を

$\hat{\mu}_w, \hat{\Sigma}_w = argmax(\prod_{i \in S_w} Pr(x_i | \mu_w, \Sigma_w))$

と求めることができます。 $S_w$ はカテゴリがwである観測データの集合であるとします。つまり二値分類でいうと、w=0の場合のパラメーターとw=1の場合のパラメーターはそれぞれのカテゴリに分類される観測データxから別々に推定されることになります。

最後に評価を行う場合は、 $Pr(w) = Bern_w(\lambda)$ を利用して、 $Pr(w=0|x)$ と $Pr(w=1|x)$ をベイズルールで求め利用します。ちなみに $\lambda$ の値も学習データを利用して、何かしらの方法（例えば頻度とか）で決めておく必要がありますが、以降でもあまり話題には上がりませんので「何か決まってるんだなー」くらいに思っておいてください。

以上で生成モデルによる分類を行う基本的な方法を説明しましたが、尤度が正規分布で表現されるため、この場合認識精度はそこまでよくはなりません。精度を上げるには一般的に尤度に複雑な分布を利用する必要がありますが、その場合、通常の方法では推定値を求めることができません。

そこで利用するのが、隠れ変数とEMアルゴリズムです。

隠れ変数の考え方は非常にシンプルです。尤度関数を

$Pr(x|\theta) = \int Pr(x,h|\theta)$ dh

と隠れ変数hをおいて定義します。この時最尤推定を行うことにより、以下のように推定値を求めます。

$\hat{\theta} = argmax(\sum_{i=1}^I log(\int Pr(x_i,h_i|\theta) dh_i ))$

このように隠れ変数を設定すると、単純なモデルを利用してモデルの表現力を上げることができます。例えば、以下で説明する正規混合分布では、h=1の時の分布、h=2の時の分布といったように複数の正規分布を利用し、それを足し合わせることにより複雑な分布を表現することができています。

しかし、上記のような積分の入った式では解析的に簡単にパラメーターの推定を行うことができません。そこで利用されるのがEMアルゴリズムと呼ばれる手法です。EMアルゴリズムでは、下界（Lower Bound）の考え方を利用して、パラメーターを段階的に更新していきます。この手法で一番重要な式は、下記の下界の定義式です。

${LB(q_i(h_i), \theta) = \sum_{i=1}^I \int q_i(h_i) log(\frac{Pr(x_i,h_i|\theta)} {q_i(h_i)}) dh_i ) \le \sum_{i=1}^I log(\int Pr(x_i,h_i|\theta) dh_i )}$

ちなみに、 $\int q_i(h_i) log(\frac{Pr(x_i,h_i|\theta)} {q_i(h_i)}) dh_i$ はKLダイバージェンスと呼ばれるもので、確率分布の非類似度を表す情報量の一つです。EMアルゴリズムはこの下界を最大化する分布 $q_i(h_i)$ とパラメーターθを交互に求めます。Eステップでは、

$\hat{q}_i(h_i) = Pr(h_i|x_i,\theta^t) = \frac{\Pr(x_i|h_i,\theta^t)Pr(h_i|\theta^t)} {Pr(x_i)}$

と $\theta^t$ を固定して $\hat{q}_i(h_i)$ を求めます。Mステップでは、

$\hat{\theta}^{t+1} = argmax(\sum_{i=1}^I \int \hat{q}_i(h_i) log(Pr(x_i,h_i|\theta)) dh_i ))$

のようにパラメーターの推定値を求めます。この計算を繰り返し、パラメーターをより元の尤度関数を最大化するパラメーターに近づけていきます。この手順は、Eステップで下界の関数を更新し、Mステップでその関数を最大にするパラメーターを求めていることに相当します。

では、混合正規分布を例に、EMアルゴリズムの手順を見ていきましょう。混合正規分布は、複数の正規分布の和で確率分布を表現するものとなります。求める確率分布は

$Pr(x|\theta) = \sum_{k=1}^K \lambda_k Norm_x(\mu_k,\Sigma_k)$

と定義できます。 $\lambda$ は合計すると1になる値です。この時、

$Pr(x|h,\theta) = Norm_x(mu_h,\Sigma_h)$

$Pr(h|\theta) = Cat_h(\lambda)$

と定義することができる。Eステップは、

$\hat{q}_i(h_i) = \frac{\lambda_k Norm_{x_i}(\mu_k,\Sigma_k)} {\sum_{j=1}^K \lambda_j Norm_{x_i}(\mu_j,\Sigma_j)} = r_{ik}$

となる。 $r_{ik}$ は負担率とも呼ばれる値で、観測データ $x_i$ をどの分布で負担するかを表すような値になっています。Mステップでは、この負担率を用いて以下のように更新する。

$\lambda_k^{t+1} = \frac{\sum_{i=1}^I r_{ik}} {\sum_{j=1}^K \sum_{i=1}^I r_{ij}}$

$\mu_k^{t+1} = \frac{\sum_{i=1}^I r_{ik} x_i} {\sum_{i=1}^I r_{ik}}$

$\Sigma_k^{t+1} = \frac{\sum_{i=1}^I r_{ik}(x_i - \mu_k^{t+1}){(x_i - \mu_k^{t+1})}^T} {\sum_{i=1}^I r_{ik}}$

このように、それぞれの正規分布の分布全体への影響度を示す $\lambda_k$ は負担率の合計値が、平均や分散は負担率を重みとして利用した値が使われており、直感的にも妥当そうな更新式になっています。

以上で混合正規分布を用いたEMアルゴリズムの推定の方法となります。EMアルゴリズムは他にも様々な分布に適用することができます。外れ値に強くするために用いられるt-分布、部分空間を用いて正規分布の表現力を上げる因子分析、そしてt-分布に部分空間を適用したものなど、それこそ適用方法は様々で、詳しくは上記テキストの7.5 - 7.7章をご確認ください。

判別モデルによる回帰

次に、判別モデルによる回帰の方法を説明していきます。判別モデルでの回帰は、上で示したように ${ Pr(w|x,\theta) = Norm_w(\phi_0 + \phi_1 x, {\sigma}^2) }$ と書くことができます。この式に対して下記のような最尤推定を行うことを考えます。

$\hat{\theta} = argmax(Pr(w|x,\theta))$

この式を解くと、推定値は解析的に求めることができ下記のようになります。

$\hat{\phi} = (XX^T )^{-1}Xw$

$\hat{\sigma}^2 = \frac{(w - X^T \phi)^T (w - X^T \phi)} {I}$

以上が判別モデルでの線形回帰の例となります。しかし、この単純なモデルで単純な推定を行うと問題があります。一つは、学習データ点のない部分についても過剰に推定の信頼性がたかくなってしまうover-confident、さらに線形関数にしか適用できない点、さらに過学習や高次元な観測値の場合パラメーターの次元が大きくなりすぎる点が挙げられます。

over-confident に対しては、ベイズ推定を用いることで、非線形化は関数の非線形化およびカーネルトリック、過学習やパラメーターの次元の削減にはスパース回帰を用いることで改善できます。これらの手法は組み合わせることにより、Gaussian Process Regression やRelative Vector Regression などより高度な手法でモデリングすることができます。

それぞれを簡単に見ていきましょう。まず、ベイズ推定による線形回帰です。これは何も難しいことはありません。ベイズ推定により観測データが与えられた時のパラメーターの確率分布を推定すればよいのです。

$Pr(\phi|X,w) = \frac{Pr(w|X,\phi)Pr(\phi)} {Pr(w|X)}$

ですので、 ${ Pr(w|x,\theta) = Norm_w(\phi_0 + \phi_1 x, {\sigma}^2) }$ で、 $Pr(\phi) = Norm_\phi(0, \sigma^2 I)$ とおくと、 $Pr(\phi|X,w)$ を解析的に求めることができます。また新しい観測値 $x^*$ が与えられた際に、

$Pr(w^*|x^*,X,w) = \int Pr(w^*|x^*,\phi) Pr(\phi|X,w) d \phi$

も解析的に求めることができます。このように推定を行うと、学習データに近い部分のxとwの組み合わせの信頼性は高く、そうでない点については信頼性があまり高くない分布を得ることができます。ちなみに上記の推定はパラメーター $\phi$ についてのものだけであり、 $\sigma^2$ については別途推定をする必要があります。このテキストでは、 $\phi$ についての周辺尤度 $Pr(w|X,\sigma^2)$ を最大化する $\sigma^2$ を用いて、説明を行っています。

次に非線形関数の例です。非線形関数では $z_i=f(x_i)$ となる非線形関数fを定義し、 $Pr(w|X,\theta) = Norm_w(Z^T\phi,\sigma^2I)$ とおいてパラメーターを推定します。関数fにはRBFやarctanなどが用いられます。いずれにおいても、最尤推定もベイズ推定も線形回帰と同様に行うことができます。

ただし、新たな観測値xに対してwを計算すると逆行列等の行列計算が非常に大変になります。そこでカーネルトリックを用いてzの内積計算を簡単になるような関数fを用います。RBFカーネルなどが有名なカーネル化関数となります。

スパースな回帰では、 $Pr(\phi)$ を正規分布をステューデントのt-分布を用いることにより、パラメーターが各軸が独立になるように学習させます。具体的には下記のように定義します。

$Pr(\phi) = \prod_{d=1}^D Stud_{\phi_d}(0,1,\nu)$

しかしこの分布を用いると単純に解析的に推定を行うことができません。そこで、 $\phi$ について事後確率分布を周辺化することによりその値を求めることにより、直接事後確率分布を計算します。

また、次元を削減するには双対モデルを解く方法もあります。双対モデルは、パラメーター $\phi$ を $\phi = X \psi$ とおくことでパラメーターを推定する方法です。このモデルの面白いところは、元の式で推定した場合と同じ結果を得ることができるところです。ここで $\phi$ の次元は変数の次元数に依存し、 $\psi$ の次元は学習データの数に依存します。つまり、双対モデルを用いることにより、求めるパラメーターの次元数を削減することができます。