ガウス過程シリーズ 1 概要

Stanのマニュアルの「Gaussian Processes」の章を実際に実行しましたので記録を残します。結論から言いますと、Stanでやる場合は回帰はよいですがクラス分類に使おうとすると計算が遅いし収束も悪いです。

まずGaussian Process（以下GPと呼ぶ）とは何ぞやということですがgpml（ぐぷむる？）として有名な次の書籍の1章が分かりやすいです。→Gaussian Processes for Machine Learning

これを咀嚼して勝手に補完してまとめたものが以下になります。

GPは教師あり学習の一手法です。教師あり学習では有限のトレーニングデータから関数を作ることになります。関数はありとあらゆる入力の値に対して予測値を返すものです。この関数を決めるにあたり、2つのアプローチがあります。1つめは関数をあるクラス（例えば線形だとか）に限定するものです。しかしこれは採用するクラスを誤るととたんに予測がうまくいかなくなる恐れがあります。強い非線形の関係があるにもかかわらず線形のクラスに限定した場合など。かといって使うクラスを豊富にするとオーバーフィッティングの恐れが出てきます。2つめはあらゆる可能性のある関数に事前分布を与えるものです。例えば滑らかな関数は他より確率が高いようにします。しかしながら無限のありうる関数を計算することはできません。ここでGPが処方箋となります。

GPは「正規分布に従う確率変数」の一般化です。processとはその確率変数を関数に拡張したものです。数学的な厳密さをさておくと、関数というものは、ベクトル内の1つ1つの要素がある入力のの値を持っているような非常に長いベクトルと考えることができます。GPは次のように定義されます。「GPは確率変数の集合であり、その集合の任意の有限個の変数は多変量正規分布に従う。」この条件によって上記の2つのアプローチの弱点をうまく避けています。多変量正規分布のパラメータは平均ベクトルと分散共分散行列の2つです。式で書くと、入力が個、があるとします。その一つ一つにスカラ値が個、が対応しています。その時、

のように長さの平均ベクトルと×のをパラメータとした多変量正規分布に従うと考える、これがGPです。解く問題の性質がパラメータに押し込まれており、これをどう決めるかがGPの大切なところです。また、GPはパラメータをいじることでBayesian線形モデル、スプラインモデル、（ある条件化で）大きなニューラルネットワーク、SVMなどのよく知られたモデルと数学的に等価とのことです。

の決め方（covariance function）は色々な選択肢がありますが最もよく使われるのは、の要素を入力との間の距離に応じて以下のように決めることです。

ここではクロネッカーのデルタと呼ばれるものでi=jの時だけ1,その他の場合は0の値を持ちます。はとが“近い”時のとの共分散（≒似てる具合）に相当し、はとがどれぐらいの距離だと“近く”なるかを決める特徴的な距離（lengthscale）、はのノイズの大きさに相当します。これらのhyper parameterは問題に応じて決まり、事前に与えることもできますが推定で求めることもできます。また、入力がD次元のベクトルになったとしても、例えばベクトルの差のdot_self（自分自身との内積）を使えば簡単に拡張できます。

ここからは実装です。さて、いまデータが以下のようにあります。

X	Y
-6.836951	-1.826294
-5.927641	0.0548737
-5.907641	0.0248737
-5.476121	0.3364621
-4.387844	-0.8187725
-3.739508	-1.2519856
-3.40779	-1.1943283
-2.537614	0.528471
-1.89467	1.592675
-0.683068	1.896029
0.6881621	-0.0568954
0.8978934	-0.7429603
1.129219	-1.0536914
2.512134	-2.835724
2.581766	-2.327798
4.268952	-1.3286307
4.427094	-1.0661972
4.944368	-1.634321
5.819585	-1.0883272
6.15919	-0.8673887

データは入力X1（1次元）がN1個、出力Y1（1次元）がN1個です。そして未知の入力値X2を入れた時にどんなy2になるか知りたいとします。つまり回帰曲線が求めたいわけです。この場合はX2も含めた大きな分散共分散行列を作って、それをパラメータとする多変量正規分布からY1とy2が生成されたと考えます。Stanコードは以下になります（マニュアルのほぼ写経です）。

