GAMをもう少し理解したい②
前回の続きです。 gam.fit()
から。
前回記事はこちら。
GAMの実装
gam.fit()
それでは gam.fit
の中身を覗いてみましょう。
### gam.fit は x, y に加えて smooth.frame を受け取る。これは gam で作った mf で、中身は平滑化に関する情報を持った data.frame function (x, y, smooth.frame, weights = rep(1, nobs), start = NULL, etastart = NULL, mustart = NULL, offset = rep(0, nobs), family = gaussian(), control = gam.control())
始めに、受け取ったオプションをもとにデータのサイズ、 control
や family
に関するパラメータの取り出しを行います。
### 列名 ynames <- if (is.matrix(y)) dimnames(y)[[1]] else names(y) xnames <- dimnames(x)[[2]] ### データの行数、列数 nobs <- NROW(y) nvars <- ncol(x) ## その他の gam.control オプション maxit <- control$maxit bf.maxit <- control$bf.maxit epsilon <- control$epsilon bf.epsilon <- control$bf.epsilon trace <- control$trace ### digits, weigths, offset digits <- -log10(epsilon) + 1 if (is.null(weights)) weights <- rep.int(1, nobs) if (is.null(offset)) offset <- rep.int(0, nobs) ### family に関するパラメータ variance <- family$variance dev.resids <- family$dev.resids aic <- family$aic linkinv <- family$linkinv mu.eta <- family$mu.eta if (!is.function(variance) || !is.function(linkinv)) stop("illegal `family' argument") valideta <- family$valideta if (is.null(valideta)) valideta <- function(eta) TRUE validmu <- family$validmu if (is.null(validmu)) validmu <- function(mu) TRUE eval(family$initialize) if (is.null(mustart)) { eval(family$initialize) } else { mukeep <- mustart eval(family$initialize) mustart <- mukeep } ### eta の初期値 eta <- if (!is.null(etastart)) etastart ### エラーチェック else if (!is.null(start)) if (length(start) != nvars) stop("Length of start should equal ", nvars, " and correspond to initial coefs for ", deparse(xnames)) else { coefold <- start offset + as.vector(if (NCOL(x) == 1) x * start else x %*% start) } else family$linkfun(mustart) ### mu の初期値 mu <- linkinv(eta) if (!(validmu(mu) && valideta(eta))) stop("Can't find valid starting values: please specify some") ### デビアンス(残差平方和) new.dev <- sum(dev.resids(y, mu, weights))
ここまでは関数に渡されたオプションを元に処理を進めていますが、以降のプロセスで gam
特有の処理が行われています。まずは smoothers
の取り出しです。
### smoothers が指定されている場合に smoother を抽出する ### 今回のケースでは s が入る a <- attributes(attr(smooth.frame, "terms")) smoothers <- a$specials
オブジェクトの属性を取り出すにあたり、 attr
は属性名を指定する必要がありますが、 attributes
では全ての属性をリスト形式で取り出します。 a
というオブジェクトには、具体的には以下のようなリストが格納されます。
> a $variables list(wage, s(year, 4), s(age, 5), education) $factors s(year, 4) s(age, 5) education wage 0 0 0 s(year, 4) 1 0 0 s(age, 5) 0 1 0 education 0 0 1 $term.labels [1] "s(year, 4)" "s(age, 5)" "education" $specials $specials$s [1] 2 3 $specials$lo NULL $specials$random NULL $order [1] 1 1 1 $intercept [1] 1 $response [1] 1 $class [1] "terms" "formula" $.Environment <environment: R_GlobalEnv> $predvars list(wage, s(year, 4), s(age, 5), education) $dataClasses wage s(year, 4) s(age, 5) education "numeric" "numeric" "numeric" "factor"
a$specials$s
には 2 3
という数字が入っていますが、これは smooth.frame
の2・3列目が平滑化の対象であるということだと思います。
今回は平滑化が含まれるので以下が処理されます:
if (length(smoothers) > 0) { ### NA ではない要素を取り出す smoothers <- smoothers[sapply(smoothers, length) > 0] ### 以下の処理はちょっとよくわからない for (i in seq(along = smoothers)) { tt <- smoothers[[i]] ff <- apply(a$factors[tt, , drop = FALSE], 2, any) smoothers[[i]] <- if (any(ff)) seq(along = ff)[a$order == 1 & ff] else NULL } }
smoothers
に格納されている数値が 1 2
になりました。
> smoothers $s [1] 1 2
次に、以下の処理ではパラメータの推定に使われるエンジンを決めます。平滑化が含まれるケースでは general.wam
または {s, lo}.wam
が、そうでなければ lm.wfit
が選ばれます。二種類以上の平滑化関数(s
, lo
, random
)が指定されている、または s
lo
以外の平滑化関数が指定されている場合には general.wam
が選ばれるようです。
条件にしたがってエンジンを決めたあとは、後ろの工程で評価するための bf.call
を作成しています。
### smoother が指定されている場合、ここが処理される if (length(smoothers) > 0) { gam.wlist = gam.smoothers()$wlist smooth.labels <- a$term.labels[unlist(smoothers)] assignx <- attr(x, "assign") assignx <- assign.list(assignx, a$term.labels) which <- assignx[smooth.labels] ### 2つ以上の smoothers が指定されている場合は general.wam が指定される ### wam は weighted additive model if (length(smoothers) > 1) bf <- "general.wam" ### 今回のケース(平滑化の関数として s を指定)はこちらで、 s.wam が指定される。 else { sbf <- match(names(smoothers), gam.wlist, FALSE) bf <- if (sbf) paste(gam.wlist[sbf], "wam", sep = ".") else "general.wam" } ### bf にオプション部分を文字列で結合 bf.call <- parse(text = paste(bf, "(x, z, wz, fit$smooth, which, fit$smooth.frame,bf.maxit,bf.epsilon, trace)", sep = ""))[[1]] s <- matrix(0, length(y), length(which)) dimnames(s) <- list(names(y), names(which)) fit <- list(smooth = s, smooth.frame = smooth.frame) } ### 平滑化しない場合、通常の lm.wfit に渡す。なお general.wam でも lm.wfit が使われる。 else { bf.call <- expression(lm.wfit(x, z, wz, method = "qr", singular.ok = TRUE)) bf <- "lm.wfit" }
なお上記コードブロックの最後に lm.wfit
が使われていますが、LMやGLMの中身がどうなっているのかについては過去記事を参照してください。
bf.call
が作られたので、反復しながら解を求めにいきます。ここからが本体となるブロックです。
### デビアンス(ループの打ち切り基準) old.dev <- 10 * new.dev + 10 ### ここから反復に入る for (iter in 1:maxit) { ### weight が 0 のデータは除外する good <- weights > 0 varmu <- variance(mu) if (any(is.na(varmu[good]))) stop("NAs in V(mu)") if (any(varmu[good] == 0)) stop("0s in V(mu)") mu.eta.val <- mu.eta(eta) if (any(is.na(mu.eta.val[good]))) stop("NAs in d(mu)/d(eta)") good <- (weights > 0) & (mu.eta.val != 0) ### z を生成。ただし bf.call で二番目の引数は y なので、この z は目的変数の意味 z <- eta - offset z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good] ### wz を生成。重み。今回のケースでは無視して良い wz <- weights wz[!good] <- 0 wz[good] <- wz[good] * mu.eta.val[good]^2/varmu[good]
上のブロックでは weights
を元に分析対象を定義し、 z
および wz
を生成しています。
ところで、今回のケースでは bf.call
は以下のようになるのでした:
> bf.call s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, bf.epsilon, trace)
z
は第二引数となっていますが、後ほど紹介するようにこの s.wam
という関数の中では二番目の引数は y
なので、この z
は目的変数となります(そのものではありませんが)。
続いて以下の処理が行われますが、これが gam.fit
の本体となります:
### ここで bf.call が評価される。 s.wam の場合、bakfit が呼ばれる。 ### bf.call で指定されている smooth.frame はこの時点では単なる data.frame fit <- eval(bf.call)
eval
によって bf.call
が評価されるので、 s.wam
が実行されます。 s.wam
を見に行くまえに、このループ内での残りの処理を簡単に見ておきましょう。
### 予測値にオフセットを加算する eta <- fit$fitted.values + offset ### eta から mu に変換する mu <- linkinv(eta) ### デビアンスを更新 old.dev <- new.dev new.dev <- sum(dev.resids(y, mu, weights)) ### ここの trace は対角和ではなく、 gam のオプションでループごとのデビアンスをモニターできる if (trace) cat("GAM ", bf, " loop ", iter, ": deviance = ", format(round(new.dev, digits)), " \n", sep = "") ### ループの打ち切り判定 ### デビアンスが NA となった場合、警告を出して打ち切る if (is.na(new.dev)) { one.more <- FALSE warning("iterations terminated prematurely because of singularities") } ### 差が十分に小さければ終了 else one.more <- abs(old.dev - new.dev)/(old.dev + 0.1) > epsilon if (!one.more) break }
評価結果を受けとり、デビアンスをもとにループを継続するかを判定しています。 gam.fit
の残りの工程は一旦無視して、 s.wam
を見に行きましょう。また次回。
GAMをもう少し理解したい
とても久しぶりの更新です。
背景
業務でモデリングを行うとき、私は大抵の場合GLMから始めます。目的変数に合わせて柔軟に分布を選択することが可能で、回帰係数という極めて解釈性の高い結果を得ることができるというのが理由です。
一方でGLMを使っていて不満に感じることの一つが、( の世界で)非線形な効果を表現できないという点です。もちろん2次・3次の項や交互作用項を追加することである程度それらの不満は解消できるのですが、もう少しデータからそれらの特徴を学習したいと思うことがあります。
今回取り上げる一般化加法モデル(Generalized Additive Model, GAM)は、そのような複雑な関連性を表現できるよう説明変数に非線形な変換を行うもので、GLMを拡張したものとなります。ちょっと古いですが、2015年にMicrosoft RearchがKDDに出した論文(PDF)では、GAMを指して「the gold standard for intelligibility when low-dimensional terms are considered」と言っており、解釈性を保ちつつ高い予測精度を得ることができるモデルとしています。なおこの論文ではGAMに2次までの交互作用を追加する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に乗じる()
- 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)
-
果たして今の時代にどれほどのニーズがあるのかわかりませんが。。。↩