PRML６章「ガウス過程による回帰」を R で試す

カーネル法わからん……
ガウス過程わからん……
そもそも今回の会場無事たどり着けるかな……

３つめの不安はとりあえず置いといて、わからんときは手を動かすしかない。
というわけで PRML 6.4.2 「ガウス過程による回帰」を R で試す。

訓練データは PRML のサンプルデータを使う。

# PRML's synthetic data set
curve_fitting <- data.frame(
	x=c(0.000000,0.111111,0.222222,0.333333,0.444444,0.555556,0.666667,0.777778,0.888889,1.000000),
	t=c(0.349486,0.830839,1.007332,0.971507,0.133066,0.166823,-0.848307,-0.445686,-0.563567,0.261502))

与えられた x_1,...,x_N について、対応する予測値を y_1,...,y_N としたとき、ガウス過程とは、y_1, ... , y_N がガウス分布に従うとしたものである。

って、いきなりちょっと待って！
y_n はどうやって求めるの！？
と焦ってしまうが、それは野暮らしい。

y_1, ... , y_N は求めるものではなくて、ガウス分布に従う確率変数であり、その分布はカーネル関数という１つの関数によって制御されるから、カーネル関数さえ決めればモデルが決まるよ、ということ(だと思う)。
PRML ３章〜５章では y_n を求める方法をなにがしか用意してそれを評価する、という手法だったわけだが、ガウス過程による回帰はまるっきり異なっていて、ちょっと戸惑う。
ん？　ガウス過程とカーネル法は独立した概念なのか？　PRML ではそこは若干曖昧だなー。

y_1, ... , y_N の事前分布は、グラム行列 K を共分散とする p(y)=N(y|0,K) とする(PRML 6.4.1)。
このとき、周辺分布 p(t) の共分散 C_N は K + β^-1 I となり、条件付き分布 p(t_N+1|t) の平均は $\bf{k}^t \bf{C}_N^{-1} \bf{t}$ となる。ただし、 $\bf{k}=(k(\bf{x}_n, \bf{x}_{N+1}))$ , (n=1,...,N)。
この平均が訓練データにフィットする関数を与える。まさに回帰。

カーネル関数が決まれば、上記計算を行うことができる。
機械学習ではもうさんざん言われることだが、モデル選択(この場合はカーネル選び)が最も重要。
6.4.2 では広く使われているという４つのパラメータを持つカーネル関数を紹介されているが、そんな難しいの、いきなりなんて使いこなせない！(苦笑)

今回は安直にガウスカーネルを使おう(ガウス過程だし！←関係ないし)。
ガウスカーネルのσはとりあえず 0.1。３章のお試しでは σ=0.1 でおおむねうまくいったので(あっちはガウス基底だけどね)。
精度もいるな……わかんないから、これも３章で使ったことのある β=25 とかにしておく。てきとー。

データが１次元であることを利用して、ちょこちょこ手抜きしつつも、R で regression function を求める実装を書いてみた。

# Gaussian kernel
sigma <- 0.1;
kernel <- function(x1, x2) exp(-(x1-x2)^2/(2*0.1^2));

# Gram matrix
K <- outer(curve_fitting$x, curve_fitting$x, kernel);

# covariance of marginal
beta <- 25
C_N <- K + diag(10)/beta

# mean of posterior = regression function
m <- function(x) (outer(x, curve_fitting$x, kernel) %*% solve(C_N) %*% curve_fitting$t)

求められた回帰関数 m のグラフを描く。

plot(m, xlim=c(0,1), ylim=c(-1,1))
par(new=T)
plot(curve_fitting, xlim=c(0,1), ylim=c(-1,1), xlab="", ylab="")

over fitting 気味ではあるけど、ちゃんと関数を推定できている。
ちんぷんかんぷんだったガウス過程の数式が、こうやって実際に値を放り込んでみれば、なるほど、少し意味を持って見えてきた。
気がする。

【追記】

σとβを変えて回帰させてみた。

f <- function(beta, sigma) {
	kernel <- function(x1, x2) exp(-(x1-x2)^2/(2*sigma^2));
	K <- outer(curve_fitting$x, curve_fitting$x, kernel);
	C_N <- K + diag(10)/beta
	m <- function(x) (outer(x, curve_fitting$x, kernel) %*% solve(C_N) %*% curve_fitting$t)

	plot(m, xlim=c(0,1), ylim=c(-1,1), xlab="", ylab="")
	par(new=T)
	plot(curve_fitting, xlim=c(0,1), ylim=c(-1,1), xlab=paste("beta=",beta,", sigma=",sigma), ylab="")
}

par(mfrow=c(2,2)) 
f(100, 0.1);
f(4, 0.1);
f(25, 0.3);
f(25, 0.03);