はじめに

前回の記事では、特定の分布や関係を仮定せず、与えられたデータに柔軟に対応することのできるモデルということで、ガウス過程回帰について取り上げました。

本記事ではガウス過程をさらに深堀りし、Stanなどでガウス過程を自在に扱うための土台を固めたいと思います。

本記事の構成は以下の通りとします。

• はじめに
• ガウス過程の導入
• ガウス過程潜在変数モデル
  • 多変量正規分布のサンプリング
  • ガウス過程潜在変数モデルとは
  • Stanでの実装
• ガウス過程の予測分布
• カテゴリカルな変数を用いたカーネル
  • 質的変数が1つの場合
  • 質的変数が2つ以上の場合
  • 制約のあるPDUDE
• まとめ

特に以下2章は力を入れました。

• 「ガウス過程潜在変数モデル」 Stan等でガウス過程を様々なモデルに活用するために不可欠
• 「カテゴリカルな変数を用いたカーネル」 質的変数が入力に含まれる場合にも、ガウス過程を適用できるようにするための方法を説明

当初本記事と次回の記事はひとつの記事として公開していましたが、余りにも長すぎる記事だったので、二つに分けました。

ガウス過程の導入

まず、ガウス過程を用いたモデルについて一般化します。
入力の要素数を数を$I$とし、入力データの1サンプルを$\mathit{x}=(x_1,\dots ,x_I{)}^T$と表記し、全入力を$N\times I$行列$\mathit{x}$を用いて$\mathit{X}=({\mathit{x}}_1,\dots ,{\mathit{x}}_N{)}^T$とおきます。また出力データを$\mathit{y}=(y_1,\dots ,y_N{)}^T$とおきます。

ガウス過程モデルの一般式は以下になります。

$$\begin{array}{}\text{(1)}& y(\mathit{x})=\mathit{\beta }\mathit{f}(\mathit{x})+ϵ(\mathit{x})\end{array}$$

ここで$\mathit{f}=(f_1(\mathit{x}),\dots ,f_l(\mathit{x}))$、$\mathit{\beta }=({\beta }_1,\dots ,{\beta }_l)$はそれぞれ、ガウス過程モデルとは別で設定された特徴ベクトルとその重みです。これらはガウス過程を適用する以前からデータにあてはめられたモデルになります。

$ϵ(\mathit{x})$は右辺第一項のモデルで説明しきれなかった残差としての扱いですが、この残差がガウス過程に従うと仮定します。ガウス過程モデルは得られたデータに柔軟に近づこうとする挙動をとりますので、この性質を利用して右辺第一項で説明できなかった部分をガウス過程で補おうという目論見です。

各入力の誤差項が同じ標準偏差$\sigma$をとると設定できたとき、$\mathit{ϵ}(\mathit{x})=(ϵ({\mathit{x}}_1),\dots ,ϵ({\mathit{x}}_N))$が以下のガウス過程に従うとします。

$$\begin{array}{}\text{(2)}& \mathit{ϵ}(\mathit{x})\sim \mathrm{Normal}(0,{\sigma }^2\mathit{K})\end{array}$$

上式では$ϵ$が残差を説明するものであることから多変量正規分布の平均を0としています。また分散共分散行列$\mathit{K}$は、ハイパーパラメータ$\varphi$を持つカーネル関数$k_{\varphi }({\mathit{x}}_m,{\mathit{x}}_n)$を要素に持つ$N$×$N$行列です($m,n=1,\dots ,N$)。

$$\begin{array}{}\text{(3)}& \mathit{K}(m,n)=k_{\varphi }(ϵ({\mathit{x}}_m),ϵ({\mathit{x}}_n))\end{array}$$

各入力の標準偏差を$\sigma$で設定している為、$\mathit{K}$の$(m,n)$成分は、$ϵ({\mathit{x}}_m)$と$ϵ({\mathit{x}}_n)$の相関に等しくなり、

$$\begin{array}{}\text{(4)}& k_{\varphi }(ϵ({\mathit{x}}_m),ϵ({\mathit{x}}_n))=\mathrm{cor}(ϵ({\mathit{x}}_m),ϵ({\mathit{x}}_n))\end{array}$$

となります。

