「ガウス過程と機械学習」第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のソースコードの表示方法が良くわかりませんでした。
おしまい。
randomForestで有効な交互作用を発見したい
背景
GLMは発想がわかりやすく解釈性も高くて良いアルゴリズム1なのですが、線形の仮定があるため変数間の交互作用を見るのが苦手です。実際のプロジェクトでGLMを使った結果を見せ、
- 変数の組み合わせ効果みたいなものは見れないの？
- この変数は条件によって効き方が違うんだよね〜
みたいな指摘を受けて困った経験があったりしないでしょうか。そんな時に使えるテクニックを同僚から教えてもらったので、備忘がてらメモしておきます。勝手に公開して怒られる可能性もありますが。。。
概要
手法の概要ですが、話としてはシンプルで「もしも有効な変数の組み合わせ（交互作用）が存在しているのであれば、Random Forestの各決定木において、ノードの分岐に使われる変数の順番として出現しやすいペアがあるのではないか」ということです。例えば変数X1とX2の間に交互作用があれば、決定木においてX1が選択された場合、続く分岐ではX2が選択されやすくなるのではないでしょうか。
実装
上記のアイディアを実現するために、以下のように実装してみます：
- Random Forestでモデルを作る
- 各決定木から分岐に用いられた変数ペアを得る
- 出現回数のカウントを取る
- 交互作用効果を確かめてみる
1. Random Forestでモデルを作る
まずはRandom Forestでモデルを作ります。randomForestパッケージを使ってサクッと作りましょう。
### libraryの読み込み library(randomForest) library(tidyverse)
データには前回記事と同じTelco Customer Churnを使いますが、前回の反省を踏まえてread.csvを使います。
前回の記事はこちら。
d <- read.csv("./Data/WA_Fn-UseC_-Telco-Customer-Churn.csv") %>% as_data_frame()
本来ならここから一つ一つの変数を観察するところですが、今回はそれが目的ではないので欠損だけ埋めておきます。
> colSums(apply(d, c(1, 2), is.na)) customerID gender SeniorCitizen Partner Dependents tenure 0 0 0 0 0 0 PhoneService MultipleLines InternetService OnlineSecurity OnlineBackup DeviceProtection 0 0 0 0 0 0 TechSupport StreamingTV StreamingMovies Contract PaperlessBilling PaymentMethod 0 0 0 0 0 0 MonthlyCharges TotalCharges Churn 0 11 0
TotalChargesに欠損があるようですね。
> summary(d$TotalCharges) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's 18.8 401.4 1397.5 2283.3 3794.7 8684.8 11
MedianとMeanに差があるので分布が偏っていそうです。ひとまずNAはMedianで埋めておきましょう。
d2 <- d %>% mutate(TotalCharges = if_else(is.na(.$TotalCharges), median(.$TotalCharges, na.rm = T), .$TotalCharges))
customerIDは変数として使えないので除外しましょう。またrandomForestはカテゴリ数が53より多い変数を扱えないので、カテゴリ数をチェックしておきます。
cat_vars <- sapply(d2, is.factor)
> apply(d2[, cat_vars], 2, function(x) length(unique(x))) customerID gender Partner Dependents PhoneService MultipleLines 7043 2 2 2 2 3 InternetService OnlineSecurity OnlineBackup DeviceProtection TechSupport StreamingTV 3 3 3 3 3 3 StreamingMovies Contract PaperlessBilling PaymentMethod Churn 3 3 2 4 2
大丈夫そうですね。customerIDだけ落としておきます。
d3 <- d2 %>% select(-customerID)
randomForestを当てはめます。目的変数はChurnです。
set.seed(123) result <- randomForest(Churn ~ ., d3, ntree = 500)
> result Call: randomForest(formula = Churn ~ ., data = d3, ntree = 500) Type of random forest: classification Number of trees: 500 No. of variables tried at each split: 4 OOB estimate of error rate: 20.3% Confusion matrix: No Yes class.error No 4643 531 0.1026285 Yes 899 970 0.4810059
精度とかは気にしません。
2. 各決定木から分岐に用いられた変数ペアを得る
次にRandom Forestから変数ペアを取得します。そのためにはRandom Forestの各決定木について、どの変数がどの順番で分岐に用いられたかを知る必要があります。
まずはRandom Forestから各決定木の結果を取ってきましょう。 getTree 関数を使います。過去記事も参考にしてください。
> getTree(result, 1, labelVar = TRUE) left daughter right daughter split var split point status prediction 1 2 3 Contract 1.000 1 <NA> 2 4 5 OnlineBackup 1.000 1 <NA> 3 6 7 TechSupport 2.000 1 <NA> 4 8 9 tenure 2.500 1 <NA> 5 10 11 InternetService 2.000 1 <NA> 6 12 13 Dependents 1.000 1 <NA> # 中略
getTree は分岐に用いられた変数とその分岐先ノードなどを返します。ここで必要なのは1~3列目なのですが、それぞれ「左の子ノード」「右の子ノード」「分岐に用いられた変数」を意味しています。例えば1行目を見るとノードの分岐に Contract が用いられたことがわかります。
このテーブルに行インデックスを列として追加しましょう。なおこれ以降のコードはこちらの記事を参考にさせて頂きました。
tree_tbl <- getTree(result, 1, labelVar = TRUE) %>% # labelVar = Fだとエラー rownames_to_column() %>% mutate(rowname = as.integer(rowname))
作成した tree_tbl では分岐に用いた変数（ split var ）はわかりますが、変数のペアはわかりません。例えば1行目の Contract で分岐された子ノード（2と3）は、次にそれぞれが OnlineBackup と TechSupport で分岐されています（2行目、3行目）ので、 Contract - OnlineBackup と Contract - TechSupport という変数ペアが出現したことがわかるような形に整形したいですね。
各行には「分岐に用いた変数」と「分岐先の子ノードの番号」がありますので、「分岐先の子ノード（左右両方）」に「分岐元の変数」を追加すれば欲しいものが得られそうです。
まずはノードと変数のマスタを用意しましょう。
var_name <- tree_tbl %>% select(rowname, "split var") %>% rename(split_var =`split var`) %>% # スペースを`_`に修正 unique() %>% filter(!is.na(.$split_var))
続けて left daughter と right daughter それぞれに rowname でJOINします。
> tree_tbl %>% + left_join(var_name, by = c("left daughter" = "rowname")) %>% + left_join(var_name, by = c("right daughter" = "rowname")) %>% + select(rowname, `split var`, `split_var.x`, `split_var.y`) %>% + na.omit() rowname split var split_var.x split_var.y 1 1 Contract OnlineBackup TechSupport 2 2 OnlineBackup tenure InternetService 3 3 TechSupport Dependents StreamingMovies 4 4 tenure InternetService InternetService 5 5 InternetService Dependents OnlineSecurity 6 6 Dependents TotalCharges tenure # 中略
良さそうですね。ただしこのままでは後の工程で使いにくいのでもう少し加工します。
> tree_tbl %>% + left_join(var_name, by = c("left daughter" = "rowname")) %>% + left_join(var_name, by = c("right daughter" = "rowname")) %>% + select(`split var`, `split_var.x`, `split_var.y`) %>% + na.omit() %>% + rename(from_var = `split var`, + left = `split_var.x`, + right = `split_var.y`) %>% + gather(key = node, value = to_var, -from_var) %>% + select(-node) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges # 中略
これで必要なアウトプットが得られました。あとは上記の加工を関数化しておき、各決定木に当てはめれば良さそうです。
get_var_pairs <- function(tree_num, rf = result) { # 決定木の結果を得る tree_tbl <- getTree(rf, tree_num, labelVar = TRUE) %>% rownames_to_column() %>% mutate(rowname = as.integer(rowname)) var_name <- tree_tbl %>% select(rowname, "split var") %>% rename(split_var =`split var`) %>% # スペースを`_`に修正 unique() %>% filter(!is.na(.$split_var)) out <- tree_tbl %>% left_join(var_name, by = c("left daughter" = "rowname")) %>% left_join(var_name, by = c("right daughter" = "rowname")) %>% select(`split var`, `split_var.x`, `split_var.y`) %>% na.omit() %>% rename(from_var = `split var`, left = `split_var.x`, right = `split_var.y`) %>% gather(key = node, value = to_var, -from_var) %>% select(-node) return(out) }
試してみましょう。
> get_var_pairs(5, result) from_var to_var 1 PaymentMethod Contract 2 Contract TechSupport 3 OnlineSecurity tenure 4 TechSupport PaperlessBilling 5 tenure TechSupport 6 TotalCharges StreamingTV # 中略
これを全ての決定木に当てはめます。 purrr の map_dfr を使ってみましょう。
var_pairs <- map_dfr(as.matrix(1:5), get_var_pairs, result) # 少しだけ実行
> dim(var_pairs) [1] 2950 2 > head(var_pairs) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges
素直に for で書いた時と同じ結果になっていますでしょうか？
tmp <- c() for (i in 1:5) { tmp <- bind_rows(tmp, get_var_pairs(i, result)) }
> dim(tmp) [1] 2950 2 > head(tmp) from_var to_var 1 Contract OnlineBackup 2 OnlineBackup tenure 3 TechSupport Dependents 4 tenure InternetService 5 InternetService Dependents 6 Dependents TotalCharges
合っているようなので全ての結果を取得します。 randomForest では作成する決定木の数を ntree で指定するので 1:ntree で全ての決定木に適用できるのですが、直接 ntree を取ってくることは出来ないようなので length(treesize()) を使います。エラーは気にしないことにします。
var_pairs <- map_dfr(as.matrix(1:length(treesize(result))), get_var_pairs, result)
30万弱の変数ペアが得られました。
3. 出現回数のカウントを取る
さっそく変数ペアのカウントを取ってみましょう。
> var_pairs %>% + group_by(from_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 19 x 2 from_var cnt <chr> <int> 1 TotalCharges 36908 2 MonthlyCharges 36592 3 tenure 34132 4 PaymentMethod 21244 5 gender 17626 6 MultipleLines 14496 # 中略
分岐の始点となった変数（分岐元）の数は19ですが、分析対象のデータセットには目的変数を含んで20列だったので、分岐元とならない変数はなかったようです。一方で分岐元としての出現頻度には大きなばらつきがあり、 TotalCharges 、 MonthlyCharges 、 tenure が選ばれやすいようですね。
ちなみに varImpPlot で変数重要度を見てみると、これらはいずれも上位に付けており、4位以下と大きな隔たりがあるようです。
varImpPlot(result)
続いて分岐の終点となった変数（分岐先）についても見てみましょう。
> var_pairs %>% + group_by(to_var) %>% + # group_by(from_var, to_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 19 x 2 to_var cnt <chr> <int> 1 TotalCharges 47553 2 MonthlyCharges 47137 3 tenure 37964 4 PaymentMethod 20138 5 gender 14242 6 Partner 12567 # 中略
同じく19変数ありますので、全ての変数は分岐元・分岐先ともに出現しています。分岐元と同じく出現頻度はばらつきがあり、出現しやすい変数としては、 TotalCharges 、 MonthlyCharges 、 tenure となっています。これは少し意外ですね。てっきり分岐元に選ばれやすい変数と分岐先に選ばれやすい変数は違うものになると思っていましたが。
せっかくなので分岐元と分岐先で選ばれやすさが異なるか、可視化してみましょう。
plt <- var_pairs %>% group_by(from_var) %>% summarise(cnt_f = n()) %>% left_join(var_pairs %>% group_by(to_var) %>% summarise(cnt_t = n()), by = c("from_var" = "to_var")) %>% gather(var_type, cnt, -from_var) %>% rename(var = from_var) ggplot(plt, aes(x = reorder(var, -cnt), y = cnt, fill = var_type)) + geom_bar(stat = "identity", position = "dodge") + # scale_color_brewer(palette = "Set2") + # facet_wrap(~var_type, nrow = 2) + theme_classic() + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + NULL
上位3変数は分岐元として選ばれやすいですが、分岐先としては更に頻度が多くなっています。また変数間の順位にはほとんど変動はないようですね。
組み合わせでも見てみましょう。
> var_pairs %>% + group_by(from_var, to_var) %>% + summarise(cnt = n()) %>% + arrange(desc(cnt)) # A tibble: 355 x 3 # Groups: from_var [19] from_var to_var cnt <chr> <chr> <int> 1 MonthlyCharges TotalCharges 6032 2 TotalCharges MonthlyCharges 6018 3 TotalCharges TotalCharges 5858 4 MonthlyCharges MonthlyCharges 5623 5 tenure MonthlyCharges 5602 6 tenure TotalCharges 5407 # 中略
上位10件の組み合わせを見ると、9行目の PaymentMethod を除いていずれも上位3変数の組み合わせになっています。同じ変数のペアも出てきているので、レンジを絞る形で分岐条件として選ばれているようですね、なるほど。なお19 * 19 = 361なので発生していない組み合わせがあるようですが、一部ですね。
4. 交互作用効果を確かめてみる
ひとまず目的としていた分析は以上となります。今回のデータセットおよび分析条件を用いた場合、 TotalCharges 、 MonthlyCharges 、 tenure の3変数が（組み合わせの意味でも、重要度の意味でも）影響の大きい変数であるようです。したがって冒頭のようなクライアントからの指摘があった場合には、特にこの3変数を中心に他の変数との交互作用を確認していくと良いのではないでしょうか。
と、ここまで書いたところで一つ疑問が浮かんできました。このような組み合わせ効果は、GLMでも発見できないでしょうか？例えば二次の交互作用項を準備しておき、Lassoで変数選択させるとこれらの組み合わせが残らないでしょうか？
やってみましょう。yとxを用意します。
library(glmnet) # glmnetを使う y <- as.matrix(ifelse(d3$Churn == "Yes", 1, 0)) tmp <- scale(model.matrix(Churn ~ .^2 , d3)) x_vars <- which(colSums(apply(tmp, c(1,2), is.nan)) == 0) x <- cbind(1, tmp[, x_vars])
> dim(x) [1] 7043 403
Lassoは回帰係数の絶対値に対して罰則が与えられるため、説明変数のスケールを揃えておく必要があります。そのため model.matrix でダミー化したあと scale で正規化しています。またその際に分散が0であるために NaN となってしまう変数は除外しています。
このデータでLassoにかけてみましょう。まずは適切な lambda を得るために cv.glmnet を使いますが、計算時間が少しかかるため nfolds は5にしておきましょう。
res_lasso_cv <- cv.glmnet(x = x, y = y, family = "binomial", alpha = 1, nfolds = 5)
> res_lasso_cv$lambda.min [1] 0.004130735
このときにDevianceが最小となる lambda は0.004130735のようです。プロットも見ておきましょう。
plot(res_lasso_cv)
ここで縦に引かれた二本の破線はそれぞれ lambda.min および lambda.1se を表しています。 lambda.1se は lambda.min の1SD以内で最も罰則を与えたときの lambda を示すようですね。詳しくは以下を参照してください。
https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html
ひとまず lambda.min を与えたときの結果を確認しましょう。
res_lasso <- glmnet(x = x, y = y, family = "binomial", alpha = 1, lambda = res_lasso_cv$lambda.min)
# このときの回帰係数の絶対値 > as.data.frame(as.matrix(res_lasso$beta)) %>% + rownames_to_column() %>% + filter(s0 != 0) %>% + mutate(abs_beta = abs(s0)) %>% + arrange(desc(abs_beta)) %>% + select(rowname, abs_beta, s0) rowname abs_beta s0 1 tenure 7.708817e-01 -7.708817e-01 2 ContractTwo year 4.658562e-01 -4.658562e-01 3 ContractOne year 2.723315e-01 -2.723315e-01 4 InternetServiceFiber optic 2.640571e-01 2.640571e-01 5 InternetServiceNo 2.202043e-01 -2.202043e-01 6 tenure:PaymentMethodMailed check 1.714972e-01 -1.714972e-01 # 中略
むむ。。Lassoにおいても tenure は影響の大きい（回帰係数の絶対値が大きい）変数として選ばれましたが、 TotalCharges と MonthlyCharges はいませんね。
っていうか、
> as.data.frame(as.matrix(res_lasso$beta)) %>% + rownames_to_column() %>% + filter(rowname %in% c("TotalCharges", "MonthlyCharges")) rowname s0 1 MonthlyCharges 0 2 TotalCharges 0
Lassoで落とされてる。。。
TotalCharges がどのような影響を示すのか partialPlot で見てみましょう。
partialPlot(result, as.data.frame(d3), "TotalCharges") # tibbleのままだとエラーになる
グニャグニャですね。特定のレンジでは影響が大きいものの、他ではそうでもないということなんでしょうか。だからRandom Forestのような非線形のアルゴリズムだと効果が認められる一方、Lassoのような線形のアルゴリズムでは拾いきれないのかもしれません2。これは素直に、Random Forestの結果から効果のありそうな組み合わせ変数を見つけ、分布を見ながら組み込んだ方が良さそうです。
しかしy軸が1を超えるのはなぜなんでしょうか。。。
終わりに
今回の分析はRandom Forestの結果から交互作用の良い候補を見つけようという趣旨でした。また同様の結果がLassoからも得られるかを検証しましたが、両者の結果は異なるものとなりました。Random Forestは非線形な効果を捉えることができるアルゴリズムなのでこちらの結果から有効な変数ペアを絞り込み、一つずつ検証していくスタイルが良さそうです。
