「ガウス過程と機械学習」第3章のグラフ（一部）を作図する②
前回の記事の続きです。
今回は図3.9と3.11（P71と75）を作成してみます。generate_vectorはそのままですが、図3.9は二次元平面なのでget_cov_matはx2を考慮できるように修正しました。またapplyについてですが、この書き方なら返り値は行ベクトルだろうとすっかり思い込んでいたら列ベクトルだったので、最後にt()で転置しています。その他、theta1とexpもapplyの中に入れました。
## 数列を適当な刻み幅で生成する generate_vector <- function(by, from = 1, to = 4) seq(by, from = from, to = to) ## x1とx2からRBFカーネルを使って共分散行列を作成する get_cov_mat <- function(x1, x2, theta1 = 1, theta2 = 1) { t(apply(as.matrix(x1), 1, function(x) theta1 * exp(-(x-x2)^2/theta2))) }
中身を見てみましょう。
x1 <- generate_vector(by = 0.5, from = 1, to = 4) x2 <- generate_vector(by = 0.5, from = 1, to = 4) K <- get_cov_mat(x1, x2)
> K [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.0000000000 0.778800783 0.36787944 0.1053992 0.01831564 0.001930454 0.0001234098 [2,] 0.7788007831 1.000000000 0.77880078 0.3678794 0.10539922 0.018315639 0.0019304541 [3,] 0.3678794412 0.778800783 1.00000000 0.7788008 0.36787944 0.105399225 0.0183156389 [4,] 0.1053992246 0.367879441 0.77880078 1.0000000 0.77880078 0.367879441 0.1053992246 [5,] 0.0183156389 0.105399225 0.36787944 0.7788008 1.00000000 0.778800783 0.3678794412 [6,] 0.0019304541 0.018315639 0.10539922 0.3678794 0.77880078 1.000000000 0.7788007831 [7,] 0.0001234098 0.001930454 0.01831564 0.1053992 0.36787944 0.778800783 1.0000000000
これだけだと合っているのか不安だったので素直にforで書いてみた結果と比較してみましょう。
n_row <- length(x1) n_col <- length(x2) tmp_mat <- matrix(0, nrow = n_row, ncol = n_col) theta1 <- theta2 <- 1 for (i in 1:n_row) { for (j in 1:n_col) { tmp_mat[i, j] <- theta1 * exp(-(x1[i] - x2[j])^2/theta2) } }
> tmp_mat [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.0000000000 0.778800783 0.36787944 0.1053992 0.01831564 0.001930454 0.0001234098 [2,] 0.7788007831 1.000000000 0.77880078 0.3678794 0.10539922 0.018315639 0.0019304541 [3,] 0.3678794412 0.778800783 1.00000000 0.7788008 0.36787944 0.105399225 0.0183156389 [4,] 0.1053992246 0.367879441 0.77880078 1.0000000 0.77880078 0.367879441 0.1053992246 [5,] 0.0183156389 0.105399225 0.36787944 0.7788008 1.00000000 0.778800783 0.3678794412 [6,] 0.0019304541 0.018315639 0.10539922 0.3678794 0.77880078 1.000000000 0.7788007831 [7,] 0.0001234098 0.001930454 0.01831564 0.1053992 0.36787944 0.778800783 1.0000000000
合っています。いくつか条件を変えながら試して同じ数値が得られたので、多分大丈夫なはずです。
ではこの共分散行列からサンプリングしてみます。まずは1回サンプリングしてみましょう。
set.seed(123) n <- 1 mu <- rep(0, n_row) y1 <- MASS::mvrnorm(n = n, mu = mu, Sigma = K)
> y1 [1] 0.8678710 0.7543409 -0.1371074 -0.3076295 0.2865098 0.7452395 1.3575729
サンプリングごとにn_row点のデータ（この場合は7点）が得られますが、このデータはそれぞれが独立ではなく、get_cov_matで生成した共分散Kにしたがっているはずです。Kは隣接する二点間で0.77、一つ飛ばすと0.36程度の相関となるような関係性でした。数値を見ても何となく隣同士は似ていて、少し飛ばすと離れているように見えますね。
ちょっと確かめてみましょう。同じようなベクトルを大量に生成して相関を取ったらどうなるでしょうか。
set.seed(123) n_rep <- 1000 tmp <- MASS::mvrnorm(n = n_rep, mu = mu, Sigma = K)
こんな感じのデータが得られると思います。
> head(tmp) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] -0.56071606 -0.5468666 0.07112351 0.7927312 0.99969485 0.9422141 0.5205130 [2,] -0.24848704 -0.5654978 -0.53304083 0.1560599 0.79484325 0.7221504 0.9147320 [3,] -0.31831368 -1.4937802 -1.53460128 -0.8919618 -0.76188288 -0.9871314 -1.3913327 [4,] 0.53077769 0.4894104 -0.69896348 -0.9056140 0.09296229 0.5327204 0.3531465 [5,] -2.76680410 -1.1244952 -0.43934824 -0.5900576 0.35234542 1.7280382 2.2408007 [6,] -0.08682172 -0.8829415 -0.87584558 -0.7091206 -1.74871286 -2.3692956 -1.3662469
相関を見てみると：
> cor(tmp) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.00000000 0.77732425 0.36703848 0.07417336 -0.03880697 -0.04802414 -0.02636937 [2,] 0.77732425 1.00000000 0.77589202 0.33492698 0.07326653 0.02051531 0.03531489 [3,] 0.36703848 0.77589202 1.00000000 0.76092428 0.35564692 0.12949585 0.06578033 [4,] 0.07417336 0.33492698 0.76092428 1.00000000 0.78343555 0.39987004 0.14778605 [5,] -0.03880697 0.07326653 0.35564692 0.78343555 1.00000000 0.79893421 0.40402673 [6,] -0.04802414 0.02051531 0.12949585 0.39987004 0.79893421 1.00000000 0.78588550 [7,] -0.02636937 0.03531489 0.06578033 0.14778605 0.40402673 0.78588550 1.00000000
合ってる！良さそうですね。
では改めて二回目のサンプリングを行います。
set.seed(456) y2 <- MASS::mvrnorm(n = n, mu = mu, Sigma = K)
ここでy1とy2は独立のサンプリングなので、図3.9におけるz方向の値はそれらの積で良いと思います（が、ここはかなり悩んだところなので正直に言えば自信ありません。間違っていたらご指摘ください）。outerを取ってプロットしてみましょう。
z <- outer(y1, y2) persp(x1, x2, z, theta = 320, phi = 20, expand = 0.4, border = "gray", shade = 0.1)
ちょっとカクカクしていますが、何となく良さそうです。ちなみにperspの引数ですが、
thetaやphi→ グラフをどの角度から見るかexpand→ グラフを縦軸方向に圧縮する比率borderやshade→ 二次元曲面の色に関する指定
となっています。この辺はお好みで。
もう少し図3.9に近づける努力をしてみます。下記を参考に色を変えつつ関数にしました。
gen_2dim_plot <- function(by, from = 1, to = 4, n = 2, theta1 = 1, theta2 = 1) { x1 <- generate_vector(by = by, from = from, to = to) x2 <- generate_vector(by = by, from = from, to = to) K <- get_cov_mat(x1, x2, theta1 = theta1, theta2 = theta2) n_row <- length(x1) mu <- rep(0, n_row) y <- MASS::mvrnorm(n = n, mu = mu, Sigma = K) z <- outer(y[1,], y[2,]) col.pal <- colorRampPalette(rev(rainbow(5))) colors <- col.pal(100) z.facet.center <- (z[-1, -1] + z[-1, -ncol(z)] + z[-nrow(z), -1] + z[-nrow(z), -ncol(z)])/4 z.facet.range <- cut(z.facet.center, 100) persp(x1, x2, z, theta = 320, phi = 20, expand = 0.4, border = "gray", shade = 0.1, col = colors[z.facet.range]) }
どうでしょう。先ほどはややカクカクしていたので、刻み幅を狭くしてみます。
set.seed(123) gen_2dim_plot(0.02, from = 0, to = 1)
良さそうですね。ちょっとパラメータを色々と変えてみましょう。
set.seed(123) gen_2dim_plot(0.05, from = 0, to = 4, theta1 = 5, theta2 = 1)
set.seed(456) gen_2dim_plot(0.05, from = 0, to = 4, theta1 = 1, theta2 = 0.5)
もちろんプロットごとにサンプリングが行われるので、同じパラメータであっても色々な曲面が生成され、こんなにも表現力があるのかと驚かされます。
続いて図3.11に挑戦してみます。ここでは様々なカーネルを当てはめるので、各カーネルを定義します。なおprd_knlについてはtheta2 = 0.5にしましたが、この理由については後述します。
## 線形カーネル lin_knl <- function(x1, x2) { outer(x1, x2) + 1 } ## RBFカーネル rbf_knl <- function(x1, x2, theta1 = 1, theta2 = 1) { t(apply(as.matrix(x1), 1, function(x) theta1 * exp(-(x-x2)^2/theta2))) } ## 指数カーネル exp_knl <- function(x1, x2, theta = 1) { exp(-apply(as.matrix(x1), 1, function(x) abs(x-x2)/theta1)) } ## 周期カーネル prd_knl <- function(x1, x2, theta1 = 1, theta2 = 0.5) { exp(theta1 * cos(apply(as.matrix(x1), 1, function(x) abs(x-x2)/theta2))) }
要はこれらを順次プロットすれば良いと思うのですが、どうせならget_cov_matの引数にカーネル名を追加して、総称関数みたいにしたいなと思い、以下のように書いてみました。
get_cov_mat <- function(kernel, x1, x2) { eval(call(kernel, x1 = x1, x2 = x2)) }
最初のKと同じものが得られていると思います。
> get_cov_mat("rbf_knl", x1, x2) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 1.0000000000 0.778800783 0.36787944 0.1053992 0.01831564 0.001930454 0.0001234098 [2,] 0.7788007831 1.000000000 0.77880078 0.3678794 0.10539922 0.018315639 0.0019304541 [3,] 0.3678794412 0.778800783 1.00000000 0.7788008 0.36787944 0.105399225 0.0183156389 [4,] 0.1053992246 0.367879441 0.77880078 1.0000000 0.77880078 0.367879441 0.1053992246 [5,] 0.0183156389 0.105399225 0.36787944 0.7788008 1.00000000 0.778800783 0.3678794412 [6,] 0.0019304541 0.018315639 0.10539922 0.3678794 0.77880078 1.000000000 0.7788007831 [7,] 0.0001234098 0.001930454 0.01831564 0.1053992 0.36787944 0.778800783 1.0000000000
本当は各カーネルでtheta1とtheta2の指定が異なるので、get_cov_matの引数に...を追加したかったのですがうまく行きませんでした。この辺はもう少しRそのものに対する知識が必要ですね。
改めてx1とx2を生成し、各カーネルによる共分散行列を作成しますが、どんな共分散が得られるのかをヒートマップで見てみましょう。ちょっと直感に合わないかもしれませんが黄色いほうが値が大きいことを示します。
x1 <- generate_vector(by = 0.1, from = -4, to = 4) x2 <- generate_vector(by = 0.1, from = -4, to = 4)
> heatmap(get_cov_mat("lin_knl", x1, x2), Rowv = NA, Colv = NA, revC = T) > heatmap(get_cov_mat("rbf_knl", x1, x2), Rowv = NA, Colv = NA, revC = T) > heatmap(get_cov_mat("exp_knl", x1, x2), Rowv = NA, Colv = NA, revC = T) > heatmap(get_cov_mat("prd_knl", x1, x2), Rowv = NA, Colv = NA, revC = T)
基本的には対角上に黄色が集まるのですが、周期カーネルはちょっと面白いヒートマップになっていますね。このような共分散行列になるのは、X1とX2の値が離れていても、ある周期に当てはまる距離であればカーネルの値は大きくなるためです（詳しくは本を確認してください）。
先ほど周期カーネルを定義したときにtheta2 = 0.5としたのは、1だとこの周期が少し大きすぎてあまり面白くなかったためです。見てみましょう。
prd_knl_1.0 <- function(x1, x2, theta1 = 1, theta2 = 1.0) { exp(theta1 * cos(apply(as.matrix(x1), 1, function(x) abs(x-x2)/theta2))) }
> heatmap(get_cov_mat("prd_knl_1.0", x1, x2), Rowv = NA, Colv = NA, revC = T)
一応周期が表現されてはいるのですが、この程度だとこの後のプロットがうまく行きませんでした。
ではこのget_cov_matを使ってyをサンプリングします。以下のような関数を作っておけば便利です。
generate_y <- function(kernel, x1, x2, n = 1) { K <- get_cov_mat(kernel, x1, x2) n_row <- length(x1) mu <- rep(0, n_row) MASS::mvrnorm(n = n, mu = mu, Sigma = K) }
図3.11に合わせて4回サンプリングし、その曲線をプロットします。
n <- 4 par(mfrow = c(2, 2)) ## 線形カーネル y <- generate_y("lin_knl", x1, x2, n = n) for (i in 1:n) { if(i == 1) { plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "f(x)", xlab = "Linear Kernel", col = i) next } par(new = T) plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "", xlab = "", col = i) } ## RBFカーネル y <- generate_y("rbf_knl", x1, x2, n = n) for (i in 1:n) { if(i == 1) { plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "f(x)", xlab = "RBF Kernel", col = i) next } par(new = T) plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "", xlab = "", col = i) } ## 指数カーネル y <- generate_y("exp_knl", x1, x2, n = n) for (i in 1:n) { if(i == 1) { plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "f(x)", xlab = "Exponential Kernel", col = i) next } par(new = T) plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "", xlab = "", col = i) } ## 周期カーネル y <- generate_y("prd_knl", x1, x2, n = n) for (i in 1:n) { if(i == 1) { plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "f(x)", xlab = "Periodic Kernel", col = i) next } par(new = T) plot(x1, y[i, ], type = "l", ylim = c(-4, 4), ylab = "", xlab = "", col = i) }
出来ました！周期カーネルを見てみると、ちゃんと周期的な曲線が色々生成されているのがよくわかります。
今回はここまです。色々と試行錯誤しながらなので間違いがあるかもしれませんが、その際はご指摘頂けると嬉しいです。