ガウス過程では、カーネル関数によって${\mathit{x}}_m$と${\mathit{x}}_n$の距離を定める空間を決定し、入力${\mathit{x}}_m$と${\mathit{x}}_n$の距離によって出力$ϵ({\mathit{x}}_m)$と$ϵ({\mathit{x}}_n)$の近さを決定します。また、カーネル関数によって決まる値は、$\mathrm{cor}(ϵ({\mathit{x}}_m),ϵ({\mathit{x}}_n))$と解釈することができる、ということです。

ガウス過程潜在変数モデル

多変量正規分布のサンプリング

$(2)$式で用いられる多変量正規分布について、ランダムにサンプルを得る方法を紹介します。

平均0、分散共分散行列$\mathbf{\Sigma }$の多変量正規分布からのサンプルを得る場合、まず

$$\begin{array}{}\text{(5)}& \mathbf{\Sigma }=\mathit{L}{\mathit{L}}^T\end{array}$$

を満たす行列$\mathit{L}$を求めます。$(4)$式のような行列の分解はコレスキー分解と呼ばれます。

次に、標準正規分布からの乱数$\mathit{x}=(x_1,\dots ,x_N)$を生成します。

$$\begin{array}{}\text{(6)}& \mathit{x}\sim \mathrm{Normal}(0,1)\end{array}$$

$\mathit{y}=\mathit{L}\mathit{x}$の分布は、

に$\mathit{x}={\mathit{L}}^{-1}\mathit{y}$を代入すると、変数変換による空間の単位当たり面積の変動を調整するヤコビアン$|\partial \mathit{y}/\partial \mathit{x}|$は定数だから、

となります。

このことから、$\mathrm{Normal}(0,\mathbf{\Sigma })$に従う乱数を生成するには、標準正規分布に従う$\mathit{x}$をランダムに生成し、$\mathit{y}=\mathit{L}\mathit{x}$と変換すればよいと分かります。

ガウス過程潜在変数モデルとは

モデルの残差$ϵ(\mathit{x})$がガウス過程に従うとした$(2)$式を、先ほど紹介した多変量正規分布の乱数生成法を用いて変形すると、

$$\mathit{K}=\mathit{L}{\mathit{L}}^T$$ $$\begin{array}{}\text{(9)}& \mathit{\eta }\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(0,1)\end{array}$$ $$\mathit{ϵ}(\mathit{x})=\mathit{L}\mathit{\eta }$$

となります。（$\mathit{\eta }=({\eta }_1,\dots ,{\eta }_n)$） このように、潜在変数$\mathit{\eta }$を用いたガウス過程は、ガウス過程潜在変数モデル(Latent variable GP)と呼ばれ、出力が正規分布でないとき等に有用です。

Stanでの実装

以下、Stanマニュアルを引用してLatent variable GPの実装について軽く触れておきます。

Latent variable GPのStanでの実装は以下のようになります。

data{
  int<lower=1> N;
  real x[N];
  vector[N] y;
}

transformed data{
  real delta = 1e-9;
}

parameters {
  real<lower=0> rho;
  real<lower=0> alpha;
  real<lower=0> sigma;
  vector[N] eta;
}

model {
  vector[N] f;
  {
    matrix[N, N] L_K;
    matrix[N, N] K = cov_exp_quad(x, alpha, rho);

    // diagonal elements
    for (n in 1:N)
      K[n, n] = K[n, n] + delta;

    L_K = cholesky_decompose(K);
    f = L_K * eta;
  }
  eta ~ std_normal();
  y ~ normal(f, sigma);
}

ここで、K = cov_exp_quad(x, alpha, rho)はガウスカーネルをつくる便利な関数で、

を要素に持つカーネル行列$\mathit{K}$を作成してくれます。

また、for (n in 1:N) K[n, n] = K[n, n] + delta;とすることで、カーネル行列の対角要素に微小な値を加えていますが、こうすることでカーネル行列の逆行列の計算を安定化させています。また、コレスキー分解はその対象が正定値行列であることが必須ですが、対角要素に微小量を加えることで、その行列が正定値行列であることを保証することができます。

式$(9)$の3式の計算は、それぞれL_K = cholesky_decompose(K);、eta ~ std_normal();、f = L_K * eta;と指定しています。

上のコードでは、残差ではなく出力$\mathit{y}$の平均値が$\mathit{f}\sim \mathrm{Normal}(0,{\mathit{K}}_{\varphi })$に従うと設定し、y ~ Normal(f, sigma)とすることで、平均$\mathit{f}$の正規分布に従うと設定していますが、この部分をポアソン分布等他の分布にすることで、様々なモデルを構築することが出来ます。