data {
   int<lower=1> N1;
   int<lower=1> N2;
   vector[N1] X1;
   vector[N1] Y1;
   vector[N2] X2;
}
transformed data {
   int<lower=1> N;
   vector[N1+N2] X;
   vector[N1+N2] Mu;
   cov_matrix[N1+N2] Cov;

   N <- N1 + N2;
   for (n in 1:N1) X[n] <- X1[n];
   for (n in 1:N2) X[N1 + n] <- X2[n];
   for (i in 1:N) Mu[i] <- 0;
   for (i in 1:N)
      for (j in 1:N)
         Cov[i,j] <- exp(-pow(X[i] - X[j],2)) + if_else(i==j, 0.1, 0.0);
}
parameters {
   vector[N2] y2;
}
model {
   vector[N] y;
   for (n in 1:N1) y[n] <- Y1[n];
   for (n in 1:N2) y[N1 + n] <- y2[n];

   y ~ multi_normal(Mu, Cov);
}

• 20行目: hyper parameterはとして与えています。これらの値とXの各値からCovはdeterministicに求まりますのでtransformed dataブロックで求めています。
• 17行目: 多変量正規分布のもう一つのパラメータのMuの方はゼロベクトルに設定することが多いです。もちろんparametricにすることもできますが、Stanでは全体の収束がとても悪くなりました。Hinton研のIain Murrayさんのプレゼン資料ではGPを使う際には機械学習の時と同様にドメイン知識を使ってあらかじめデータに対して標準化（中央化&スケーリング）することをすすめています。とが1ぐらいになるようにするのが良いみたいです。これらによってMuをゼロベクトルのままで実行したいところです。

キックするRコードは以下の通りです。

library(rstan)
library(ggplot2)

d <- read.delim("input/GPbook-Fig2_05.data.txt", sep="\t")
N1 <- nrow(d)
x2 <- seq(from = -8, to = 8, by = 0.1)
N2 <- length(x2)

data <- list(
   N1 = N1,
   X1 = d$X,
   Y1 = d$Y,
   N2 = N2,
   X2 = x2
)

fit <- stan(
   file = 'model/GP1.stan',
   data = data,
   iter = 2000,
   warmup = 200,
   seed = 123,
   chains = 3
)


# plot
la <- extract(fit)
N.mcmc <- length(la$lp__)
N.y2 <- length(x2)

y2.med <- apply(la$y2, 2, median)
y2.low <- apply(la$y2, 2, quantile, probs = 0.05)
y2.high <- apply(la$y2, 2, quantile, probs = 0.95)

d.est <- data.frame(x = x2, y = y2.med, ymax = y2.high, ymin = y2.low)
d.obs <- data.frame(x = d$X, y = d$Y, ymax = rep(0,N1), ymin = rep(0,N1))
p <- ggplot(d.est, aes(x = x, y = y, ymax = ymax, ymin = ymin))
p <- p + xlab("input, x") + ylab("output, y")
p <- p + ylim(-4.5, 3)
p <- p + theme(legend.position = "none")
p <- p + geom_ribbon(alpha = 1/3)
p <- p + geom_line(aes(y = y), color　=　"blue")
p <- p + geom_point(data = d.obs, aes(shape = as.factor(1)))
p <- p + scale_shape_manual(values = c(1))
ggsave(file = "output/GPbook-Fig2_05-1.png", plot = p, dpi = 300, width = 4, height = 3)

• 6行目: 予測したい入力のX2を作っています。

計算には1chainあたり28sほどかかりました。推定結果の図（GPbook-Fig2_05-1.png）は以下になりました。

GP本のFig2.5とほぼ同じデータを使ったのですが、図を比べるとこちらの方が95%信用区間が広めになっています。うーん。 次回は高速化およびフルベイズを行います。