Rと格闘 半日

研究で多次元尺度法を使うのだがRのソースがどうしても読めずに困ってしまってニャンニャンニャニャーンだったのが解決した。以下cmdscaleのsource

function (d, k = 2, eig = FALSE, add = FALSE, x.ret = FALSE) 
{
    browser()//一行ずつ実行する為の関数を追記して追っていた跡
    if (any(is.na(d))) 
        stop("NA values not allowed in 'd'")
    if (is.null(n <- attr(d, "Size"))) {
        if (add) 
            d <- as.matrix(d)
        x <- as.matrix(d^2)
        if ((n <- nrow(x)) != ncol(x)) 
            stop("distances must be result of 'dist' or a square matrix")
        rn <- rownames(x)
    }
    else {
        x <- matrix(0, n, n)
        if (add) 
            d0 <- x
        x[row(x) > col(x)] <- d^2
        x <- x + t(x)
        if (add) {
            d0[row(x) > col(x)] <- d
            d <- d0 + t(d0)
        }
        rn <- attr(d, "Labels")
    }
    if ((k <- as.integer(k)) > n - 1 || k < 1) 
        stop("'k' must be in {1, 2, ..  n - 1}")
    storage.mode(x) <- "double"
    .C(R_dblcen, x, as.integer(n), DUP = FALSE)
    /*R_dblcenというcのfileをコンパイルして使う部分。
     *ここが二重中心化しているのに全然気づけなかった。
     *後から考えれば名前がまんまじゃねーか！*/
    if (add) {
        i2 <- n + (i <- 1L:n)
        Z <- matrix(0, 2L * n, 2L * n)
        Z[cbind(i2, i)] <- -1
        Z[i, i2] <- -x
        Z[i2, i2] <- .C(R_dblcen, x = 2 * d, as.integer(n))$x
        e <- eigen(Z, symmetric = FALSE, only.values = TRUE)$values
        add.c <- max(Re(e))
        x <- matrix(double(n * n), n, n)
        non.diag <- row(d) != col(d)
        x[non.diag] <- (d[non.diag] + add.c)^2
    }
    e <- eigen(-x/2, symmetric = TRUE)
    ev <- e$values[1L:k]
    if (any(ev < 0)) 
        warning(gettextf("some of the first %d eigenvalues are < 0", k), domain = NA)
    points <- e$vectors[, 1L:k, drop = FALSE] %*% diag(sqrt(ev), k)
    dimnames(points) <- list(rn, NULL)
    if (eig || x.ret || add) {
        evalus <- e$values[-n]
        list(points = points, eig = if (eig) ev, x = if (x.ret) x, 
            ac = if (add) add.c else 0, GOF = sum(ev)/c(sum(abs(evalus)), 
                sum(evalus[evalus > 0])))
    }
    else points
}

ちなみにR_dblcenは自分で使う時は.C(stats:::R_dblcen,...)のように指定してやらなきゃならないようです。なぜならcmdscaleがstatsのメンバでありstatsの外からR_dblcenは呼べないからです。多分。
追記
/opt/local/var/macports/sources/rsync.macports.org/release/ports/math/R/work/R-2.10.0/src/library/stats/srcにR_dblcenがdblcen.cとして存在する。