例えば、0か1のみをとる出力$\mathit{y}$にベルヌーイ分布を仮定し、$y_n=1$となる確率$\mathit{p}=(p_1,\dots ,p_N)$を、入力${\mathit{x}}_n$のガウス過程を用いて説明する場合、

$$\mathit{y}\sim \mathrm{Binomial}(\mathit{p})$$ $$\begin{array}{}\text{(11)}& \mathit{p}=\mathrm{inverselogit}(\mathit{f}) \end{array}$$ $$\mathit{f}\sim \mathrm{Normal}(\mu ,\mathit{K})$$

と表現でき、これをStanで実行する場合、以下のようになります。

data{
  int<lower=1> N;
  real x[N];
  vector[N] y;
  ...
}

parameters{
  real mu; //muはpの期待値 観測データが近くに無い場合に漸近する値
  ...
}

transformed parameters{
  vector[N]<lower=0, upper=1> p;
  ...
  p = mu + f;
  ...
}
...
model {
  mu ~ std_normal()
  ...
  y ~ bernoulli(p);
}

ガウス過程の予測分布

ガウス過程モデルにおいて、入力$\mathit{x}$に含まれない値${\mathit{X}}^{\ast }=({\mathit{x}}_1^{\ast },\dots ,{\mathit{x}}_N^{\ast })$に対応する出力の値${\mathit{y}}^{\ast }=(y_1^{\ast },\dots ,y_M^{\ast })$を予測したい場合、$\mathit{y}$と${\mathit{y}}^{\ast }$の同時分布を次のようにすればよいです。

$$\begin{array}{}\text{(12)}& \left(\begin{array}{c}{y}_1\\ ⋮\\ y_N\\ y_1^{\ast }\\ ⋮\\ y_M^{\ast }\end{array}\right)\sim \mathrm{N}\mathrm{o}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}(\overrightarrow{0},\left(\begin{array}{cc}\mathit{K}& {\mathit{k}}^{\ast }\\ {\mathit{k}}^{\ast T}& {\mathit{k}}^{\ast \ast }\end{array}\right))\end{array}$$

ここで、${\mathit{k}}^{\ast }(n,m)=k_{\varphi }({\mathit{x}}_n,{\mathit{x}}_m^{\ast })$、${\mathit{k}}^{\ast \ast }(m,m)=k_{\varphi }({\mathit{x}}_m^{\ast },{\mathit{x}}_m^{\ast })$です。

カテゴリカルな変数を用いたカーネル

前回記事も含めこれまでは連続型変数を扱うことを前提にしていましたが、連続的な値をとらず、絶対的な大小関係ももたない質的変数（カテゴリカルな変数）が入力に含まれる場合、量的変数と質的変数の両方の性質を考慮した空間を定義することのできるカーネルを設定する必要があります。しかし、質的変数には「距離」の概念が無いため、ガウスカーネル等のように、各入力の「近接性」を再現するカーネルで対応することはできません。では、どのようにカーネル関数を設定すればよいのでしょうか。

質的変数が1つの場合

まず、量的変数と質的変数を含んだ入力を$\mathit{w}=({\mathit{x}}^T,{\mathit{z}}^T{)}^T$とし、$\mathit{x}=(x_1,\dots ,x_I)$を量的変数、$\mathit{z}=(z_1,\dots ,z_J)$を質的変数とします。それに従い、$(1)$式、$(2)$式を

$$\begin{array}{}\text{(1’)}& y(\mathit{w})=\mathit{\beta }\mathit{f}(\mathit{w})+ϵ(\mathit{w})\end{array}$$

$$\begin{array}{}\text{(2’)}& ϵ(\mathit{z})\sim \mathrm{Normal}(0,{\sigma }^2\mathit{K})\end{array}$$

としておきます。

簡略のため、$m_1$個の値をとる一つの質的変数$z_1$について考えます。 $\mathit{w}=({\mathit{x}}^T,z_1)=({\mathit{x}}^T,u)$($u=1,\dots ,m_1$)における$\mathit{ϵ}(\mathit{x})$を、

$$\begin{array}{}\text{(13)}& {\mathit{ϵ}}^{\ast }(\mathit{x})=\left(\begin{array}{c}{ϵ}_1(\mathit{x})\\ ⋮\\ {ϵ}_{m_1}(\mathit{x})\end{array}\right)\end{array}$$

