swh:1:snp:3a44eb759780145deea094ac2a25c5049546a085
Tip revision: ef9d905377e5c3b729fb66e8733c6fe57c61a6da authored by Han Lin Shang on 15 March 2010, 00:00:00 UTC
version 1.6
version 1.6
Tip revision: ef9d905
rapca.r
rapca = function (x, FUN = Qn, order = 4, mean = TRUE)
{
if (order < 1)
stop("Order must be positive")
X <- t(x)
n <- nrow(X)
p <- ncol(X)
if (mean) {
med <- colMeans(X)
xx <- sweep(X, 2, med)
}
else xx <- X
tmp <- La.svd(xx)
r = sum(tmp$d > (max(n, p) * max(tmp$d) * 1e-12))
P <- t(tmp$vt)[, 1:r]
tmp2 <- rstep(t(xx %*% P), order = order, r = tmp$r, mean = mean)
tmp <- P %*% tmp2$basis
if (mean) {
med <- c(med + tmp[, 1])
xx <- sweep(X, 2, med)
basis <- cbind(med, tmp[, (1:order) + 1])
coef <- cbind(rep(1, n), xx %*% basis[, -1])
}
else {
basis <- tmp
coef <- xx %*% basis
}
return(list(basis = basis, coeff = coef, X = xx))
}