統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

「ガウス過程と機械学習」第3章のグラフ(一部)を作図する③

前回の記事の続きです。

ushi-goroshi.hatenablog.com

今回は図3.16と3.20に挑戦します。データはサポートサイト(http://chasen.org/~daiti-m/gpbook/data/gpr.dat)から取得しました。

dat <- read.table("http://chasen.org/~daiti-m/gpbook/data/gpr.dat", sep = "\t",
                  col.names = c("x", "y"))
train <- head(dat, 5)

まずは図3.16から。取得したデータをプロットしてみます。

plot(train, ylim = c(-1, 3), xlab = "x", pch = 4, cex = 2)

f:id:ushi-goroshi:20190328141819p:plain

本ではX軸の数値が書かれていませんが、多分これで合っています。

このデータからRBFカーネルを使ってxの共分散行列を作成しますが、前回定義したRBFカーネルは入力としてスカラを想定していました。しかしxはベクトルとなることもあるため、RBFカーネルを少し修正します。

## 前回定義したRBFカーネル
rbf_knl_old <- function(x1, x2, theta1 = 1, theta2 = 1) {
   t(apply(as.matrix(x1), 1, function(x) theta1 * exp(-(x-x2)^2/theta2)))
}

## ベクトルに対応するために修正
rbf_knl <- function(x1, x2, theta1 = 1, theta2 = 1) {
   theta1 * exp(-norm(x1 - x2, "2")^2/theta2)
}

ベクトル同士の距離を計算するためにnorm関数を使いました。ここでは距離としてL2ノルム(ユークリッド距離)を使用していますが、その後に2乗していることに注意してください。

ベクトルでの入力に対応した代わりに共分散行列を一度に求めることができなくなってしまったので、get_cov_matも以下のように定義し直します:

## d次元のベクトルを要素とするx1[n_1, d]およびx2[n_2, d]について、x1_nとx2_n'の共分散を計算する
get_cov_mat <- function(x1, x2, theta1 = 1, theta2 = 1) {
   ## matrixに変換
   x1 <- as.matrix(x1)
   x2 <- as.matrix(x2)
   
   ## 行の組み合わせを作成する
   n <- nrow(x1)
   m <- nrow(x2)
   d <- ncol(x1)
   tmp <- cbind(kronecker(x1, matrix(1, m)), kronecker(matrix(1, n), x2))
   
   ret <- apply(tmp, 1, function(x, d, theta1, theta2) {
      rbf_knl(x[(1:d)], x[(d+1):ncol(tmp)], theta1, theta2)
   }, d, theta1, theta2)
   return(matrix(ret, n, m, byrow = T))
}

ここでは

  1. kronecker関数を使い、元の行列x1x2について全ての行の組み合わせを持つ(n*m, 2*d)のサイズとなる行列を新たに作成する(tmp
  2. 各行について、前半d個の要素を持つベクトルと後半d個のベクトルを抽出してカーネルを計算する(ret

という少しややこしい手順を踏んでいます。ここは本来なら

x1 <- matrix(seq(0, 1, length.out = 15), 5, 3, byrow = T)
x2 <- matrix(seq(1, 2, length.out = 12), 4, 3)
> apply(x1, 1, rbf_knl, x2[1, ]) # x2の1行目だけでなく、各行を順に渡したい!
[1] 0.005627029 0.025822249 0.089961445 0.237939323 0.477775042

このような形でapplyを使ってすっきり書きたいのですが、applyを二重で渡す方法がわからなかったため、一度すべての組み合わせを作成することにしました。

ひとまずこれでサンプルデータの共分散行列を見てみましょう。

> get_cov_mat(train$x, train$x)
             [,1]        [,2]       [,3]       [,4]         [,5]
[1,] 1.000000e+00 0.367879441 0.10539922 0.02705185 4.785117e-06
[2,] 3.678794e-01 1.000000000 0.77880078 0.44485807 1.930454e-03
[3,] 1.053992e-01 0.778800783 1.00000000 0.85214379 1.831564e-02
[4,] 2.705185e-02 0.444858066 0.85214379 1.00000000 7.730474e-02
[5,] 4.785117e-06 0.001930454 0.01831564 0.07730474 1.000000e+00

また、forを使って素直に書いた場合の結果と一致するかも確認しておきましょう:

n <- nrow(train)
K <- matrix(0, n, n)
for (i in 1:n) {
   for (j in 1:n) {
      if (i > j) next
      K[i, j] <- K[j, i] <- rbf_knl(train$x[i], train$x[j])
   }
}
> K
             [,1]        [,2]       [,3]       [,4]         [,5]
[1,] 1.000000e+00 0.367879441 0.10539922 0.02705185 4.785117e-06
[2,] 3.678794e-01 1.000000000 0.77880078 0.44485807 1.930454e-03
[3,] 1.053992e-01 0.778800783 1.00000000 0.85214379 1.831564e-02
[4,] 2.705185e-02 0.444858066 0.85214379 1.00000000 7.730474e-02
[5,] 4.785117e-06 0.001930454 0.01831564 0.07730474 1.000000e+00

合っているようです。

共分散行列が数字のままだと結果がわかりにくいのでヒートマップを作成してみます。

tmp <- get_cov_mat(train$x, train$x)
heatmap(tmp, Rowv = NA, Colv = NA, revC = T)

f:id:ushi-goroshi:20190328142133p:plain

おや、なんか色が変ですね、対角要素は同じ色になるはずなのですが。gplotsheatmap.2関数を使ってみます。

library(gplots)
heatmap.2(get_cov_mat(train$x, train$x), Rowv = NA, Colv = NA, dendrogram = "none", 
          trace = "none")

f:id:ushi-goroshi:20190328142213p:plain

こっちは良さそうですね。ではこのまま進めます。共分散行列を作成し、そこからデータをサンプリングしてみましょう:

K <- get_cov_mat(train$x, train$x)
n <- 1
mu <- rep(0, nrow(train))

set.seed(123)
y <- MASS::mvrnorm(n = n, mu = mu, Sigma = K)
plot(y, type = "l")

f:id:ushi-goroshi:20190328142255p:plain

良さそうです。

ところでこのグラフを何回かプロットしてみるとわかりますが、y[1]y[2]およびy[4]y[5]の間には他と比べて大きな差が生じやすい傾向にあります。これは、train$xの数値が離れているため共分散が小さい("似ていない")ことを反映しているのだと思います。そのためyのプロットも少し滑らかさに欠けるものとなっています。

図3.17のアルゴリズムを参考に、yの平均と分散を推定します。なおtheta1theta2およびtheta3はそれぞれ10.40.1とのことですが、作図の印象を図3.16に近づけるためにtheta30.02としました。

test <- seq(-1, 3.5, 0.05)
theta1 <- 1
theta2 <- 0.4
theta3 <- 0.02 # 本では0.1

## x * xの共分散行列を計算する
K <- get_cov_mat(train$x, train$x, theta1 = theta1, theta2 = theta2)

## 対角要素に誤差分散を加える
diag(K) <- diag(K) + theta3

## x * new_xの共分散行列を計算する
k   <- get_cov_mat(train$x, test, theta1 = theta1, theta2 = theta2)

## new_x * new_xの共分散行列を計算する
s   <- get_cov_mat(test, test, theta1 = theta1, theta2 = theta2)
diag(s)   <- diag(s) + theta3

## 平均の推定
Mu  <- t(k) %*% solve(K) %*% train$y # yyはここでまとめて計算した

## 分散の推定
Var <- s - t(k) %*% solve(K) %*% k
# Var <- s

## 信頼区間
CI_u  <- Mu + 2 * sqrt(diag(Var))
CI_l  <- Mu - 2 * sqrt(diag(Var))
plot(train, xlim = c(-0.5, 3.5), ylim = c(-1, 3), type = "n")
polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA)
points(train$x, train$y, xlim = c(-0.5, 3.5), ylim = c(-1, 3), xlab = "x", pch = 4, cex = 2)
lines(test, Mu, xlim = c(-0.5, 3.5), ylim = c(-1, 3), type = "l", ylab = "")

f:id:ushi-goroshi:20190328142511p:plain

ちょっと本のグラフとは異なりますが、同様のプロットを作成することが出来ました。

ここで大事なポイントとして、これは本にも記載されていることですが、上記の計算では一般の機械学習アルゴリズムで認められる「反復的な処理によって点推定値を求める工程」(いわゆるフィッティングや最適化)は一切行われていません。平均も分散も、行列計算によって得られています。ガウス過程回帰においては(ハイパーパラメータを除いて)学習は存在しないということを押さえておく必要があると思います。

さて、上記のプロットを見るとデータが観測されている付近では分散(グレーの帯)が小さくなっていることに気付きます。特に左端のデータと右端のデータで顕著ですね。この分散の値はVar <- s - t(k) %*% solve(K) %*% kで計算されているのですが、それぞれの対角要素をプロットしてみると以下のようになります:

par(mfrow = c(3, 1), mar = c(3, 3, 1, 1))
plot(cbind(test, diag(s)), type = "l", xlab = "", ylab = "diag(Var)")
plot(cbind(test, diag(t(k) %*% solve(K) %*% k)), type = "l", xlab = "", ylab = "diag(Var)")
plot(cbind(test, diag(Var)), type = "l", xlab = "", ylab = "diag(Var)")

f:id:ushi-goroshi:20190328142533p:plain

sの対角要素は一定の値を取り、そこからt(k) %*% solve(K) %*% k)の値を減じるためVarとはちょうど反転した形状になっていますね。このktrain$xtestの共分散なので、元のデータが存在する付近では共分散が大きくなっているはずです。確かめてみましょう:

heatmap.2(k, Rowv = NA, Colv = NA, dendrogram = "none",
          trace = "none")

f:id:ushi-goroshi:20190328142551p:plain

横軸がわかりにくいですが、左から右まで-1から3.5の範囲を取っています。train$xの各点付近で数値が大きくなっていることがわかります。この数値(を二乗してKで除したもの)をsから減じているため、Varでは元のデータが存在する付近において分散が小さくなっているということがわかります。

ちなみに元のデータが存在する点について、前後5ポイントずつのVarの値を取ってくると以下のような感じになります:

> data.frame(
+    "x1" = diag(Var)[(-5:5)+which(round(seq(-1, 3.5, 0.05), 3) == train$x[1])],
+    "x2" = diag(Var)[(-5:5)+which(round(seq(-1, 3.5, 0.05), 3) == train$x[2])],
+    "x3" = diag(Var)[(-5:5)+which(round(seq(-1, 3.5, 0.05), 3) == train$x[3])],
+    "x4" = diag(Var)[(-5:5)+which(round(seq(-1, 3.5, 0.05), 3) == train$x[4])],
+    "x5" = diag(Var)[(-5:5)+which(round(seq(-1, 3.5, 0.05), 3) == train$x[5])]
+ )
           x1         x2         x3         x4         x5
1  0.29935028 0.20437153 0.06164572 0.04390402 0.30258031
2  0.21438419 0.14757055 0.05877749 0.04509549 0.21725814
3  0.14172533 0.10076425 0.05308955 0.04409054 0.14389957
4  0.08615993 0.06698715 0.04664421 0.04129083 0.08741468
5  0.05140771 0.04700753 0.04148571 0.03860137 0.05178516
6  0.03960411 0.03938423 0.03894127 0.03922129 0.03960784
7  0.05093298 0.04093911 0.03921663 0.04706448 0.05178599
8  0.08347275 0.04755360 0.04140554 0.06596218 0.08742076
9  0.13328870 0.05511537 0.04390402 0.09885069 0.14392206
10 0.19476888 0.06039537 0.04509549 0.14713939 0.21732025
11 0.26116772 0.06164572 0.04409054 0.21039485 0.30272710

元のデータが存在する点(6行目)が相対的に小さくなっているのがわかります。


続いて図3.20に取り掛かります。先ほどのガウス過程回帰ではハイパーパラメータであるtheta1theta2theta3に事前に決めた数値を入力していましたが、今度はこれをデータから推定しましょう。

その前にrbf_knlについて引数theta1theta2expを取ってから計算するように修正します。これは本の脚注に記載のある通り、パラメータの取りうる値域に制約がある場合に、その制約を回避するためのテクニックです。なお以降の処理で最適化にはoptimを使用しますが、その際にmethodとしてL-BFGS-Bを指定することでパラメータに対して矩形制約を与えることも可能です(が、今回はパスします)。

## theta > 0の制約を外すためにexpを付ける
rbf_knl_exp <- function(x1, x2, theta1 = 1, theta2 = 1) {
   exp(theta1) * exp(-norm(x1 - x2, "2")^2/exp(theta2))
}

関数を新しく定義したのでget_cov_matも新たに定義しましょう。

get_cov_mat_exp <- function(x1, x2, theta1 = 1, theta2 = 1) {
   
   ## matrixに変換
   x1 <- as.matrix(x1)
   x2 <- as.matrix(x2)
   
   ## 行の組み合わせを作成する
   n <- nrow(x1)
   m <- nrow(x2)
   d <- ncol(x1)
   tmp <- cbind(kronecker(x1, matrix(1, m)), kronecker(matrix(1, n), x2))
   
   ret <- apply(tmp, 1, function(x, d, theta1, theta2) {
      rbf_knl_exp(x[(1:d)], x[(d+1):ncol(tmp)], theta1, theta2) # rbf_knl_expを使う
   }, d, theta1, theta2)
   return(matrix(ret, n, m, byrow = T))
}

続いて、本の式(3.92)に従い対数尤度関数を定義します:

L <- function(param, x, y) {
   # C <- -nrow(train)/2*log(2*pi)
   C <- 0
   K <- get_cov_mat_exp(x, x, theta1 = param[1], theta2 = param[2])
   diag(K) <- diag(K) + exp(param[3])
   
   return(-log(det(K)) - t(as.matrix(y)) %*% solve(K) %*% as.matrix(y) + C)
}

ここでCは定数部です。本にしたがってCを定義すると上のCになると思うのですが、結果の対数尤度が今いち合わなかったので一旦0を与えています。パラメータの推定には影響しないはずです。

ではoptimを使って推定してみましょう。最大化したいのでfnscale = -1としています。

param <- c(1, 1, 1)
res_par_01 <- 
   optim(par = optim(par = param, fn = L, x = train$x, y = train$y, 
                     control = list(fnscale = -1))$par,
         fn = L, x = train$x, y = train$y, control = list(fnscale = -1))
> print(exp(res_par_01$par), digits = 3)
[1] 1.596 6.560 0.082

> L(res_par_01$par, train$x, train$y)
         [,1]
[1,] -1.73877

theta1theta2theta3は小数点以下3桁まで合っていますね!ただし対数尤度が少し違います(本では-1.788)。。。

optimは初期値によって影響を受けることがあるので、値を変えてもう一度。

param <- c(0.1, 0.1, 0.1)
res_par_02 <- 
   optim(par = optim(par = param, fn = L, x = train$x, y = train$y, 
                     control = list(fnscale = -1))$par,
         fn = L, x = train$x, y = train$y, control = list(fnscale = -1))
> print(exp(res_par_02$par), digits = 3)
[1] 1.597 6.560 0.082

> L(res_par_02$par, train$x, train$y)
         [,1]
[1,] -1.73877

良さそうです。

ここで得られたハイパーパラメータを使ってプロットしてみましょう。これは図3.20の(a)に当たります。

gen_gpreg_plot <- function(param, x, y) {
   
   test <- seq(-1, 3.5, 0.05)
   
   theta1 <- exp(param[1])
   theta2 <- exp(param[2])
   theta3 <- exp(param[3])
   
   K       <- get_cov_mat(x, x, theta1 = theta1, theta2 = theta2)
   diag(K) <- diag(K) + theta3
   k       <- get_cov_mat(x, test, theta1 = theta1, theta2 = theta2)
   s       <- get_cov_mat(test, test, theta1 = theta1, theta2 = theta2)
   diag(s) <- diag(s) + theta3
   Mu      <- t(k) %*% solve(K) %*% y
   Var     <- s - t(k) %*% solve(K) %*% k
   CI_u    <- Mu + 2 * sqrt(diag(Var))
   CI_l    <- Mu - 2 * sqrt(diag(Var))
   
   plot(x, y, xlim = c(-0.5, 3.5), ylim = c(-1, 3), type = "n")
   polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA)
   points(x, y, xlim = c(-0.5, 3.5), ylim = c(-1, 3), xlab = "x", pch = 4, cex = 2)
   lines(test, Mu, xlim = c(-0.5, 3.5), ylim = c(-1, 3), type = "l", ylab = "")
   
}
gen_gpreg_plot(res_par_02$par, train$x, train$y)

f:id:ushi-goroshi:20190328143010p:plain

良さそうです。最後に、datのデータを全て使って推定してみましょう。このようなデータになります。

plot(dat, ylim = c(-1, 3), xlab = "x", pch = 4, cex = 2)

f:id:ushi-goroshi:20190328143029p:plain

param <- c(1, 1, 1)
res_par_03 <- 
   optim(par = optim(par = param, fn = L, x = dat$x, y = dat$y, 
                     control = list(fnscale = -1))$par,
         fn = L, x = dat$x, y = dat$y, control = list(fnscale = -1))
> print(exp(res_par_03$par), digits = 3)
[1] 1.524 0.689 0.067
> L(res_par_03$par, dat$x, dat$y)
          [,1]
[1,] -2.509299

theta1が0.001ずれましたがパラメータとしてはうまく推定できているようです。ただし対数尤度はやっぱり一致しません。なぜだ。。。

推定されたパラメータを使ってプロットを描いてみましょう(図3.20の(b))。

gen_gpreg_plot(res_par_03$par, dat$x, dat$y)

f:id:ushi-goroshi:20190328143117p:plain

出来ました!