RのMASSパッケージのlm.ridge関数のソースコードを見たのだが,ひと目では何をしているのかさっぱりわからなかった.

ので,解読した結果をメモしておく.

リッジ回帰

$\mathrm{X} \in \mathbb{R}^{n \times p}$を説明変数を集めた行列,$\boldsymbol{y} \in \mathbb{R}^n$ を従属変数とする.プログラム中では $\mathrm{X}$ は X に,$\boldsymbol{y}$ は Y に対応する.このとき,リッジ回帰による回帰係数は以下のように書ける.

以下のように記号を定義する.

$\mathrm{X}^m, \boldsymbol{y}^m$はそれぞれXm,Ymに相当する.

次に $\mathrm{X}$ をスケーリングした新しいデータ行列 $\tilde{\mathrm{X}}$ を以下のように定義する.

プログラム中の X <- X/rep(Xscale, rep(n, p)) の時点での X が $\tilde{\mathrm{X}}$ に相当する.

$\tilde{\mathrm{X}}$ のSVDを $\tilde{\mathrm{X}}=\mathrm{UDV}^\top$と定義する.$\mathrm{U,D,V}$はそれぞれ,$n \times p$, $p \times p$, $p \times p$行列であり,以下のように書けるとする.

このときリッジ回帰の解は以下のように書ける.

$\hat{\boldsymbol{\beta}}_{\mathrm{ridge}}$がソースコード中のcoefに相当する.

lm.ridgeのソースコード

MASSパッケージ中のlm.ridge関数のソースコードを以下に載せておく.(バージョンは7.3-29)

> lm.ridge
function (formula, data, subset, na.action, lambda = 0, model = FALSE, 
    x = FALSE, y = FALSE, contrasts = NULL, ...) 
{
    m <- match.call(expand.dots = FALSE)
    m$model <- m$x <- m$y <- m$contrasts <- m$... <- m$lambda <- NULL
    m[[1L]] <- quote(stats::model.frame)
    m <- eval.parent(m)
    Terms <- attr(m, "terms")
    Y <- model.response(m)
    X <- model.matrix(Terms, m, contrasts)
    n <- nrow(X)
    p <- ncol(X)
    offset <- model.offset(m)
    if (!is.null(offset)) 
        Y <- Y - offset
    if (Inter <- attr(Terms, "intercept")) {
        Xm <- colMeans(X[, -Inter])
        Ym <- mean(Y)
        p <- p - 1
        X <- X[, -Inter] - rep(Xm, rep(n, p))
        Y <- Y - Ym
    }
    else Ym <- Xm <- NA
    Xscale <- drop(rep(1/n, n) %*% X^2)^0.5
    X <- X/rep(Xscale, rep(n, p))
    Xs <- svd(X)
    rhs <- t(Xs$u) %*% Y
    d <- Xs$d
    lscoef <- Xs$v %*% (rhs/d)
    lsfit <- X %*% lscoef
    resid <- Y - lsfit
    s2 <- sum(resid^2)/(n - p - Inter)
    HKB <- (p - 2) * s2/sum(lscoef^2)
    LW <- (p - 2) * s2 * n/sum(lsfit^2)
    k <- length(lambda)
    dx <- length(d)
    div <- d^2 + rep(lambda, rep(dx, k))
    a <- drop(d * rhs)/div
    dim(a) <- c(dx, k)
    coef <- Xs$v %*% a
    dimnames(coef) <- list(names(Xscale), format(lambda))
    GCV <- colSums((Y - X %*% coef)^2)/(n - colSums(matrix(d^2/div, 
        dx)))^2
    res <- list(coef = drop(coef), scales = Xscale, Inter = Inter, 
        lambda = lambda, ym = Ym, xm = Xm, GCV = GCV, kHKB = HKB, 
        kLW = LW)
    class(res) <- "ridgelm"
    res
}
<bytecode: 0x450f980>
<environment: namespace:MASS>

その他

kHKBkLWはリッジ回帰のチューニングパラメータ $\lambda$ の候補である.詳しくは3つ目の参考文献の論文(とその参考文献)に書いてある.

参考文献



blog comments powered by Disqus

Published

23 September 2014

Tags