GAMをもう少し理解したい③
前回の続きです。過去記事はこちらから。
ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com
GAMの実装
s.wam
s.wam は gam.fit の中で bf.call として定義されたフィッティングのための関数で、 s.wam の wam はweighted additive modelの意味のようです。名前からして gam の本体のような気がしてきますが、先に言ってしまうとこの関数内で bakfit というFortranの関数を呼び出しており、少なくともRの「 gam というライブラリの本体」という意味では合っているのではないでしょうか。
s.wam の説明に入る前に、見通しを良くするために関数全体をさっと眺めてみましょう。 s.wam はざっくりと以下の４つの処理に分けることができそうです：
- 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
- 後でFortranに渡すために必要な指定を行う
- Fortran で書かれた bakfit を呼び出す（本体）
- 後処理
### gam::s.wam function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, trace = FALSE, se = TRUE, ...) { ### 1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する if (is.data.frame(smooth.frame)) { first <- TRUE data <- smooth.frame[, names(which), drop = FALSE] smooth.frame <- gam.match(data) dx <- as.integer(dim(x)) smooth.frame$n <- dx[1] smooth.frame$p <- dx[2] oldClass(data) <- NULL smooth.frame$spar <- unlist(lapply(data, attr, "spar")) smooth.frame$df <- unlist(lapply(data, attr, "df")) } else first <- FALSE ### 2. 後でFortranに渡すために必要な指定を行う storage.mode(tol) <- "double" storage.mode(maxit) <- "integer" which <- unlist(which) storage.mode(which) <- "integer" storage.mode(y) <- "double" storage.mode(w) <- "double" p <- smooth.frame$p n <- smooth.frame$n for (ich in which) x[, ich] = signif(x[, ich], 6) ### 3. Fortran で書かれた bakfit を呼び出す（本体） fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), df = as.double(smooth.frame$df), as.integer(smooth.frame$o), as.integer(smooth.frame$nef), etal = double(n), s = s, eta = double(n), beta = double(p), var = s, tol, qr = x, qraux = double(p), qpivot = as.integer(1:p), effects = double(n), double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 15 * n + 15 + length(which)), PACKAGE = "gam") ### 4. 後処理 #### 収束に関する警告 nit <- fit$npetc[5] qrank <- fit$npetc[7] if ((nit == maxit) & maxit > 1) warning(paste("s.wam convergence not obtained in ", maxit, " iterations")) #### ループの２回目以降は序盤の処理をスキップさせる if (first) { smooth.frame$spar <- fit$spar first <- FALSE } #### 最終的な返り値となる rl に情報を追加していくための準備 names(fit$df) <- dimnames(s)[[2]] names(fit$beta) <- labels(x)[[2]] qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, pivot = fit$qpivot, tol = 1e-07), class = "qr") effects <- fit$effects r1 <- seq(len = qrx$rank) dn <- colnames(x) if (is.null(dn)) dn <- paste("x", 1:p, sep = "") names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank)) #### rl を生成 rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, fitted.values = fit$eta, effects = effects, weights = w, rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, nl.df = fit$df - 1) rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 0) if (se) rl <- c(rl, list(var = fit$var)) c(list(smooth.frame = smooth.frame), rl) }
以下、一つずつ見ていきましょう。
1. 平滑化の対象変数について元のデータを順位に置き換え、smooth.frameを再定義する
s.wam では始めに smooth.frame に対する処理を行います。具体的には
- 平滑化の対象変数を抽出
- 指定のdigitsで丸める
- ソートした上でユニークなデータ数を得る
- 元データをラベル（数値の大小で並べかえたときの順位）に変換
- 行数・列数を付与
- 平滑化パラメータおよび自由度を付与
という加工を、下記の処理にて実施します。
### smooth.frame はループ初回は data.frame なのでここが評価される if (is.data.frame(smooth.frame)) { first <- TRUE data <- smooth.frame[, names(which), drop = FALSE] # 平滑化対象変数を抽出 smooth.frame <- gam.match(data) # 下で解説 dx <- as.integer(dim(x)) # 行列のサイズ smooth.frame$n <- dx[1] # number of row smooth.frame$p <- dx[2] # number of column oldClass(data) <- NULL ### 各列に対する平滑化パラメータおよび自由度 smooth.frame$spar <- unlist(lapply(data, attr, "spar")) smooth.frame$df <- unlist(lapply(data, attr, "df")) } else first <- FALSE
上のブロックのポイントは gam.match という関数で、 gam.match(data) とすると以下の結果が返ります：
> gam.match(data) $o s(year, 4) s(age, 5) [1,] 4 1 [2,] 2 7 [3,] 1 28 [4,] 1 26 [5,] 3 33 [6,] 6 37 ### 中略 $nef s(year, 4) s(age, 5) 7 61
ちなみに data は元データです：
> head(data) s(year, 4) s(age, 5) 231655 2006 18 86582 2004 24 161300 2003 45 155159 2003 43 11443 2005 50 376662 2008 54
これは gam.match という関数において下記の部分が評価された結果です。
### gam.match から xr <- signif(as.vector(x), 6) # signif 関数は指定の有効桁で丸める。デフォルトは6 sx <- unique(sort(xr)) # それをソートしてユニークにする nef <- as.integer(length(sx)) # ユニークなデータ数 if (nef <= 3) stop("A smoothing variable encountered with 3 or less unique values; at least 4 needed") o <- match(xr, sx, nef + 1) # 元のデータを順位に変更する。 +1 は欠損用。 o[is.na(o)] <- nef + 1
match の部分が少しわかりにくいかもしれませんが、以下のような結果が得られます：
> head(xr) # 元の値 [1] 2006 2004 2003 2003 2005 2008 > match(head(xr), sx, nef + 1) # 順位に変換 [1] 4 2 1 1 3 6
これは xr の値が sx の何番目に位置するかを表わしていますが、参照先の sx は元の値をソートしてユニークにしたものなので、順位に変換していることになります。
また gam.match の結果を smooth.frame に代入されているので、この時点で smooth.frame は数値の大小関係のみを持つことになり、また一連の処理で smooth.frame は data.frame ではなくなるため、ループの２回目以降ではこの処理はスキップされます。
2. 後でFortranに渡すために必要な指定を行う
次のブロックは後ろの工程でFortranの関数に渡すための準備となります。
### 後でFortranに渡すために必要な指定 storage.mode(tol) <- "double" storage.mode(maxit) <- "integer" which <- unlist(which) storage.mode(which) <- "integer" storage.mode(y) <- "double" storage.mode(w) <- "double" p <- smooth.frame$p n <- smooth.frame$n
また以下のコードでは、先程の smooth.frame 同様に有効桁を6桁に丸めます。
### 平滑化対象の変数を6桁に丸める for (ich in which) x[, ich] = signif(x[, ich], 6)
3. Fortran で書かれた bakfit を呼び出す（本体）
以上で準備が整いました。以下が s.wam の本体になります。 .Fortran でFortranで書かれた bakfit というモジュールを呼び出します。
### ここが本体。Fortran で書かれた bakfit を呼び出す fit <- .Fortran("bakfit", x, npetc = as.integer(c(n, p, length(which), se, 0, maxit, 0)), y = y, w = w, which, spar = as.double(smooth.frame$spar), df = as.double(smooth.frame$df), as.integer(smooth.frame$o), as.integer(smooth.frame$nef), etal = double(n), s = s, eta = double(n), beta = double(p), var = s, tol, qr = x, qraux = double(p), qpivot = as.integer(1:p), effects = double(n), double((10 + 2 * 4 + 5) * (max(smooth.frame$nef) + 2) + 15 * n + 15 + length(which)), PACKAGE = "gam")
さてこの bakfit 、ソースコードとしてはおそらくこちら
が正しいのかなぁと思うのですが、 gam をインストールする際に同時にダウンロードされる以下のファイル：
/Users/hogehoge/Library/R/3.6/library/gam/ratfor/backfit.r
の方が説明もあってわかりやすいので、こちらベースに進めたいと思います。なお前回記事でも述べたことですが、 s.wam は第二引数として y を受け取ります：
### s.wam function (x, y, w, s, which, smooth.frame, maxit = 30, tol = 1e-07, trace = FALSE, se = TRUE, ...)
ただしこの s.wam にはYそのものではなく、以下のように z として渡されるのでした：
### 以下は gam.fit から再掲 z <- eta - offset z[good] <- z[good] + (y - mu)[good]/mu.eta.val[good] ### 中略 # eval(bf.call) ← これは下の表現となる s.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit, bf.epsilon, trace)
要するにリンク関数で変換したあとのデータを渡している、ということです。
では bakfit に移りたいのですが、その前にまずは残りの工程を見ておきましょう。
4. 後処理
s.wam の最後の部分です。
nit <- fit$npetc[5] # number of iteration qrank <- fit$npetc[7] # qrank ### nit が maxit と同じであった場合、収束していないことを警告 if ((nit == maxit) & maxit > 1) warning(paste("s.wam convergence not obtained in ", maxit, " iterations")) ### first == T の場合、spar を加えて FALSE にして返す ### ループ２回目以降は評価されない if (first) { smooth.frame$spar <- fit$spar first <- FALSE } names(fit$df) <- dimnames(s)[[2]] # 平滑化対象変数の列名 names(fit$beta) <- labels(x)[[2]] # 説明変数 qrx <- structure(list(qr = fit$qr, qraux = fit$qraux, rank = qrank, pivot = fit$qpivot, tol = 1e-07), class = "qr") # 属性を指定しながらオブジェクトを作成 effects <- fit$effects r1 <- seq(len = qrx$rank) dn <- colnames(x) if (is.null(dn)) dn <- paste("x", 1:p, sep = "") names(effects) <- c(dn[qrx$pivot[r1]], rep.int("", n - qrx$rank)) #### 最終的な返り値である rl を生成 rl <- list(coefficients = fit$beta, residuals = fit$y - fit$eta, fitted.values = fit$eta, effects = effects, weights = w, rank = qrank, assign = attr(x, "assign"), qr = qrx, smooth = fit$s, nl.df = fit$df - 1) rl$df.residual <- n - qrank - sum(rl$nl.df) - sum(fit$w == 0) if (se) rl <- c(rl, list(var = fit$var)) #### return c(list(smooth.frame = smooth.frame), rl)
基本的には s.wam という関数の返り値となる rl を生成するための工程と言えそうです。 bakfit はまた次回。
応用統計学フロンティアセミナー「データサイエンスと応用統計学」参加記録
10/19に開催された応用統計学フロンティアセミナーに参加してきましたので、そのメモを共有しておきます。話を聴きながらのメモなので単語しか書けず意味がわかりにくいところもありますが、ご参考まで。なおセミナーの様子は以下のtogetterでまとめられていますので、そちらも合わせてご覧ください。
- データサイエンスにおける応用統計人材の育成
- 楽天技術研究所におけるデータサイエンスおよび統計の様々な応用事例
- 自動車とデータサイエンス
- 企業ビッグデータから捉える企業活動と未来の活用可能性
- 早稲田大学におけるデータサイエンス人材育成への取り組み
- 学校における「統計」教育の課題 我々は木に縁りて魚を求めてはいないか
- 岩崎先生あいさつ
データサイエンスにおける応用統計人材の育成
横浜市立大学 岩崎学 氏
これまでとこれからの統計的データ解析
- 研究目的の設定
- データ収集法の立案（実験、観察、調査）
- データの収集（モニタリング）
- データの電子化
- データのチェック・マージ
- データの集計とグラフ化（記述統計）
- 統計的推測ないしは予測（推測統計）
- プレゼン・文章化、意思決定（終了もしくは最初に戻る） → 役に立たないと意味がない
いまのデータサイエンス、データありきになっているかもしれない
- 推測統計の部分にフォーカスが当たっている
- 全体の流れを見ることが大事
データサイエンスとは？
横浜市立大学での状況
- 線形代数の学部生の試験成績、文系と理系でほとんど差がない（120点満点で76 vs 72とか）
- 2020年、大学院
- ヘルスデータサイエンス専攻（医学部があるので）
- データサイエンス専攻（M20名、D3名）
楽天技術研究所におけるデータサイエンスおよび統計の様々な応用事例
楽天 森正弥 氏
1.3Bユーザー、48Kの店舗（ビジネスパートナー）
250Mのアイテムについて需要予測
- 価格、在庫の調整
- 非線形回帰。基盤を作ったロシア人が手作業でフィッティング。DLに勝つ。
金融市場予測
- Now-Cast
- Forecastに対してNowcast
- 今の状態を推測する
- Hal Varian
- Bayesian variable selection for nowcasting economic time series
- https://www.nber.org/chapters/c12995
- 独自の景気動向指数
- 会社の株価
- 売上と株価、２ヶ月のラグ
- 相関が0.96
- でも予測精度は芳しくない
- 化粧品会社の株価
- 多段のDoc2Vecでメーカー名のリンキング
- 楽天トラベルのデータ、インバウンド
- Now-Cast
自動車とデータサイエンス
日産自動車 上田哲郎 氏
- Connected Carのエキスパート、ITのエキスパート
- 自動車、AI/ITに関しては新卒の人たちの方が専門性高い
- 今はリアルタイムで車の車種まで特定できる
- BigData、 AI、IoT
- 箱根登って下るとバッテリー回復するので沼津まで行ける
- シミュレーションではなく過去のリーフのデータから、どこまで行けるか示せる
- ProPilot 2.0
- マジでやばい
- ハンドルから手を離して良いと言っている（高速に限る）
- ドライバーを監視している、前方を見ている限り
- Adaptive Cruise Control
- 単眼可視光カメラで前方の車との距離を測る
- ナンバープレートのサイズがわかっているのでカメラの画素数から車間距離がわかる
- AI
- オイルフィルターの真贋判定
- 空気抵抗を計算で求める
- スパコンでも数日かかる
- CADとCdの関係性
- Voxel dataとCdの関係性を学習
- 精度はまだまだ
- R値で0.7（R2のこと？）
- GANで車のデザインを生成
- 画像の生成を三次元に拡張する
- データアナリスト
- 問題のモデル化
- 人材育成
- 自動車業界、サプライヤーが持ってくる
- うまくいってない
- データを使っていて車に興味を持っている人に来てもらいたい
- 自動車がなくなっていく？
- MaaS、CASE
- 痛し痒しと思っている
- 体力のあるうちに全部やっておこう
企業ビッグデータから捉える企業活動と未来の活用可能性
帝国データバンク 中川みゆき 氏
- 帝国データバンクの保有するデータ
- 調査員が現地に足を運んで入手した、ネットでは得られないデータベース
- 100近い項目
- 定量、定性
- 中小企業が99%
- 年間60万件超の調査
- 活用事例
- 守りの活用
- 信用リスク
- 攻めの活用
- 取引関係からつながりを見る
- バリューチェーン全体のマッピング
- データはファクト、データサイエンスとは価値を創り出すこと
- 人間が本来もつ、マシンではできないことに注力
- 「未来の力」はマシンインテリジェンスに対してのリーダーシップ
- 人がやる部分
- 内訳
- 複雑な問題解決
- クリティカルシンキング
- 想像力
- マネジメント
- 人間関係調整
- 情緒的知性
- 決断
- サービスディレクション
- 交渉力
- 認識の柔軟性
- 守りの活用
早稲田大学におけるデータサイエンス人材育成への取り組み
早稲田大学 松嶋敏泰 氏
- 大隈重信と統計学
- 統計院を設立
- 統計センターの前身
- 統計院を設立
- DSとは
- Interdisciplinary Field
- concept to unify
- 日本だとSexiest Jobで広まったので、ビジネスの面が強い
- Fourth paradigm of science
- メタサイエンス
- 知覚と思考の拡大
- 早稲田のフォーカス領域
- 専門性（応用先）
- データサイエンス
- 教育
学校における「統計」教育の課題 我々は木に縁りて魚を求めてはいないか
代々木ゼミナール 西岡康夫 氏
- 初等・中等教育の現状
- サッカー、世界で活躍できるようになった理由は多くの子供がサッカーをやるようになったから
- 棟梁レベルのDSを多く育てるなら裾野を広げないといけない
- 国によってはエリート制を採用しているが、日本のカルチャーに合わないのでは
- 教育改革
- センター試験の廃止
- 混沌の極み
- 大学入学選抜は暗記、再生
- ドーリットル、甚兵衛、小ピピン
- 理系尚もて数学す、況や文系をや
- 統計が入ってきた
- 産業界の要請
- 教員研修や教職課程のカリキュラム改革
- まったく進んでない
- 統計のイロハのイもダメ
- 教員研修や教職課程のカリキュラム改革
- 産業界の要請
- 若者の強い自己否定感
- センター試験の廃止
岩崎先生あいさつ
- Technologyは変わっていくがPrincipleは変わらない
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] ### ２つ以上の 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を使っていて不満に感じることの一つが、（ の世界で）非線形な効果を表現できないという点です。もちろん２次・３次の項や交互作用項を追加することである程度それらの不満は解消できるのですが、もう少しデータからそれらの特徴を学習したいと思うことがあります。
今回取り上げる一般化加法モデル（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) }
出来ました！周期カーネルを見てみると、ちゃんと周期的な曲線が色々生成されているのがよくわかります。
今回はここまです。色々と試行錯誤しながらなので間違いがあるかもしれませんが、その際はご指摘頂けると嬉しいです。
