GAMをもう少し理解したい
とても久しぶりの更新です。
背景
業務でモデリングを行うとき、私は大抵の場合GLMから始めます。目的変数に合わせて柔軟に分布を選択することが可能で、回帰係数という極めて解釈性の高い結果を得ることができるというのが理由です。
一方でGLMを使っていて不満に感じることの一つが、（ の世界で）非線形な効果を表現できないという点です。もちろん２次・３次の項や交互作用項を追加することである程度それらの不満は解消できるのですが、もう少しデータからそれらの特徴を学習したいと思うことがあります。
今回取り上げる一般化加法モデル（Generalized Additive Model, GAM）は、そのような複雑な関連性を表現できるよう説明変数に非線形な変換を行うもので、GLMを拡張したものとなります。ちょっと古いですが、2015年にMicrosoft RearchがKDDに出した論文(PDF)では、GAMを指して「the gold standard for intelligibility when low-dimensional terms are considered」と言っており、解釈性を保ちつつ高い予測精度を得ることができるモデルとしています。なおこの論文ではGAMに２次までの交互作用を追加するGA2Mという手法を提案しています。
このGAMの実装について調べた内容を書き留めておきます。GAMがどういうものかとか、平滑化に関する説明は、他に良いページがありますのでそちらを参照してください。例えば：
- https://logics-of-blue.com/%E5%B9%B3%E6%BB%91%E5%8C%96%E3%82%B9%E3%83%97%E3%83%A9%E3%82%A4%E3%83%B3%E3%81%A8%E5%8A%A0%E6%B3%95%E3%83%A2%E3%83%87%E3%83%AB/
- http://testblog234wfhb.blogspot.com/2014/07/generalized-additive-model.html
- https://blog.datarobot.com/jp/2017/10/24/ga2m-and-rating-table
GAMの実行結果
始めに、GAMを使うとどのような結果を得ることができるのか確認しましょう。ちなみに検証した環境は以下の通りです。
> sessionInfo() R version 3.6.0 (2019-04-26) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Mojave 10.14.6 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib locale: [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] compiler_3.6.0 tools_3.6.0 knitr_1.23 xfun_0.7
GAMの実装としてRでは {mgcv}が使われることが多いようですが、今回は「Rによる統計的学習入門」を参考に{gam}を使用しました。ちなみに、この本は統計・機械学習の主要な手法を網羅的に押さえつつ章末にRでの実行方法が紹介されており、大変勉強になる良書です。
- 作者: Gareth James,Daniela Witten,Trevor Hastie,Robert Tibshirani,落海浩,首藤信通
- 出版社/メーカー: 朝倉書店
- 発売日: 2018/08/03
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
同書で用いているサンプルデータ Wage を使用するため、{ISLR}を同時に読み込みます。この{ISLR}パッケージは上記の本の原著であるIntroduction of Statistical Learning with Applications in Rから来ているようです。またこのデータは、アメリカの大西洋岸中央部における男性3000人の賃金、および年齢や婚姻状況、人種や学歴などの属性が記録されています。
library(gam) library(ISLR)
データの中身を見てみましょう。
> head(Wage) year age maritl race education region jobclass 231655 2006 18 1. Never Married 1. White 1. < HS Grad 2. Middle Atlantic 1. Industrial 86582 2004 24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic 2. Information 161300 2003 45 2. Married 1. White 3. Some College 2. Middle Atlantic 1. Industrial 155159 2003 43 2. Married 3. Asian 4. College Grad 2. Middle Atlantic 2. Information 11443 2005 50 4. Divorced 1. White 2. HS Grad 2. Middle Atlantic 2. Information 376662 2008 54 2. Married 1. White 4. College Grad 2. Middle Atlantic 2. Information health health_ins logwage wage 231655 1. <=Good 2. No 4.318063 75.04315 86582 2. >=Very Good 2. No 4.255273 70.47602 161300 1. <=Good 1. Yes 4.875061 130.98218 155159 2. >=Very Good 1. Yes 5.041393 154.68529 11443 1. <=Good 1. Yes 4.318063 75.04315 376662 2. >=Very Good 1. Yes 4.845098 127.11574
このデータを用いて早速フィッティングしてみましょう。 glm と同様に関数 gam でモデルを当てはめることができます。
res_gam <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
ここで s(age, 5) は、説明変数 year を平滑化した上でモデルに取り込むことを意味し、5は平滑化の自由度です。なお s() の時点ではまだ平滑化は行われておらず、平滑化に必要な情報を属性として付与しているだけのようです。
> head(s(Wage$age, 5)) [1] 18 24 45 43 50 54 attr(,"spar") [1] 1 attr(,"df") [1] 5 attr(,"call") gam.s(data[["s(Wage$age, 5)"]], z, w, spar = 1, df = 5) attr(,"class") [1] "smooth"
ではフィッティングした結果を見てみましょう。
> summary(res_gam) Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage) Deviance Residuals: Min 1Q Median 3Q Max -119.43 -19.70 -3.33 14.17 213.48 (Dispersion Parameter for gaussian family taken to be 1235.69) Null Deviance: 5222086 on 2999 degrees of freedom Residual Deviance: 3689770 on 2986 degrees of freedom AIC: 29887.75 Number of Local Scoring Iterations: 2 Anova for Parametric Effects Df Sum Sq Mean Sq F value Pr(>F) s(year, 4) 1 27162 27162 21.981 2.877e-06 *** s(age, 5) 1 195338 195338 158.081 < 2.2e-16 *** education 4 1069726 267432 216.423 < 2.2e-16 *** Residuals 2986 3689770 1236 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Anova for Nonparametric Effects Npar Df Npar F Pr(F) (Intercept) s(year, 4) 3 1.086 0.3537 s(age, 5) 4 32.380 <2e-16 *** education --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
パッと見ると、 glm を当てはめたときと同様の結果が得られるようです。実際、 （2019/10/17修正 今回のケースでは lm.wfit は使われないので正しくありませんでした。）gam の本体（と思われる） gam.fit の中では stats::lm.wfit が呼ばれます。
しかし、 glm の結果と大きく異なる点として、各説明変数の回帰係数が出ていません。これはどういうことでしょう。以下のように説明変数の効果をプロットしてみます。
par(mfrow = c(1, 3)) plot(res_gam, se = TRUE, col = "blue")
真ん中の age が分かりやすいですが、この変数は34まではy軸の値が0未満となっており、 age は（ year や education を固定した下で）35歳程度まで wage の平均を下回っていることがわかります。その後45歳をピークに減少傾向に入り、65歳を過ぎると再び平均を下回るようになり、以降は急激に減少していきます。このような現象は、年齢とともに役職が上がることで賃金が増加し、退職によって減少することを考えれば非常に納得感のあるものだと思います。
なお上のプロットに必要なデータは以下のように取れるので、説明変数の各点における目的変数に対する影響を数値でも確認できます（必要な関数がエクスポートされていないので gam::: で直接呼び出しています）。
tmp <- gam:::preplot.Gam(res_gam, terms = gam:::labels.Gam(res_gam)) age_sm <- cbind(tmp$`s(age, 5)`$x, tmp$`s(age, 5)`$y) age_sm_uniq <- unique(age_sm[order(age_sm[, 1]), ])
プロットしてみましょう。同じ線が描けます。
plot(age_sm_uniq, type = "l", xlab = "age", ylab = "s(age, 5)")
さて、上記の結果は例えば age の二次の項を含めることで lm でも再現できるかもしれません。例えば以下のようになります：
### lmでフィッティング res_lm <- lm(wage ~ year + poly(age, 2) + education, Wage) ## poly(age, 2)で二次の多項式とする ### 予測用のデータ作成。ageだけを変化させ、yearとeducationは固定する。 x_lm <- seq(min(Wage$age), max(Wage$age), length.out = 100) nd <- data.frame(year = rep(2003, 100), age = x_lm, education = rep("2. HS Grad")) prd_lm <- predict(res_lm, nd, se.fit = T) ### 予測値から平均を引いてからプロット prd_m <- 90 plot(x_lm, prd_lm$fit - prd_m, type = "l", col = "blue", ylim = c(-40, 10)) lines(x_lm, prd_lm$fit - prd_m + 2*prd_lm$se.fit, lty = "dashed") lines(x_lm, prd_lm$fit - prd_m - 2*prd_lm$se.fit, lty = "dashed")
大体同じようなプロットを作成する事ができました。しかしピークを過ぎてからの緩やかな減少は表現できていませんし、一つ一つの変数について何次までの多項式を含めるかを検討していくのは少し手間がかかります。GAMを使えばデータに存在する細やかな変化を自動的に捉えることができます（もちろん代償もあります）。
GAMの実装
gam()
それでは gam がどのようにフィッティングを行っているのかを見ていきましょう。本体は最後の方に出てくる gam.fit なのですが、途中も少し細かく追ってみます。コンソールで gam を実行すると、以下のように関数の中身を見ることができます。
まず以下のブロックでは、 gam にオプションとして指定した内容に沿った処理を実行しています。
function (formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, control = gam.control(...), model = TRUE, method = "glm.fit", x = FALSE, y = TRUE, ...) { ### 関数の引数を名前付きで確定。gam(wage ~ s(year, 4) + education, Wage)として与えた場合、 ### formula = と data = がそれぞれ保持される。 ### match.call returns a call in which all of the specified arguments are specified by their full names. call <- match.call() ### familyの判定 if (is.character(family)) family <- get(family, mode = "function", envir = parent.frame()) if (is.function(family)) family <- family() if (is.null(family$family)) { print(family) stop("`family' not recognized") } ### データが指定されていない場合 if (missing(data)) data <- environment(formula) ### 指定されている引数の取り出し mf <- match.call(expand.dots = FALSE) m <- match(c("formula", "data", "subset", "weights", "etastart", "mustart", "offset"), names(mf), 0L) mf <- mf[c(1L, m)] ### 指定されていない引数を指定し、stats::model.frame()の形式に仕立てる mf$na.action = quote(na.pass) mf$drop.unused.levels <- TRUE mf[[1L]] <- quote(stats::model.frame)
次に、ここで一つ gam ならではの処理として平滑化に使う関数を取り出しています。
### 平滑化の関数を取ってくる（s, lo, random) gam.slist <- gam.smoothers()$slist
gam.smoothers として新しい平滑化関数を指定することも出来るようですが、デフォルトは s 、 lo および random で、それぞれ 平滑化スプライン、局所回帰、ランダム効果としての指定を意味しているようです。 3つめの random がわからなかったので調べてみたところ、これはカテゴリ変数に対する指定で、パラメータの推定においていわゆる縮小推定を行なうもののようでした。
https://www.rdocumentation.org/packages/gam/versions/1.16.1/topics/random
これらに接頭語として gam を加えた（e.g. gam.s）関数が平滑化のための関数として実行されます。 gam.s については後述します。
次のブロックですが、 call クラスであった mf を評価することで data.frame に変換しています。
### term を mf$formula に渡す mt <- if (missing(data)) terms(formula, gam.slist) else terms(formula, gam.slist, data = data) mf$formula <- mt ### ここで mf 、つまり model.frame が実行されて data.frame になる ### ただし平滑化は実行されず、平滑化のパラメータは attribute として持っている mf <- eval(mf, parent.frame()) if (missing(na.action)) { naa = getOption("na.action", "na.fail") na.action = get(naa) } mf = na.action(mf) mt = attributes(mf)[["terms"]]
ここまでは mf は call クラス、すなわち未評価の関数およびその引数を要素に持つオブジェクトでした。それが eval で評価されたため stats::model.frame が実行され、 formula にしたがい data.frame が生成された、という流れのようです（間違っていたらすみません）。
### method の指定によって処理を分ける。 glm.fit または glm.fit.null 以外の場合はエラー switch(method, model.frame = return(mf), glm.fit = 1, glm.fit.null = 1, stop("invalid `method': ", method))
ここでは method によって処理を分けていますが、現状は glm.fit または model.frame のみ受け付けているようです。 なおhelpを参照すると、 model.frame を指定した場合フィッティングは行われないようですね。以下のようになります。
> gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage, method = "model.frame") wage s(year, 4) s(age, 5) education 231655 75.04315 2006 18 1. < HS Grad 86582 70.47602 2004 24 4. College Grad 161300 130.98218 2003 45 3. Some College 155159 154.68529 2003 43 4. College Grad 11443 75.04315 2005 50 2. HS Grad
この通りデータが返ってきます。
以下では、YおよびXをそれぞれ抽出し、さらにフィッティングに必要なオプションを指定しています。
### Y を取り出す Y <- model.response(mf, "any") ### X を matrix で取り出す。 ### gam を実行したときのエラーメッセージ（ `non-list contrasts argument ignored` ）はここで出ている。 contrasts の指定が良くない様子。 X <- if (!is.empty.model(mt)) model.matrix(mt, mf, contrasts) else matrix(, NROW(Y), 0) ### その他パラメータ（weights, offset, mustart, etastart） weights <- model.weights(mf) offset <- model.offset(mf) if (!is.null(weights) && any(weights < 0)) stop("Negative wts not allowed") if (!is.null(offset) && length(offset) != NROW(Y)) stop("Number of offsets is ", length(offset), ", should equal ", NROW(Y), " (number of observations)") mustart <- model.extract(mf, "mustart") etastart <- model.extract(mf, "etastart")
以上で準備が完了し、 gam.fit で当てはめを行います。
### ここが本体。 gam.fit を呼び出している fit <- gam.fit(x = X, y = Y, smooth.frame = mf, weights = weights, start = start, etastart = etastart, mustart = mustart, offset = offset, family = family, control = control)
いったん後続の部分は無視して gam.fit に移りたいと思いますが、長くなったので今回はここまでにします。
「ガウス過程と機械学習」第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の求め方については本当にこれで合っているのか自信がないため、推定できているかについての判断は保留です。ご存知の方がいらっしゃれば教えてください。
