orizuru

つながる.見える.わかる IoTソリュ-ション

Automatic Relevance Determination(ARD)

約 14 分
Automatic Relevance Determination(ARD)

はじめに

 今回は、Automatic Relevance Determination(ARD)について解説する。ARDとは、目的変数に対する個々の説明変数の寄与の大きさを見積もる手法である。今回は以下の2つの文脈でARDを行う。

  • ベイズ推定を用いた回帰によりARDを行う。
  • ガウス過程を用いた回帰によりARDを行う。
  • ベイズ推定を用いた回帰によるARD

    定式化

     目的変数をX=(x_0,\cdots,x_{N-1})、説明変数をY=(y_0,\cdots,y_{N-1})とする。ここで、x_n\in\mathbb{R}^M,y_n\in\mathbb{R}である。潜在変数\thetaを導入し、同時確率分布p(X,Y,\theta)を考える。この分布にベイズの定理を適用すると次式を得る。

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

    尤度p(Y|X,\theta)として

    (2)    \begin{eqnarray*} p(Y|X,\theta)&=&\prod_{n=0}^{N-1}p(y_n|x_n,\theta) \\ p(y_n|x_n,\theta)&=&\mathcal{N}(y_n|w^T x_n,\alpha^{-1}) \end{eqnarray*}

    を仮定する。\mathcal{N}(y|\mu,\sigma^2)は平均\mu、標準偏差\sigmaを持つ正規分布を表す。\theta=(w,\alpha)とした。w\in \mathbb{R}^M,\alpha\in\mathbb{R}である。これらのパラメータに対する事前確率p(\theta)として次式を仮定する。

    (3)    \begin{eqnarray*} p(\theta) &=& p(w,\alpha)\\ &=& p(w|\lambda)p(\lambda)p(\alpha) \end{eqnarray*}

    ここで、新たなパラメータ\lambda=(\lambda_0,\cdots,\lambda_{M-1})を導入した。各確率分布を次式で定義する。

    (4)    \begin{eqnarray*} p(w|\lambda)&=&\mathcal{N}(w|0_M,A^{-1})\\ &=&\prod_{m=0}^{M-1}\mathcal{N}(w_m|0,\lambda_m^{-1}) \\ A&=&diag(\lambda_0,\cdots,\lambda_{M-1}) \\ p(\alpha)&=&{\rm Gam}(\alpha|a,b) \\ p(\lambda_m)&=&{\rm Gam}(\lambda_m|c,d) \\ \end{eqnarray*}

    {\rm Gam}(\lambda|c,d)はハイパーパラメータc,dを持つガンマ分布、0_MM次元のゼロベクトル、AM\times Mの対角行列である。w,\alpha,\lambdaは推定すべきパラメータであり、a,b,c,dはエビデンスp(Y|X)が最大になるように決めるハイパーパラメータである。

     さて、式(2)は回帰式を

    (5)    \begin{eqnarray*} y_n&=&w^T x_n \\ &=&w_0 x_{n0} + \cdots + w_{M-1} x_{nM-1} \end{eqnarray*}

    と仮定していることを意味する(ノイズ部分は無視した)。つまり、w_mは説明変数x_{nm}の寄与の大きさを表す。これがベイズ推定によるARDの仕組みである。

     仮定した全ての確率分布を式(1)に代入すると、次の事後確率分布を厳密に(解析的に)求めることができる。

  • p(w|X,Y)
  • p(\alpha|X,Y)
  • p(\lambda|X,Y)

  • 事後確率分布p(w|X,Y)から、個々の説明変数の寄与の大きさを知ることができる。また、次式により未知データx_{*}に対するy_{*}を予測することができる。

    (6)    \begin{eqnarray*} p(y_*|x_*,X,Y)=\int dw\;d\alpha\;d\lambda\;p(y_*|x_*,w,\alpha,\lambda)p(w|X,Y)p(\alpha|X,Y)p(\lambda|X,Y) \end{eqnarray*}

    データ作成

     回帰に用いるデータを次の手順で作成した。

  • N=130,M=3とした。つまり、3次元ベクトルx_n=(x_{n0},x_{n1},x_{n2})を130個用意する。
  • [0,1)の範囲で乱数を生成し、これをx_{n0}とした。
  • [0,100)の範囲で乱数を生成し、これをx_{n1}とした。
  • [0,0.1)の範囲で乱数を生成し、これをx_{n2}とした。
  • 次式によりy_nを生成した。

    (7)    \begin{equation*} y_n=\cos\left({\frac{x_{n0}\;x_{n2}}{2}\right)+\frac{x_{n0}\;x_{n2}}{10} \end{equation*}

     
    説明変数x_{n1}は目的変数y_nに全く寄与せず、各説明変数のスケールは大きく異なる状況をシミュレートしたものである。訓練データを100個、テストデータを30個とした。前者で学習を行い、後者で検証を行う。

    学習

     学習の前に以下の前処理を行う。

    (8)    \begin{eqnarray*} x_{nm}^{\prime}&=&\frac{x_{nm}-\mu_{m}}{\sigma_{m}}\\ \mu_{m}&=&\frac{1}{N}\sum_{n=0}^{N-1} x_{nm}\\ \sigma_{m}^2&=&\frac{1}{N}\sum_{n=0}^{N-1} (x_{nm}-\mu_m)^2 \end{eqnarray*}

    各軸(各次元)ごとに、平均値を0に標準偏差を1に正規化する。このあと、ARDを行う。Pythonライブラリsklearnに、上の「定式化」で示したロジックを実装したクラスARDRegressionがある。これを用いると、学習部分のコードは以下のようになる(ソースはここのsample_sklearn.pyである)。

    3行目でARDRegressionオブジェクトを作成し、4行目で学習する。次に結果を示す。

    結果

     最初に示すのは、(w_0,w_1,w_2)の平均値と標準偏差である。これはp(w|X,Y)から求まる量である。

    期待通りの結果である。x_0,x_2からの寄与が大きく、x_1からの寄与はほとんどない。標準偏差は小さ過ぎるので上図では判別できなくなっている。次に示すのは30個のテストデータを予測した結果である。

    |{\rm Truth - Prediction}|の平均値は0.266となった。

    ガウス過程を用いた回帰によるARD

    定式化

     出発点は先と同じである。

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

     
    パラメータ\thetaをガウス過程に使うカーネルに関わる部分\theta_Kとそうでない部分\theta_Lの2つに分離する。さらに、XYの間に次式を仮定する。

    (10)    \begin{eqnarray*} y_n&=&f(x_n)+\epsilon_n \\ \epsilon_n&\sim&\mathcal{N}(\epsilon_n|0,\sigma^2) \end{eqnarray*}

    ここで、\mathcal{N}(\epsilon_n|0,\sigma^2)は平均0、分散\sigma^2の正規分布を表す。各観測値は独立同分布から生成されると仮定すると、尤度は

    (11)    \begin{eqnarray*} p(Y|X,\theta_L) &=&\prod_{n=0}^{N-1} p(y_n|x_n,\theta_L) \\ &=&\prod_{n=0}^{N-1} \mathcal{N}(y_n|f(x_n),\sigma^2)  \end{eqnarray*}

    と書くことができる。ここで、\theta_L=(f(\cdot),\sigma)とした。\sigmaの事前分布としてp(\sigma)を、関数f(\cdot)の事前分布として次のガウス過程を仮定する。

    (12)    \begin{equation*} p(f|\theta_K)=\mathcal{N}(f|0_N,K_{NN}^{\theta_K}) \end{equation*}

    ここで、f=\left(f(x_0),\cdots,f(x_{N-1})\right)とした。0_NN次元のゼロベクトル、K_{NN}^{\theta_K}はパラメータ\theta_Kを持つN\times Nの共分散行列、その成分(K_{NN}^{\theta_K})_{n,n^{\prime}}=k_{NN}^{\theta_K}(x_n,x_{n^{\prime}})がカーネル関数となる。カーネル関数として次のものを用いる。

    (13)    \begin{equation*} k_{NN}^{\theta_K}(x_n,x_{n^{\prime}})=\exp{\left(-\frac{1}{2}\sum_{m=0}^{M-1}\left(\frac{x_{nm}-x_{n^{\prime}m}}{l_m}\right)^2\right)} \end{equation*}

    ここで、\theta_K=(l_0,\cdots,l_{M-1})とした。l_mは各軸のスケール長を表す。この逆数が説明変数の重みに相当する。さらに、\theta_Kの事前分布p(\theta_K)も導入する。ここまでの結果をまとめると、式(9)は次のようになる。

    (14)    \begin{equation*} p(f,\sigma,\theta_K|X,Y)=\frac{p(Y|X,\sigma,f)\;p(\sigma)\;p(f|\theta_K)\;p(\theta_K)}{p(Y|X)} \end{equation*}

    fの事前分布p(f|\theta_K)がガウス過程であることに注意する。上式をベイズ推定により評価することもできるが、ここでは右辺分子を取り出しMAP推定を行う。右辺分子の対数を取ると

    (15)    \begin{equation*} L=\log{\left(p(Y|X,\sigma,f)\;p(\sigma)\;p(f|\theta_K)\;p(\theta_K)\right)} \end{equation*}

    を得る。これを最大にするような\theta_K\sigmaを求める。

    (16)    \begin{equation*} (\theta_K^{MAP},\sigma^{MAP})={\rm arg}\max_{\theta_K,\sigma}L \end{equation*}

    この2つのパラメータについては点推定であることに注意する(確率を求めるのではない)。未知の値x_{\star}に対するf_\star=f(x_\star)の確率分布は次式から計算される。

    (17)    \begin{equation*} p(f_\star|\theta_K^{MAP})=\int df\;p(f_\star|f,\theta_K^{MAP})\;p(f|\theta_K^{MAP}) \end{equation*}

    右辺の積分内の最初の項p(f_\star|f,\theta_K^{MAP})p(f(x_\star)|f(x_0),\cdots,f(x_{N-1}),\theta_K^{MAP})を表し、式(12)から解析的に求めることができる条件付き確率分布である。上式からf_\starの分布が分かれば、未知の値x_{\star}に対するy_\starの値は次式から決定することができる。

    (18)    \begin{eqnarray*} f_\star&\sim&p(f_\star|\theta_K^{MAP}) \\ \epsilon_\star&\sim&\mathcal{N}(\epsilon_\star|0,\left(\sigma^{MAP}\right)^2) \\ y_\star&=&f_\star+\epsilon_\star \end{eqnarray*}

    以上で定式化できたので、実際に計算を行う。

    学習

     上の定式化で示したロジックは、ガウス過程のためのPythonライブラリGPyに実装されている。先と同じデータを使用し、先と同じ前処理を施したあと、以下を実行する(ソースはここのsample_gp_multi_dim.ipynbである)。

    2行目でカーネル関数(式(13)に相当する)を定義し、6行目と9行目でガウス過程による回帰を行う(ベイズ推定ではなくMAP推定が行われる)。

    結果

     最初に重み(1/l_0,1/l_1,1/l_2)の値を示す。

    期待通りの結果である。x_0,x_2からの寄与が大きく、x_1からの寄与はほとんどない。l_mは点推定される量であるから、「ベイズ推定を用いた回帰によるARD」の場合とは異なり標準偏差は存在しないことに注意する。次に、テストデータを予測した結果を示す。

    標準偏差は小さ過ぎて上図では判別できなくなっている。|{\rm Truth - Prediction}|の平均値は0.212となった。「ベイズ推定を用いた回帰によるARD」のときは0.266であったから改善されている。ガウス過程の場合、非線形回帰となるため精度が上がるのだろう。

    まとめ

     今回は、2つの文脈でAutomatic Relevance Determination(ARD)を見た。

    1. ベイズ推定を用いた回帰によるARD
    2. ガウス過程を用いた回帰によるARD

    1の回帰は線形回帰であるが、2の回帰は非線形回帰である。今回のデータに対しては、予測精度は2の方が高い。また、1の計算ではM\times Mの行列を、2の計算ではN\times Nの行列を扱う。したがって、今回のようなM<Nの場合は、1の方が高速に計算できる。

  • About The Author

    IoT/AIソリューション事業部(深層/機械学習・画像処理エンジニア)KumadaSeiya
    深層/機械学習と画像処理などを担当。物性理論で博士号を取得。
    http://seiya-kumada.blogspot.jp/
    https://twitter.com/seiya_kumada

    Leave A Reply

    *
    *
    * (公開されません)