「ガウス過程と機械学習」第3章のグラフ(一部)を作図する③
前回の記事の続きです。
今回は図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)
本では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)) }
ここでは
kronecker
関数を使い、元の行列x1
とx2
について全ての行の組み合わせを持つ(n*m, 2*d)
のサイズとなる行列を新たに作成する(tmp
)- 各行について、前半
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)
おや、なんか色が変ですね、対角要素は同じ色になるはずなのですが。gplots
のheatmap.2
関数を使ってみます。
library(gplots) heatmap.2(get_cov_mat(train$x, train$x), Rowv = NA, Colv = NA, dendrogram = "none", trace = "none")
こっちは良さそうですね。ではこのまま進めます。共分散行列を作成し、そこからデータをサンプリングしてみましょう:
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")
良さそうです。
ところでこのグラフを何回かプロットしてみるとわかりますが、y[1]
とy[2]
およびy[4]
とy[5]
の間には他と比べて大きな差が生じやすい傾向にあります。これは、train$x
の数値が離れているため共分散が小さい("似ていない")ことを反映しているのだと思います。そのためy
のプロットも少し滑らかさに欠けるものとなっています。
図3.17のアルゴリズムを参考に、y
の平均と分散を推定します。なおtheta1
、theta2
およびtheta3
はそれぞれ1
、0.4
、0.1
とのことですが、作図の印象を図3.16に近づけるためにtheta3
を0.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 = "")
ちょっと本のグラフとは異なりますが、同様のプロットを作成することが出来ました。
ここで大事なポイントとして、これは本にも記載されていることですが、上記の計算では一般の機械学習アルゴリズムで認められる「反復的な処理によって点推定値を求める工程」(いわゆるフィッティングや最適化)は一切行われていません。平均も分散も、行列計算によって得られています。ガウス過程回帰においては(ハイパーパラメータを除いて)学習は存在しないということを押さえておく必要があると思います。
さて、上記のプロットを見るとデータが観測されている付近では分散(グレーの帯)が小さくなっていることに気付きます。特に左端のデータと右端のデータで顕著ですね。この分散の値は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)")
s
の対角要素は一定の値を取り、そこからt(k) %*% solve(K) %*% k)
の値を減じるためVar
とはちょうど反転した形状になっていますね。このk
はtrain$x
とtest
の共分散なので、元のデータが存在する付近では共分散が大きくなっているはずです。確かめてみましょう:
heatmap.2(k, Rowv = NA, Colv = NA, dendrogram = "none", trace = "none")
横軸がわかりにくいですが、左から右まで-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に取り掛かります。先ほどのガウス過程回帰ではハイパーパラメータであるtheta1
、theta2
、theta3
に事前に決めた数値を入力していましたが、今度はこれをデータから推定しましょう。
その前にrbf_knl
について引数theta1
とtheta2
はexp
を取ってから計算するように修正します。これは本の脚注に記載のある通り、パラメータの取りうる値域に制約がある場合に、その制約を回避するためのテクニックです。なお以降の処理で最適化には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
theta1
、theta2
、theta3
は小数点以下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)
良さそうです。最後に、dat
のデータを全て使って推定してみましょう。このようなデータになります。
plot(dat, ylim = c(-1, 3), xlab = "x", pch = 4, cex = 2)
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)
出来ました!