「ガウス過程と機械学習」第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を使うなどしており、ちょっと理解に時間がかかりそうなので一旦ここで終わりにしておきます。
「ガウス過程と機械学習」第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)
出来ました！
「ガウス過程と機械学習」第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) }
出来ました！周期カーネルを見てみると、ちゃんと周期的な曲線が色々生成されているのがよくわかります。
今回はここまです。色々と試行錯誤しながらなので間違いがあるかもしれませんが、その際はご指摘頂けると嬉しいです。
「ガウス過程と機械学習」第3章のグラフ（一部）を作図する
※ 3/13 一部修正
いま「ガウス過程と機械学習」を読んでいるのですが、第3章のグラフが非常に面白かったので自分でも作図してみたところちょっと腹落ちするものがありました。せっかくなので記事にしておきます。今回グラフを作成したのは図3.7と3.8（P69、70）の2つです。
ガウス過程と機械学習 (機械学習プロフェッショナルシリーズ)
- 作者: 持橋大地,大羽成征
- 出版社/メーカー: 講談社
- 発売日: 2019/03/09
- メディア: 単行本（ソフトカバー）
- この商品を含むブログを見る
順番は前後してしまうのですが、図3.8から作図します。まずは必要な関数を定義しておきましょう。本にしたがい、共分散行列を計算するRBFカーネル（ガウスカーネル）のtheta1とtheta2はともに1にしてあります。
## 1~4までの数列を適当な刻み幅で生成する generate_vector <- function(by, from = 1, to = 4) seq(by, from = from, to = to) ## xからRBFカーネルを使って共分散行列を作成する get_cov_mat <- function(x, theta1 = 1, theta2 = 1) { vec <- x theta1 * exp(-(apply(as.matrix(x), 1, function(x) x-vec))^2/theta2) }
※ 3/13 apply中のfunctionの引数を修正しました
上記の関数を使って1から4までの数列xを生成し、その共分散行列Kから生成されるyとのプロットを描く関数を以下のように定義します。なるべく本のグラフと近づけるために、共分散行列のヒートマップはcorrplotパッケージのcorrplot関数を使い、method = "shade"を指定しました。
## xからKとyを生成し、そのプロットを描く gen_y_plot <- function(by) { x <- generate_vector(by) K <- get_cov_mat(x) mu <- rep(0, length(x)) y <- MASS::mvrnorm(n = 1, mu = mu, Sigma = K) par(mfcol = c(1, 2)) plot(x, y, xlim = c(0, 5)) corrplot::corrplot(K, method = "shade") }
数列の刻み幅を変えながらいくつかプロットしてみます。
set.seed(123) gen_y_plot(1) gen_y_plot(0.5) gen_y_plot(0.2) gen_y_plot(0.1)
本にある通り、刻み幅を小さくすることで順に滑らかな曲線になっていくのがよくわかります。
続いて図3.7を作成します。先ほどとほぼ同じです。
gen_x_plot <- function(by) { x <- generate_vector(by, -4, 4) K <- get_cov_mat(x) mu <- rep(0, length(x)) y <- MASS::mvrnorm(n = 1, mu = mu, Sigma = K) plot(x, y, type = "l", ylim = c(-4, 4)) }
plotを重ねるためにpar(new = T)を使っているのですが、一回目だけは作図を0から始める必要があるので処理を分けています。もっと良い方法がある気がしていますが、ひとまず欲しいグラフは得られたので良しとします。
set.seed(123) for (i in 1:4) { if(i == 1) { gen_x_plot(0.05) next } par(new = T) gen_x_plot(0.05) }
自分で書いてみると理解できているかの確認になりますし、パラメータを変えたときの挙動を見ることで理解が深まりますので本当にオススメです。この後もなるべく自分で作図してみようと思います。
データを小集団に分割しながら線形回帰の解を推定する
背景
突然ですが、一般に線形回帰と言えば以下の正規方程式：
※ 7/3 上の式を修正
をbについて解くことで得られると教わり、そのまま理解していることが多いのではないでしょうか。
これ自体は決して間違っていないのですが、実装においては計算精度の問題から、逆行列ではなくQR分解を活用して解いている場合があります。例えばRでも、以前の記事においてlmのソースコードをたどっていった結果、ハウスホルダー変換によってQR分解が行われていることを確認しました。
過去の記事はこちら。 ushi-goroshi.hatenablog.com
ここでlmと逆行列およびQR分解による解の推定値をちょっと見てみましょう。適当にデータを作成します。
set.seed(123) n <- 100 b <- c(1, 1.5) # 切片と回帰係数 x <- cbind(1, rnorm(n)) y <- x %*% b + rnorm(n)
また、それぞれによる解の推定方法を以下のように定義します。
## lm(.fit)を使う my_lm <- function() { coef(lm.fit(x, y)) } ## 逆行列で解く my_solve <- function() { solve(crossprod(x, x)) %*% crossprod(x, y) } ## QR分解で解く my_qr <- function() { solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y }
上で定義した関数は、いずれも同じ解を返します：
> cbind(my_solve(), my_qr(), my_lm()) [,1] [,2] [,3] x1 0.8971969 0.8971969 0.8971969 x2 1.4475284 1.4475284 1.4475284
一緒の値になっていますね。少し脱線しますが、ついでに計算時間も見てみましょう：
time_1000 <- data.frame(microbenchmark::microbenchmark(my_solve(), my_qr(), my_lm(), times = 1000))
> library(ggplot2) > ggplot(time_1000, aes(x = expr, y = log(time), group = expr)) + geom_violin() + coord_flip() + labs(x = "functions") + NULL
逆行列を用いた場合が一番早く、QR分解を用いたものが最も遅いようでした。なおこのグラフは横軸が対数となっていることに注意してください。
さて、このようにして線形回帰の解はQR分解を使って求めることができますが、実は計算を工夫することで、Xを小集団に分割した上でそれぞれのデータからX全体の解を得ることができます。これが何を意味するかというと、メモリに全部載せきれないような大きいデータであっても解を推定したり、あるいは線形回帰であっても並列に計算を回すことができる、ということです1。
もともと今回の記事を書こうと思ったのは、以前に「線形回帰はデータを分割して並列計算できる」という話を知人から聞いたことをふと思い出したのがきっかけです。当時は何を言っているのか今いち理解できなかったのですが、大変わかりやすい下記の記事を見つけたため、写経した内容をメモしておきます。
freakonometrics.hypotheses.org
手順
実装に取り掛かる前に手順について簡単に理解しておきましょう。まずXをQR分解すると、冒頭に示した正規方程式から得られるは以下のようになります：
QR分解によって得られる行列Qは直交行列であるため、となります。またここで積の逆行列は
という性質があることから、
となります。すなわちQR分解によって得られた行列Rの逆行列と、行列Qの転置があれば良いことになります。先ほどmy_qrを定義したときは説明なく示しましたが、これは下のように書けます：
## my_qrの定義（再掲） solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y
問題は、このおよび
をどのようにして小集団から再構成するか、ということになりますが、これは以下の手順で計算できるようです：
- 共通処理
- X、yをそれぞれ小集団に分割する
- 各小集団のXをQR分解する
を計算する
を計算する
- 1-2で得られたQを2-2で得たQに乗じる（
- 1-2で得られたQを2-2で得たQに乗じる（
- 2と3の結果およびyにより解を得る
- 3-1で得たQ'にyを乗じる
- 両者を乗じる
なおこの手順で確かにや
が再構成できることは確認できたのですが、これがなぜ上手くいくのかについては残念ながら調べても分かりませんでした。もしご存知でしたら誰か教えてください。
実装
それでは実装に入りますが、先にデータをすべて使った時の回帰係数を確認しておきましょう。サンプルデータにはcarsを使い、目的変数をdist、説明変数をspeedとした単回帰を回してみます。
> lm(dist ~ speed, data = cars)$coefficients (Intercept) speed -17.579095 3.932409
切片とspeedの回帰係数がそれぞれ-17.579、3.932と推定されました。冒頭でも確認した通り、lmの結果は下記の方法と一致します。
y <- cars$dist x <- cbind(1, cars$speed)
> cbind( + solve(crossprod(x, x)) %*% crossprod(x, y), + solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y + ) [,1] [,2] [1,] -17.579095 -17.579095 [2,] 3.932409 3.932409
この数値を、分割した小集団に対する計算結果から再び得ることが目標となります。
1. 共通処理
1. X、yをそれぞれ小集団に分割する
それではxを小集団に分割した上で解を推定していきます。今回はデータを5つに分割しましょう。xは50行のデータなので各データセットには10行ずつ割り当てられます。各データをlist形式で保存しておきます。
# 分割するデータの数 m <- 5 n_per_d <- nrow(x) / m # 割り切れなかった場合用 if (nrow(x) %% m != 0) m <- m + 1 xlist <- list() # 各xの保存用リスト ylist <- list() # 各yの保存用リスト for (i in 1:m) { if(i == m) { xlist[[i]] = x[((i-1) * n_per_d + 1):nrow(x), ] ylist[[i]] = y[((i-1) * n_per_d + 1):nrow(x)] } xlist[[i]] = x[(i-1) * n_per_d + 1:n_per_d, ] ylist[[i]] = y[(i-1) * n_per_d + 1:n_per_d] }
このような形でデータが保存されます：
> head(xlist[[1]]) [,1] [,2] [1,] 1 4 [2,] 1 4 [3,] 1 7 [4,] 1 7 [5,] 1 8 [6,] 1 9
2. 各小集団のXをQR分解する
次に各小集団をQR分解し、その結果として得られる行列QおよびRをそれぞれ保存しておきましょう。リストの各要素は、更にそれぞれQとRを要素に持つリストとなります。
QR1 <- list() # 各データセットに対するQR分解の結果を保存するリスト for (i in 1:m) { QR1[[i]] <- list(Q = qr.Q(qr(xlist[[i]])), R = qr.R(qr(xlist[[i]]))) }
この時点でQR1は、10行2列の行列Qと2行2列の上三角行列Rを要素に持つリストになっています。
> str(QR1) List of 5 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -25.3 7.48 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -39.53 2.55 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -48.38 3.48 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -58.82 2.9 $ :List of 2 ..$ Q: num [1:10, 1:2] -0.316 -0.316 -0.316 -0.316 -0.316 ... ..$ R: num [1:2, 1:2] -3.16 0 -71.47 5.87
2. R^{-1}を計算する
1. 各小集団からのRを統合する
続いてQR1に格納された行列Rを、rbindで一つにまとめます。
R1 <- c() for(i in 1:m) { R1 <- rbind(R1, QR1[[i]]$R) }
2. 再度QR分解してRを得る
このR1を再度QR分解し、 その行列Rを得ます。
R2 <- qr.R(qr(R1))
3. Rの逆行列を求める（R^{-1}）
この逆行列が、当初求めようとしていたものの1つになります。
R_inv <- solve(R2)
では、このR_invがデータ全体を使って求めたを同じ値になっているかを確認してみましょう。
> R_inv [,1] [,2] [1,] 0.1414214 0.41606428 [2,] 0.0000000 -0.02701716 > solve(qr.R(qr(x))) [,1] [,2] [1,] -0.1414214 -0.41606428 [2,] 0.0000000 0.02701716
あれ？符号が反転していますね。。
3. Q'を計算する
ひとまず置いておいて、先に進みましょう。
1. 1-2で得られたQを2-2で得たQに乗じる（Q'）
先ほどR2を計算したときと同じQR分解で、今度は行列Qを得ます。
Q1 <- qr.Q(qr(R1))
さらに説明変数の数（今回は2）ごとにデータを分割します。
## 説明変数の数 p <- ncol(x) Q2list <- list() for(i in 1:m) { Q2list[[i]] <- Q1[(i-1) * p + 1:p, ] }
このQ2listに、最初にQR分解した結果の行列Q（QR1$Q）を掛け合わせます。
Q3list <- list() for(i in 1:m) { Q3list[[i]] <- QR1[[i]]$Q %*% Q2list[[i]] }
ここで得られたQ3listはデータ全体を使ってQR分解したときのになっているはずです。確認してみましょう：
> head(cbind( + do.call("rbind", Q3list), + qr.Q(qr(x)))) [,1] [,2] [,3] [,4] [1,] 0.1414214 0.3079956 -0.1414214 -0.3079956 [2,] 0.1414214 0.3079956 -0.1414214 -0.3079956 [3,] 0.1414214 0.2269442 -0.1414214 -0.2269442 [4,] 0.1414214 0.2269442 -0.1414214 -0.2269442 [5,] 0.1414214 0.1999270 -0.1414214 -0.1999270 [6,] 0.1414214 0.1729098 -0.1414214 -0.1729098
また符号が反転してますね。。。
原因はわかりませんが、も符号が反転していたので、結果的には元に戻るはずです。気にしないで進めましょう。
4. 2と3の結果およびyにより解を得る
1. 3-1で得たQ'にyを乗じる
上で計算された行列をylistと乗じ、結果を要素ごとに足し合わせます。
Vlist <- list() for(i in 1:m) { Vlist[[i]] <- t(Q3list[[i]]) %*% ylist[[i]] } sumV <- Vlist[[1]] for(i in 2:m) { sumV <- sumV + Vlist[[i]] }
2. 両者を乗じる
最後に、2-3で得たと
sumVを掛け合わせれば解が得られるはずです。どうでしょうか？
> cbind( + R_inv %*% sumV, + solve(crossprod(x, x)) %*% crossprod(x, y), + solve(qr.R(qr(x))) %*% t(qr.Q(qr(x))) %*% y + ) [,1] [,2] [,3] [1,] -17.579095 -17.579095 -17.579095 [2,] 3.932409 3.932409 3.932409
同じですね！
終わりに
今回はデータを小集団に分割しながら線形回帰の解を推定する、ということを紹介しました。今の時代にどうしても必要な知識かと言えばそんなこともありませんが、知っておくとと良いこともあるよ、ということで。
なおこの記事の冒頭で紹介したこちらのページでは、さらに「データソースが複数に分かれている」条件でも線形回帰の解が推定できることを示しています（例えばデータを格納しているサーバーが複数に分かれており、しかもデータのコピーが難しい状況を想定しているようです）。こちらはなかなか実用的なのではないでしょうか？
freakonometrics.hypotheses.org
おまけ
上記の工程をtidyに実行しようとすると、以下のようになるようです（こちらから）
library(tidyverse) X <- data_frame(intercept = 1, speed = cars$speed) %>% as.matrix() y <- cars$dist mats <- X %>% as_data_frame() %>% mutate( id = rep(1:5, each = 10) , y = y ) %>% ## this is where partitioning happens nest(-id) %>% mutate( X = map(data, ~ .x %>% select(-y) %>% as.matrix()), y = map(data, ~ .x %>% pull(y)) ) %>% ## We calculate QR decomposition for each partition independently mutate( Q2 = map(X, ~ .x %>% qr() %>% qr.Q()), R1 = map(X, ~ .x %>% qr() %>% qr.R()) ) df_collect <- mats$R1 %>% do.call(what = 'rbind', args = .) data.frame(dimension = c('rows', 'columns'), cbind(X %>% dim(), df_collect %>% dim())) ## Number of groups for nesting can be automatically inferred m2 <- dim(mats$R1[[1]])[2] ## The map-stage QR-decomposition Q1 = df_collect %>% qr %>% qr.Q R2 = df_collect %>% qr %>% qr.R ## For some reason this did not work with a `mutate` command... mats$Q1 = Q1 %>% as_data_frame() %>% mutate(id = ceiling(row_number() / m2)) %>% nest(-id) %>% mutate(data = map(data, ~ as.matrix(.x))) %>% pull(data) v_sum = mats %>% mutate(Q3_t = map2(.x = Q2, .y = Q1, .f = ~ t(.x %*% .y))) %>% mutate(V = map2(.x = Q3_t, .y = y, .f = ~ .x %*% .y)) %>% pull(V) %>% reduce(`+`) t(solve(R2) %*% v_sum)
-
果たして今の時代にどれほどのニーズがあるのかわかりませんが。。。↩
過小分散なカウントデータを扱いたい
背景
カウントデータをモデリングしようと思ったとき、まず思い浮かべる分布といえばポアソン分布だと思います。しかしポアソン分布は期待値と分散が等しいという性質があるため、実際のデータに当てはめようとすると、（ポアソン分布から期待されるものよりも）分散が大きい（過分散）または小さい（過小分散）という問題に直面することがあります。
そのようなとき、過分散ならともかく1、過小分散の場合にどんな分布を当てはめれば良いのか知らなかったのですが、先日某所でダブルポアソン分布というものを教えてもらったので試してみます。
doublepoissonを触ってみる
いきなりですが、ダブルポアソン分布で調べてみたところRではrmutilパッケージでdoublepoisson関数が使えるようなので早速インストールして使ってみます。
install.packages("rmutil") library(rmutil)
rmutilでは、runifやdnormなどと同じようにr/d/p/qを頭に付けたdoublepois関数が使えます。ダブルポアソン分布に従う乱数を生成してみましょう。Overdispersion parameterとしてsが指定できるのでひとまず小さめの値にしてみます。
set.seed(123) n <- 500 m <- 5 # 平均 s <- 2 # Overdispersion parameter d <- data.frame(y = rdoublepois(n, m, s))
> head(d) y 1 4 2 6 3 5 4 7 5 8 6 3
ぱっと見た感じでは通常のポアソンと変わりありません。このときの平均と分散は：
> mean(d$y) [1] 5.016 > var(d$y) [1] 2.452649
でした。うーん、分散がsと結構違ってますが、これ合っているんですかね。。。
ひとまず次に進めます。rpoisと比べ分散が小さいデータを生成できているか確認してみましょう。カウントごとの頻度を描き2、その上に（通常の）ポアソン分布の確率関数を引いてみます。
## プロット用に最大値を得る max_count <- max(d$y) ## 平均がmのポアソン分布における確率密度を得る nrm_poi <- dpois(0:max_count, m) ## データの分布を得る dens <- table(d$y) / nrow(d) ## プロット用に最大密度を得る。まるめの単位はよしなに。 max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 ## 欠損に対応するために0からmax_countまでのベクトルを用意して値を代入する dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens ## ヒストグラムではなく棒グラフで作図するので、座標を指定する xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) ## 棒グラフに確率密度関数を載せる plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "") rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue")
通常のポアソン分布で期待される密度（青線）と比較するとダブルポアソン分布では平均値周りの密度が高くなっており、全体としてバラツキが小さくなっているようです。
ではsの値を変えるとどうなるでしょうか？簡単に変更できるように関数化しておきます。
plot_double_poisson <- function(n, m, s) { d <- data.frame(y = rdoublepois(n, m, s)) max_count <- max(d$y) nrm_poi <- dpois(0:max_count, m) dens <- table(d$y) / nrow(d) max_dens <- ceiling(max(c(dens, nrm_poi)) * 20) / 20 dens_vec <- rep(0, length(0:max_count)) names(dens_vec) <- c(0:max_count) dens_vec[names(dens)] <- dens xleft <- -1:max_count - 0.2 xright <- -1:max_count + 0.2 ybottom <- rep(0, length(-1:max_count)) ytop <- c(0, dens_vec) plot(-1:max_count, xlim = c(-1, max_count), ylim = c(0, max_dens), type = "n", xlab = "", ylab = "", main = paste0("s = ", s), cex.main = 0.7) rect(xleft = xleft, ybottom = ybottom, xright = xright, ytop = ytop, col = "gray") lines(c(0:max_count), nrm_poi, col = "blue") }
sを1, 3, 5, 10とした場合のグラフを描いてみます。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 5, 1) plot_double_poisson(500, 5, 3) plot_double_poisson(500, 5, 5) plot_double_poisson(500, 5, 10)
おや？ sが小さい方がデータのバラツキが大きくなっていますね。sが10のときにはデータの半数以上が5となっているようです。またsが1のときにはほぼポアソン分布と一致しているように見えます3。この傾向はmを変更しても変わらないようです。
set.seed(123) par(mfrow = c(2, 2), mar = c(3, 3, 2, 2)) plot_double_poisson(500, 10, 1) plot_double_poisson(500, 10, 3) plot_double_poisson(500, 10, 5) plot_double_poisson(500, 10, 10)
想像していた結果とは異なるのですが、このOverdispersion parameterを変更することでポアソン分布では扱えなかった過小分散なカウントデータに対応できそうです。
切片を推定してみる
glmで推定してみる
それでは続いて、この過小分散なデータに対してGLMを当てはめたときの結果を確認してみましょう。改めてサンプルデータを作成しますが、十分に過小分散となるようsは10としておきましょう。
set.seed(123) n <- 500 m <- 5 s <- 10 d <- data.frame(y = rdoublepois(n, m, s))
glmを使ってフィッティングしますが、ここで切片として5が推定されれば問題ありません。またリンク関数はlogとします。
res_glm <- glm(y ~ 1, d, family = poisson("log"))
> coef(res_glm) (Intercept) 1.609838
切片は1.609と表示されていますが、逆リンク関数expで戻してみると、
> exp(coef(res_glm)) (Intercept) 5.002 > mean(d$y) [1] 5.002
一致しています。過小分散であっても切片は正しく推定されるようです。しかしポアソン分布であると仮定した場合、その期待値と分散は等しい（E(Y) = Var(Y)）はずなのですが、実際の分散を見てみると：
> var(d$y) [1] 0.4709379
大きく異なります。
dglmで推定してみる
それではこのようなデータに対してどのようにモデリングすれば良いかということなのですが、調べてみたところdglmというパッケージ・関数があるようなので使ってみましょう。
install.packages("dglm") library(dglm)
> res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log")) stats::glm.fit(x = Z, y = d, weights = rep(1/2, N), offset = doffset, でエラー: 'x' の中に NA/NaN/Inf があります
あらら、何らかのエラーで止まってしまいました。原因を探りたいところですが、ひとまずsを小さくして先に進みます。
set.seed(123) n <- 500 m <- 5 s <- 3 d <- data.frame(y = rdoublepois(n, m, s))
res_dglm <- dglm(y ~ 1, ~1, data = d, family = poisson("log"))
今度は大丈夫なようですので結果を見てみましょう。
> exp(res_dglm$coefficients[1]) (Intercept) 5.01 > mean(d$y) [1] 5.01
今度もmは推定できているようです。分散はどうでしょうか。
> exp(res_dglm$dispersion.fit$coefficients) (Intercept) 0.3372774 > var(d$y) [1] 1.685271
なんか、やっぱり分散は違いますね。。。
実はここまで触れてこなかったのですが、ダブルポアソン分布についてはこちらのブログを参考にさせて頂きました：
この記事ではダブルポアソン分布の分散として
と紹介されていますので、もしかしたら分散はmとの積を求める必要があるかもしれません。
> sqrt(exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients)) (Intercept) 1.299908 > var(d$y) [1] 1.685271
さっきよりもだいぶ近い値が得られるようになりましたが、右辺が平方根になっているのが気になります（分散なので）。取ってみましょう。
> exp(res_dglm$dispersion.fit$coefficients) * exp(res_dglm$coefficients) (Intercept) 1.68976 > var(d$y) [1] 1.685271
だいぶ近い！…のですが、分散の定義…。
回帰係数を推定してみる
glm、dglmで推定してみる
なかなか良さそうな感触が得られたので、切片だけでなく回帰係数も推定してみましょう。以下のようにサンプルデータを作成します。sは3としました。
set.seed(123) n <- 500 b <- 1.5 x <- runif(n, 3, 5) m <- b * x s <- 3 y <- rdoublepois(n, m, s) d <- data.frame(y, x)
平均mを与えたときのダブルポアソン分布を見てみましょう。
hist(d$y, col = "gray", ylim = c(0, 0.5), probability = T, main = "", xlab = "", breaks = "scott")
また平均と分散を確認すると：
> mean(d$y) [1] 5.982 > var(d$y) [1] 2.751178
このようになりました。
それぞれのデータを、glmとdglmそれぞれでフィッティングしてみます。まずはglmから。
res_glm <- glm(y ~ x, d, family = poisson("log"))
> exp(coef(res_glm)) (Intercept) x 2.057136 1.302974
回帰係数は真の値(1.5)にまぁまぁ近いものが推定されています。続いてdglmでフィッティングすると：
res_dglm <- dglm(y ~ x, ~1, data = d, family = poisson("log"))
> exp(res_dglm$coefficients) (Intercept) x 2.057136 1.302974
glmと同じ結果が得られました。
では分散を見てみましょう。mにはd$xの平均値を用いることにします。
m_hat <- exp(res_dglm$coefficients %*% c(1, mean(d$x)))
> m_hat [,1] [1,] 5.914556 > mean(m) [1] 5.985851
> exp(res_dglm$dispersion.fit$coefficients) * m_hat [,1] [1,] 1.947131
うーん、合いません（真の値は3）。。。
ちなみにこのデータをlmでフィッティングするとどうなるかと言うと：
res_lm <- lm(y ~ x, d)
> var(res_lm$residuals) [1] 1.935673
あまり変わりませんね。
optimで推定してみる
もう少しチャレンジしてみましょう。dglmの代わりにoptimで推定してみます。目的関数を以下のように定義します：
my_obj_fun <- function(b) { b1 <- b[1] b2 <- b[2] b3 <- b[3] ret <- sum(ddoublepois(d$y, exp(b1 + d$x * b2), b3, log = TRUE)) ret }
パラメータのベクトルを適当に用意して、optimにかけてみましょう。
b <- c(0, 1, 3) res_opt <- optim(b, my_obj_fun, control = list(fnscale = -1))
結果を確認します：
> exp(res_opt$par)[-3] # 回帰係数 [1] 2.050437 1.303756 > res_opt$par[3] # Overdispersion [1] 3.064928
おぉ！回帰係数はglm、dglmの結果と大体一致しています。またsの推定値は3.06となっており、これも真の値と近しいものになっています。optimいいですね。
終わりに
今回は過小分散なカウントデータを扱うためにrmutilとdglmというパッケージを使ってみました。optimによる推定は思いの外うまくいったようですが、dglmの結果は少しパッとしませんでした。ただしsの求め方については本当にこれで合っているのか自信がないため、推定できているかについての判断は保留です。ご存知の方がいらっしゃれば教えてください。
FindBestSplitを書いてみる
背景
前回、前々回の記事でrandomForestを使ってみたのですが、ソースコードを読んでいるとノードの分割においてfindbestsplitというサブルーチンが使われていることに気が付きました。このサブルーチン自体はこちらのL191に定義されているのでそれを読めばわかる（はずな）のですが、もう少しわかりやすい説明はないかなーと探してみたところ、こんな解説記事を見つけました。
http://dni-institute.in/blogs/cart-algorithm-for-decision-tree/
これによると、どうやらfindbestsplitは
というステップによって最良の閾値を探しているようです。それほど難しくなさそうなので、これをRで書いてみましょう。
実装
findbestsplitを実装するためには、以下のような関数が必要となりそうです。
- データ、説明変数を与えると閾値の候補を返す（return_threshold_values）
- データ、目的変数、説明変数、閾値を与えるとGini係数を返す（return_gini_index）
- 現在のGini係数との差分が最大となる（最良な）閾値を返す（return_best_value）
まずは１つ目から書いてみましょう。
1. データ、説明変数を与えると閾値の候補を返す関数
先ほど紹介したページでは閾値の候補を生成するための方法として
One of the common approach is to find splits /cut off point is to take middle values
との説明がありましたので、これに倣います。以下のように書いてみました：
- 説明変数列を
uniqueする sortで並び替える- 各要素について1つ前の値との差分（
diff）を取る - 差分を2で割る
sort後の列（最初の要素は除く）に加える
return_threshold_values <- function(dat, col) { uniq_vals <- sort(unique(dat[, col])) diffs <- diff(uniq_vals) / 2 thre_vals <- uniq_vals[-1] - diffs return(thre_vals) }
2. データ、目的変数、説明変数、閾値を与えるとGini係数を返す関数
続いてGini係数を求める関数を定義します。ここでは引数として閾値も与え、後ほどapplyで閾値候補をまとめて並列に処理することを考えました。Gini係数の求め方はここを参考に、以下のように書きました：
- 説明変数を所与の閾値で1/0のカテゴリに振り分ける
- 全体および各カテゴリのサンプルサイズを求める
- 目的変数 × 説明変数による混同行列の各要素の割合（の二乗）を計算する
- Gini係数を求める
return_gini_index <- function(dat, target, col, val) { d <- dat[, c(target, col)] d$cat <- ifelse(d[, col] > val, 1, 0) n0 <- nrow(d) n1 <- sum(d$cat) n2 <- n0 - n1 p11 <- (nrow(d[d$cat == 1 & d$target == 1, ]) / n1)^2 p12 <- (nrow(d[d$cat == 1 & d$target == 0, ]) / n1)^2 p21 <- (nrow(d[d$cat == 0 & d$target == 1, ]) / n2)^2 p22 <- (nrow(d[d$cat == 0 & d$target == 0, ]) / n2)^2 gini_val <- (n1/n0) * (1 - (p11 + p12)) + (n2/n0) * (1 - (p21 + p22)) return(gini_val) }
3. 現在のGini係数との差分が最大となる（最良な）閾値を返す関数
上記の処理によってある閾値におけるGini係数を求めることが出来ましたので、これをapplyで並列化します。また実際のところ必要な値はGini係数そのものではなく、現時点におけるGini係数との差分なので、それも計算しましょう。これまでに定義した関数を使って以下のように書きます：
- 閾値候補のベクトルを
listにする - 閾値候補リストを
return_gini_indexにsapplyで渡す - 現時点のGini係数を計算する
- 差分の最大値、Gini係数の最小値を得る
- 差分が最大となる閾値を得る
- 変数名と合わせて返す
差分の最大値やGini係数の最小値は別に必要ないのですが、参考のために取っておきます。
return_best_val <- function(dat, target, col) { thre_vals <- as.list(return_threshold_values(dat, col)) gini_vals <- sapply(thre_vals, return_gini_index, dat = dat, target = target, col = col) p1 <- sum(dat$target) / nrow(dat) p2 <- 1 - p1 current_gini <- 1 - ((p1)^2 + (p2)^2) max_gap <- max(current_gini - gini_vals) min_gini <- min(gini_vals) max_thre_val <- thre_vals[[which(max_gap == current_gini - gini_vals)]] return(c(col, max_thre_val, round(max_gap, 2), round(min_gini, 2))) }
では、このようにして定義した関数を実際に当てはめまめてみます。irisを使って適当に以下のようなデータを用意します。
my_iris <- iris my_iris$target <- ifelse(my_iris$Species == "setosa", 1, 0) dat <- my_iris[, -5]
まずはSepal.Lengthの最良な閾値を取得してみます。
> return_best_val(dat, "target", "Sepal.Length") [1] "Sepal.Length" "5.45" "0.3" "0.14"
これが合っているのかはわかりませんが、続いて全部の変数に同時に当てはめてみましょう。lapplyとdo.callを使います。
> do.call("rbind", + lapply(as.list(colnames(dat)[1:4]), return_best_val, dat = dat, target = "target")) [,1] [,2] [,3] [,4] [1,] "Sepal.Length" "5.45" "0.3" "0.14" [2,] "Sepal.Width" "3.35" "0.17" "0.28" [3,] "Petal.Length" "2.45" "0.44" "0" [4,] "Petal.Width" "0.8" "0.44" "0"
各変数について最良な閾値を取得することが出来たようです。ところで4列目（Gini係数の最小値）を確認すると、Petal.LengthとPetal.Widthで0になっていますが、Gini係数が0ということは完全に分離されていることを意味します。確かめてみましょう。
> plot(dat$Petal.Length, col = dat$target + 1) > abline(h = 2.45) > plot(dat$Petal.Width, col = dat$target + 1) > abline(h = 0.8)
確かに完全に分離しているようです。
終わりに
今回はrandomForestの中でも使われているfindbestsplitというアルゴリズムをRで書いてみました。実際には決定木やRandom Forestは一連の処理を再帰的に繰り返しているのですが、最も重要なポイントはこちらになるのだと思います。なおPythonのDecisionTreeClassifierでも同じようなアルゴリズムとなっているのか確認したかったのですが、Pythonのソースコードの表示方法が良くわかりませんでした。
おしまい。
