5 線形回帰モデルと確率モデル

　モノづくりにおいて、品質に大きな影響を及ぼす重要なプロセスが設計プロセスであり、品質指標のバラツキである不確かさを定量化することが求められる。ここまで、設計対象の振舞いを代理モデルと表現するために線形回帰モデルを紹介してきた。観測データ集合は、試作モデルの測定結果でもよいが、より多くのデータを得るためにはCAEツールによる解析結果を利用することができる。
　線形回帰モデルは、入力データを基底関数によって写像した特徴ベクトルの線形結合で表される。この線形結合の重み係数を求めるために最小二乗法によって正規方程式を導き、重み係数と特徴ベクトルの内積によってスカラーである回帰値が得られた。
　今回のゴールは、線形回帰モデルを確率モデルとして捉えることである。

5.1 確率密度関数

　回帰モデルでは、出力は連続変数であるので、確率分布も確率変数の確率密度で表される。確率変数を\(\ X\ \)とするとき、確率変数がある値\(\ x\ \)となる確率密度を \[f\,(X=x)\ \tag{5.1}\] と表し、関数\(\ f\ \)を確率密度関数（以下、確率密度という）と呼ぶ。確率変数\(\ X\ \)を明示する必要がないときには \[p\,(x)=p\,(X=x)\ \tag{5.2}\] 　バラツキを確率モデルとして表現するには、標本サイズが小さい場合にも安定な近似が得られる正規分布（ガウス分布）を適用することが多い。
　確率変数が１次元の場合、平均\(\ \mu\ \)、分散\(\ \sigma_n^{\,2}\ \)の正規分布の確率密度は、 \[\ND(x/\mu/\sigma_n^{\,2}) = \frac{1}{\sqrt{2\pi \sigma_n^{\,2}}}\,\exp\, \left( -\frac{( {x-\mu)^2}}{2\sigma_n^{\,2}}\right) \tag{5.3}\] で与えられる。さらに、確率変数が\(\ N \ \)次元の場合には、平均\(\ \BS(\mu) \ \)、共分散行列\(\ \Sigma \ \)の正規分布に従っているとき、多変量正規分布の確率密度は、 \[\ND({\bf x}/\BS(\mu) /\Sigma ) = \frac{1}{(\sqrt{2\pi})^N\, \sqrt{|\Sigma |}}\,\exp\, \left( -\frac{1}{2}({\bf x}-\BS(\mu) )\TP \Sigma ^{-1}({\bf x}-\BS(\mu) )\right) \tag{5.4}\] で与えられ、\(\ |\Sigma |\ \)は共分散行列の行列式である。ここで、平均値ベクトルと\(\ N\times N\ \)共分散行列は、それぞれ次式の期待値で求められる。 \begin{cases} \ \BS(\mu) = \EV({\bf x}) \\ \lower 1em {\ \Sigma = \EV({\bf x}{\bf x}\TP ) – \EV({\bf x}) \EV({\bf x}) \TP } \tag{5.6} \end{cases} 　文脈から明らかな場合には、式 (5.3) と (5.4) をそれぞれ\(\ x \sim {\cal N}(\mu, \sigma_n^{\,2})\ \)、\(\ {\bf x} \sim {\cal N}(\BS(\mu) , \Sigma )\ \)と略記することがある。

5.2 線形回帰モデルを確率モデルとして表すための前提

　線形回帰モデルは、入力\(\ {\bf x} \in \RR^D\ \)を\(\ \BS(\phi) ({\bf x})\text{：　} \RR^D \mapsto \RR^N\ \)によって特徴量空間への写像し、特徴空間における線形結合\(\ {\bf w}\TP \BS(\phi) ({\bf x}) \ \)と出力\(\ y\in \RR \ \)から重みベクトル（回帰係数）\(\ {\bf w} \in \RR^N\ \)を推定した。そして、回帰値を\(\ \hat{y} = {\bf w}\TP \BS(\phi) ({\bf x})\ \)で求めた。この回帰モデルにおいて、\(\ y \ \)と\(\ {\bf w} \ \)を確率変数とするモデルとして表すために、以下の３つの前提条件を仮定する。

5.2.1 回帰値と観測値の差の確率モデル