と定義します。すると、量的変数については${\mathit{ϵ}}^{\ast }(\mathit{x})$の各要素内でのガウス過程で完結させ、質的変数の影響については、${\mathit{ϵ}}^{\ast }(\mathit{x})$の各要素間の分散共分散行列を決定すれば、量的変数・質的変数双方の差異を考慮した相関関数$\mathrm{cor}(ϵ({\mathit{w}}_m),ϵ({\mathit{w}}_n))$を指定できそうです。そこで、${\mathit{ϵ}}^{\ast }(\mathit{x})$を

$$\begin{array}{}\text{(14)}& {\mathit{ϵ}}^{\ast }(\mathit{x})=\mathit{A}\mathit{\eta }(\mathit{x})\end{array}$$

と推定することにします。

ここで、$\mathit{\eta }(\mathit{x})=({\eta }_1(\mathit{x}),\dots ,{\eta }_{m1}(\mathit{x}){)}^T$は量的変数における「距離」を考慮する部分で、各要素がそれぞれ独立に、標準偏差$\sigma$、相関関数$\mathit{K}$のガウス過程に従って生成されるものとします。

また、$m_1\times m_1$行列$\mathit{A}$は質的変数の影響を考慮する部分で、単位行ベクトル${\mathit{a}}_u$(${\mathit{a}}_u{\mathit{a}}_u^T=1$、$u=1,\dots ,m_1$)を用いて

$$\begin{array}{}\text{(15)}& \mathit{A}=\left(\begin{array}{c}{\mathit{a}}_1\\ ⋮\\ {\mathit{a}}_{m1}\end{array}\right)\end{array}$$

とし、${ϵ}_i(\mathit{x})$($i=1,\dots ,m_1$)を単位行ベクトル${\mathit{a}}_i$の指定する重みに基づく$\mathit{\eta }(\mathit{x})$の要素の線形和で表現することにします。

すると、

(${\eta }_i({\mathit{x}}_n)$、${\eta }_{i^{\prime }}({\mathit{x}}_m)$は$i=i^{\prime }$のときのみ$k_{\varphi }({\mathit{x}}_m,{\mathit{x}}_n)$の相関をとる)だから、2つの入力${\mathit{w}}_i=({\mathit{x}}_i^T,z_i{)}^T$($i=m,n$)の相関関数$\mathrm{c}\mathrm{o}\mathrm{r}(ϵ({\mathbf{w}}_{\mathrm{m}}),ϵ({\mathbf{w}}_{\mathrm{n}}))$について、

が成り立ちます($zm,zn=1,\dots ,m_1$)。

ここで、${\tau }_{r,s}={\mathit{a}}_r^T{\mathit{a}}_s$($r=zm,s=zn,r,s=1,\dots ,m_1,$)とおくと、半正定値行列$\mathit{T}=\mathit{A}{\mathit{A}}^T$は、

と、${\tau }_{r,s}$を$(r,s)$成分にもつ対角成分が1の行列(positive definite matrix with unit diagonal elements)になります。

以降、この性質をもつ行列をPDUDEと書きます。

以上のことから、

$$\begin{array}{}\text{(19)}& \mathit{T}(r,s)k_{\varphi }({\mathit{x}}_m,{\mathit{x}}_n)=\mathrm{c}\mathrm{o}\mathrm{r}(ϵ({\mathbf{w}}_{\mathrm{m}}),ϵ({\mathbf{w}}_{\mathrm{n}}))\end{array}$$

は、質的変数と量的変数の影響を考慮することのできる相関関数ととらえることができます。

質的変数が2つ以上の場合

一般的なケースとして、$J$個の質的変数$\mathit{z}=(z_1,\dots ,z_J{)}^T$の場合を考えます。

ここで、$z_j(j=1,\dots ,J)$は$1,\dots ,m_j$の値をとるものとします。すると、式$(19)$の拡張により、$ϵ(\mathit{w})$の相関は関数は以下のようになります。

$$\begin{array}{}\text{(20)}& \prod _{j=1}^J({\tau }_{z_{j,r,s}}{k}_{\varphi }({\mathit{x}}_m,{\mathit{x}}_n))=\mathrm{c}\mathrm{o}\mathrm{r}(ϵ({\mathbf{w}}_{\mathrm{m}}),ϵ{\mathbf{w}}_{\mathrm{n}})\end{array}$$

