「ガウス過程と機械学習」第3章のグラフ（一部）を作図する④
前回の記事の続きです。
今回は図3.23に挑戦します。データはこちらから該当する部分を取ってきました。
まずはプロットしてみます。一部のデータは除外しました。
dat <- read.table("./Data/World_Record.dat", sep = "\t", col.names = c("name", "time", "date")) dat <- dat[-c(6, 17:21), -1] # ドーピング！ dat$date <- as.integer(gsub(".*, ", "", dat$date)) colnames(dat) <- c("y", "x") plot(dat$x, dat$y, xlab = "year", ylab = "time", xlim = c(1960, 2020), ylim = c(9.5, 10.1), pch = 4, cex = 2)
本で使われているデータとは少し違っている様子ですが、このまま進めます。
本にしがたい、Xおよびyを平均0、分散1にスケーリングします。scaleを使いますが、datはdata.frameで持っておきたいのでスケーリングのパラメータは別に取っておきましょう。
dat <- scale(dat) scale_pars <- unlist(attributes(dat)[c("scaled:center", "scaled:scale")]) dat <- as.data.frame(dat)
それではこのデータを使ってガウス過程回帰を実行します。まずはRBFカーネルを使って推定してみましょう。
前回の記事で定義した関数を色々使い回すことにします。まずはL、これはハイパーパラメータを与えたときに対数尤度を返す関数でした。
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) }
Lの中で使われているget_cov_mat_expも定義しましょう。これは共分散行列を得る関数です。
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)) }
さらにrbf_knl_expも使います。
rbf_knl_exp <- function(x1, x2, theta1 = 1, theta2 = 1) { exp(theta1) * exp(-norm(x1 - x2, "2")^2/exp(theta2)) }
上記の関数は、パラメータ探索における効率化のためにパラメータのlogを取ってから与えることを前提としていましたが、可視化用に元の関数も定義しておきましょう。
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)) }
rbf_knl <- function(x1, x2, theta1 = 1, theta2 = 1) { theta1 * exp(-norm(x1 - x2, "2")^2/theta2) }
ではこれらの関数を使い、今回のデータでハイパーパラメータの最適化を実行してみます。
param <- c(1, 1, 1) res_par <- 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$par), digits = 3) [1] 2.877 6.887 0.101 > L(res_par$par, dat$x, dat$y) [,1] [1,] 6.582661
教科書P94に示された値とは結構異なるようです（それぞれ1.62、0.44、0.06）。特にtheta2が大きく推定されていますね。
このパラメータを使って可視化します。やっつけですが、以下のように関数を定義しました。
gen_gpreg_plot <- function(param, x, y) { min_y <- 9.4 max_y <- 10.1 min_x <- -2 max_x <- 1.7 test <- seq(min_x, max_x, 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 * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] Var <- (s - t(k) %*% solve(K) %*% k) * (scale_pars["scaled:scale.y"])^2 CI_u <- Mu + 2 * sqrt(diag(Var)) CI_l <- Mu - 2 * sqrt(diag(Var)) x <- x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] y <- y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] test <- test * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] min_x <- min_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] max_x <- max_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] xlim <- c(min_x, max_x) ylim <- c(min_y, max_y) plot(x, y, xlim = xlim, ylim = ylim, type = "n") polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA) points(x, y, xlim = xlim, ylim = ylim, xlab = "x", pch = 4, cex = 2) lines(test, Mu, xlim = xlim, ylim = ylim, type = "l", ylab = "") }
上の関数の中で、ガウス過程によりK、kおよびsを求めるところまでは標準化後の値を使っていますが、Muを求めるところからは元のスケールに戻しています。
プロットしてみましょう。
gen_gpreg_plot(res_par$par, dat$x, dat$y)
おや、随分とスムーズな線になってしまいました。試しに教科書のパラメータを使ってみましょう。
book_par <- c(log(1.62), log(0.44), log(0.06)) gen_gpreg_plot(book_par, dat$x, dat$y)
教科書と同じようなプロットが得られました。
ガウス過程回帰の予測値（Mu）が滑らかであるということは、隣接する値が似ていることを意味しています。つまりそれらの共分散が（相対的に）大きい状態です。
RBFカーネルではtheta1とtheta2がパラメータとして与えられますが、theta2はexpの中で使われるため、特に共分散に対する影響が大きいのではないでしょうか。今回の例で言えば（expの中で分母として使われる）theta2が教科書の値よりも10倍以上大きな値となっているので、expの項が0に近づきやすかったのではないかと思います。
ちょっと確認してみましょう。theta2を教科書の値と差し替えた場合の共分散行列の変化をプロットします。
library(gplots) a <- get_cov_mat(dat$x, dat$x, res_par$par[1], res_par$par[2]) b <- get_cov_mat(dat$x, dat$x, res_par$par[1], 0.44)
heatmap.2(a, Rowv = NA, Colv = NA, dendrogram = "none", trace = "none") heatmap.2(b, Rowv = NA, Colv = NA, dendrogram = "none", trace = "none")
これを見ると、データから推定した値による共分散行列（a）はdat$xの1~5番目と6~11番目の値の共分散を大きく評価していますが、教科書の値を使った場合（b）では「似ていない」と判断しています。その結果、データから推定した方では相対的にグネグネした線が描かれたのだと思います。
続いてRBFカーネルに線形カーネルを追加して当てはめてみたいと思います。対数尤度関数を以下のように修正します：
L2 <- 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]) K <- K + exp(param[4]) * outer(x, x) # 線形カーネルを追加 return(-log(det(K)) - t(as.matrix(y)) %*% solve(K) %*% as.matrix(y) + C) }
optimによる最適化を実行しましょう。
param <- c(1, 1, 1, 1) res_par_2 <- optim(par = optim(par = param, fn = L2, x = dat$x, y = dat$y, control = list(fnscale = -1))$par, fn = L2, x = dat$x, y = dat$y, control = list(fnscale = -1))
> print(exp(res_par_2$par), digits = 3) [1] 0.147 1.234 0.104 0.772 > L2(res_par_2$par, dat$x, dat$y) [,1] [1,] 9.436611
相変わらずtheta2が教科書の値と外れますが、先ほどよりは小さな値となりました。しかし対数尤度は一致しないですね。。。
可視化用の関数を少し修正します。
gen_gpreg_plot_2 <- function(param, x, y) { min_y <- 9.4 max_y <- 10.1 min_x <- -2 max_x <- 1.7 test <- seq(min_x, max_x, 0.05) theta1 <- exp(param[1]) theta2 <- exp(param[2]) theta3 <- exp(param[3]) theta4 <- exp(param[4]) K <- get_cov_mat(x, x, theta1 = theta1, theta2 = theta2) diag(K) <- diag(K) + theta3 K <- K + theta4 * outer(x, x) k <- get_cov_mat(x, test, theta1 = theta1, theta2 = theta2) k <- k + theta4 * outer(x, test) s <- get_cov_mat(test, test, theta1 = theta1, theta2 = theta2) diag(s) <- diag(s) + theta3 s <- s + theta4 * outer(test, test) Mu <- t(k) %*% solve(K) %*% y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] Var <- (s - t(k) %*% solve(K) %*% k) * (scale_pars["scaled:scale.y"])^2 CI_u <- Mu + 2 * sqrt(diag(Var)) CI_l <- Mu - 2 * sqrt(diag(Var)) x <- x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] y <- y * scale_pars["scaled:scale.y"] + scale_pars["scaled:center.y"] test <- test * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] min_x <- min_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] max_x <- max_x * scale_pars["scaled:scale.x"] + scale_pars["scaled:center.x"] xlim <- c(min_x, max_x) ylim <- c(min_y, max_y) plot(x, y, xlim = xlim, ylim = ylim, type = "n") polygon(c(test, rev(test)), c(CI_u, rev(CI_l)), col = 'gray', border = NA) points(x, y, xlim = xlim, ylim = ylim, xlab = "x", pch = 4, cex = 2) lines(test, Mu, xlim = xlim, ylim = ylim, type = "l", ylab = "") }
プロットしてみます。
gen_gpreg_plot_2(res_par_2$par, dat$x, dat$y)
依然としてかなりスムーズな線が引かれています。教科書の値を使ってみましょう。
book_par <- c(log(0.07), log(0.02), log(0.06), log(0.92)) gen_gpreg_plot_2(book_par, dat$x, dat$y)
同じようなプロットが得られました！
以上です。 今回まで「ガウス過程と機械学習」第3章のグラフをいくつか作図してきましたが、これ以降はMCMCやGPyを使うなどしており、ちょっと理解に時間がかかりそうなので一旦ここで終わりにしておきます。