　観測値\(\ y \ \)と回帰値\(\ \hat{y} \ \)の差\(\ \epsilon \ \)が平均\(\ 0\ \)、分散\(\ \sigma_n^{\,2}\ \)の正規分布に従うとする。 \begin{cases} \ y_i = {\bf w}\TP \BS(\phi) ({\bf x})+\epsilon \\ \lower 0.5em {\ p(\epsilon) = {\cal N}(\epsilon\,|\,0 , \sigma_n^{\,2}) } \tag{5.7} \end{cases} 　このとき、観測値 \(\ y \ \)が得られる条件付確率分布は、平均値\(\ {\bf w}\TP \BS(\phi) ({\bf x})\ \)にノイズ\(\ \epsilon\ \)が重畳したと見做せ、第１の前提条件 \[\CD(y/{\bf w},{\bf x}) = \ND(y/{\bf w}\TP \BS(\phi) ({\bf x}) /\sigma_n^{\,2}) \tag{5.8}\] が得られる。

5.2.2 重みベクトルの事前分布

　重みベクトルの各要素\(\ w_i\ \)の事前分布が、平均\(\ 0\ \)、分散\(\ \lambda^2\ \)の正規分布に従うとし、それぞれが独立であるとする。このとき、共分散行列は\(\ N\times N\ \)単位行列\(\ I_N\ \)の\(\ \lambda^2\ \)倍となり、第２の前提条件として\(\ {\bf w}\ \)の事前分布が導かれる。 \[\PD({\bf w}) =\prod_{i=1}^{N} \ND(w_i/0/\lambda^2) = \ND({\bf w}/\BS(0) /\lambda^2 I_N) \tag{5.9}\] 　一方で、\(\ {\bf w}\ \)の事後分布の確率密度\(\ \CD({\bf w}/y,{\bf x}) \ \) は、ベイズの定理によって \[\CD({\bf w}/y,{\bf x}) =\frac{\CD(y/{\bf x},{\bf w}) \,p({\bf w})}{\CD(y/{\bf x}) } \] と表される。右辺の\(\ \CD(y/{\bf x},{\bf w}) \ \)を尤度（likelihood）、\(\ \CD(y/{\bf x}) \ \)を周辺尤度(marginal likelihood)というが、周辺尤度は正規化定数である。従って、左辺の事後確率密度を最大化する（いわゆるMAP推定）ことと、右辺の分子の同時確率密度\(\ \CD(y/{\bf x},{\bf w}) \,p({\bf w})\ \)を最大化することは等価である。

5.2.3 出力と重みベクトルの同時確率分布

　線形回帰モデルを表す２つの確率変数である出力\(\ y\ \)と重みベクトル\(\ {\bf w}\ \)の確率分布を仮定したので、２つの確率分布の同時確率分布を次式で仮定し、第３の前提条件とする。 \[\CD(y,{\bf w}/{\bf x}) = \CD(y/{\bf w},{\bf x}) \, \PD({\bf w}) \tag{5.10}\]

5.3 確率モデルから得られる回帰モデル

