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
はざっくりと以下の4つの処理に分けることができそうです:
- 平滑化の対象変数について元のデータを順位に置き換え、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")) #### ループの2回目以降は序盤の処理をスキップさせる 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回目以降ではこの処理はスキップされます。
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 にして返す ### ループ2回目以降は評価されない 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
はまた次回。