${\tau }_{j,r,s}$は、$J$番目のPDUDE${\mathit{T}}_j$の$(r,s)$成分です。

特に、$k_{\varphi }({\mathit{x}}_m,{\mathit{x}}_n)$に式$(10)$のガウスカーネルをもちいた場合、$(20)$式は以下のようになります。

式$(21)$は、量的変数の影響は$exp(-\frac{1}{2{\rho }^2}\sum _{i=1}^I(x_{im}-x_{in}{)}^2)$で考慮し、質的変数の影響はそれとは独立に$\left(\prod _{j=1}^J{\tau }_{z_{j,r,s}}\right)$で考慮する、という構造を持っています。パラメータ${\tau }_{z_{j,r,s}}$は、質的変数$z_j$について$r=z_{j}m$をとる入力${\mathit{w}}_m$と、$s=z_{j}n$をとる入力${\mathit{w}}_n$の、$z_j$のみによる共通性(相関)への影響を反映する役割を担っています。なお、式$(21)$ではshapeパラメータ${\alpha }^2$を考慮していませんが、これは$(2)$式でshapeパラメータの役割を$\sigma$が受け持っている為です。

モデリングにおいては、${\tau }_{z_{j,r,s}}$は正の値をとるように設定し、任意の2点の入力が無相関もしくは正の相関のみをとるようにします。

制約のあるPDUDE

前節では、PDUDEについて制約を設けない相関行列を用いていました。柔軟なモデリングにおいてはこれで問題ないのですが、質的変数が順序尺度であったり、カテゴリカルな変数であったりするということがあらかじめ自明な場合は、PDUDEに制約を持たせることで、その情報をモデルに反映させることができます。ここでは、sesがカテゴリカルな変数であることから、質的変数がカテゴリカルな場合にPDUDEに設ける制約について説明します。

結論からですが、$m$個の値をとる$z$がカテゴリカルな場合、下記の${\tau }_{r,s}$を$(r,s)$成分に持つ等方性を持った$m\times m$相関行列$\mathit{T}$が用いられます。

$$\begin{array}{}\text{(22)}& {\tau }_{r,s}=\{\begin{array}{l}c(0<c<1)(r\ne s)\\ 1(r=s)\end{array} \end{array}$$

$r=s$のとき、入力間の相関は$z$に関しては1、$r\ne s$のとき、c(一定)にする、ということです。 $\mathit{T}$は以下のように分解できます。

このとき、任意の$m\times 1$ベクトル$\mathit{a}$について、 $$\begin{array}{}\text{(24)}& {\mathit{a}}^T\mathit{T}\mathit{a}=(1-c){\mathit{a}}^T\mathit{a}+c({\mathit{a}}^{T}1{)}^2>0\end{array}$$ だから、$\mathit{T}$は正定値行列なので、PDUDEです。$\mathit{T}$が正定値行列であることは結構重要で、各要素の値を出力するカーネル関数$k({\mathit{z}}_m,{\mathit{z}}_n)$が何かしらの特徴ベクトル空間の内積を表現するために必要な条件です。このあたりの詳細はこちらなどを見てください。また上記$\mathit{T}$の各要素を出力する関数はisotropic correlation functionと呼ばれ、isotropic correlation functionによる出力を要素に持つ行列はcompound symmetric correlation matrixと呼ばれています。

上記の$\mathit{T}$をPDUDEに用いる場合、式$(21)$は下記のように変形できます。

ここで、$I[r\ne s]$は下記の関数になります。

$$\begin{array}{}\text{(26)}& I[r\ne s]=\{\begin{array}{l}1(r\ne s)\\ 0(r=s)\end{array} \end{array}$$

式$(25)$最後の項は、対数をとると

となります。よって、対数スケールにおいて、量的変数についてはL2距離を、質的変数については0~1の値をとる距離を使用していることが分かります。

以上、カテゴリカルな変数を用いたカーネルについて説明しました。参考文献はこちらになります。制約のあるPDUDEについては、今回紹介したもののほかにも順序尺度に対するPDUDE,グループ相関に対するPDUDE等紹介されています。

まとめ

今回はガウス過程みまつわる理論を深堀りしてみました。ガウス過程潜在変数モデル、カテゴリカルな変数への対応の2点が重要です。ここで紹介した内容は、次回の記事で実際のデータへ応用し、期待通り機能することを確かめたいと思います。