　観測値集合\(\ {\cal D}=\{\,({\bf x}_i,y_i)\,|\,{\bf x}_i \in \RR^D \ ,\ y_i\in \RR\ ,\ i=1,\cdots ,n\,\} \ \)が与えられたときの同時確率密度は次式に式 (5.8) を代入すると次式が得られる。 \begin{eqnarray} {\cal P}({\bf w})&=&p({\bf w})\prod_{i=1}^n\,\CD(y_i/{\bf w},{\bf x}_i) \\ &=&p({\bf w})\prod_{i=1}^n\,\ND(y_i/{\bf w}\TP \BS(\phi) ({\bf x}_i) /\sigma_n^{\,2}) \end{eqnarray} \(\ n \ \)個のペア\(\ ({\bf x}_i,y_i)\ \)が互いに独立であれば、互いに独立な正規分布と多変量正規分布の関係 \[\prod_{i=1}^n\ND(y_i/{\bf w}\TP \BS(\phi) ({\bf x}_i)/\sigma_n^{\,2}) = \ND({\bf y}/\Phi{\bf w}/\sigma_n^{\,2}I_n) \] が成り立つので、 　 \[{\cal P}({\bf w})=p({\bf w})\,\ND({\bf y}/\Phi{\bf w}/\sigma_n^{\,2}I_n) \tag{5.11} \] が得られる。
　同時確率密度\(\ {\cal P}({\bf w})\ \)を最大化する\(\ {\bf w}\ \)を求めるために、対数関数が単調増加関数であることを利用すると \begin{eqnarray} \AG(max/{\bf w} \in \RR ^N / \log{\cal P}({\bf w})) &=& \AG(max/{\bf w} \in \RR ^N / \log{\ND({\bf y}/\Phi{\bf w}/\sigma_n^{\,2}I_n) }+ \log\ND({\bf w}/\BS(0) /\lambda^2I_N ) ) \\ &=& \AG(max/{\bf w} \in \RR ^N / -\frac{1}{2\sigma_n^{\,2}}({\bf y}-\Phi{\bf w})\TP ({\bf y}-\Phi{\bf w})-\frac{1}{2\lambda^2}{\bf w}\TP {\bf w} ) +const. \\ &=& \AG(min/{\bf w} \in \RR ^N / ({\bf y}-\Phi{\bf w})\TP ({\bf y}-\Phi{\bf w}) +\frac{\sigma_n^{\,2}}{\lambda^2}{\bf w}\TP {\bf w} ) \tag{5.12} \end{eqnarray} 　式 (5.12) は、観測値集合が与えられたとき、\(\ {\bf y}\ \)と\(\ {\bf y}\ \)の同時確率密度を最大化する\(\ {\bf w}\ \)を求めることと、誤差の平方和と重みベクトルのノルムの\(\ \alpha=\sigma_n^{\,2}/\lambda^2\ \)倍の和を最小化する\(\ {\bf w}\ \)を求めることが等価であることを示している。式 (5.12) は線形回帰を拡張したモデルでありリッジ回帰（ridge regression）と呼ばれる。

5.4 リッジ回帰モデル

　リッジ回帰の正規方程式を導出することとする。式 (5.12) から重回帰モデルの誤差の平方和を拡張すると \begin{eqnarray} \varDelta &=& ({\bf y}-\Phi{\bf x})\TP ({\bf y}-\Phi{\bf x})+\alpha\IP(x/x) \\ \frac{\partial \varDelta}{\partial {\bf w}}&=& -2\Phi\TP {\bf y}+2\Phi\TP \Phi{\bf w} +2\alpha I_N{\bf w} \tag{5.13} \end{eqnarray} が得られ、この結果からリッジ回帰の正規方程式が導かれる。 \[(\Phi\TP \Phi+\alpha I_N){\bf w}=\Phi\TP {\bf y} \tag{5.14}\] 　この正規方程式は、線形回帰モデルに現れた\(\ N\times N\ \)行列\(\ \Phi\TP \Phi\ \)に、\(\ \alpha I_N\ \)を加えた形式となっている。\(\ \Phi\TP \Phi\ \)の対角成分に\(\ \alpha\gt 0\ \)を加えたとも見做せる。半正定値行列\(\ \Phi\TP \Phi\ \)に正定値行列である正則化項を加えた行列は、正則な正定値行列となるので、\(\ (\Phi\TP \Phi+\alpha I_N)\ \)の逆行列は必ず存在し、重みベクトルを次式で求めることができる。極端な例として、観測データが１点の場合でも式 (5.15) で\(\ {\bf w}\ \)を求めることができる。
　行列\(\ \alpha I_N\ \)を正則化項（あるいは、ペナルティ項）と呼ぶ。 \[{\bf w}=(\Phi\TP \Phi+\alpha I_N)^{-1}\Phi\TP {\bf y} \tag{5.15}\]
　ここまで、単回帰→重回帰→線形回帰とモデルを拡張し、さらに重みと出力を確率変数として扱うモデルを構築した。重みについて事前分布を想定し、事後分布の確率密度を最大化することと等価な式 (5.10) で与えられる同時確率密度を最大化することによってリッジ回帰が導かれ、観測データ集合から重み\( \ {\bf w}\ \)が必ず求められることを示してきた。
　次回以降では、確率過程について解説する。

[参考] 線形回帰モデルのまとめ

　ここまで紹介した重要な関係を別ページにまとめて紹介する。
　式 (5.9) で示す重みベクトルの事前分布は同時同一分布ではなく、平均ベクトルが\(\ {\bf 0}\ \)、共分散行列が\(\ \Sigma_p\ \)の多変数正規分布に従うものとし一般化している。また導出は示していないが、テスト入力に対する出力の確率分布を求める式も紹介する。