patch to mva:prcomp to use La.svd instead of svd (PR#2227)
gregory_r_warnes@groton.pfizer.com
gregory_r_warnes@groton.pfizer.com
Tue, 29 Oct 2002 18:34:30 +0100 (MET)
Per the discussion about the problems with prcomp() when n << p, which
boils down to a problem with svd() when n << p,
here is a patch to prcomp() which substitutes La.svd() instead of svd().
-Greg
(This is really a feature enhancement, but submitted to R-bugs to make sure
it doesn't get lost. )
*** R-1.6.0/src/library/mva/R/prcomp.R Mon Aug 13 17:41:50 2001
--- R-1.6.0-GRW//src/library/mva/R/prcomp.R Tue Oct 29 11:57:23 2002
***************
*** 1,18 ****
prcomp <- function(x, retx = TRUE, center = TRUE, scale. = FALSE,
! tol = NULL) {
x <- as.matrix(x)
x <- scale(x, center = center, scale = scale.)
! s <- svd(x, nu = 0)
if (!is.null(tol)) {
rank <- sum(s$d > (s$d[1]*tol))
if (rank < ncol(x))
! s$v <- s$v[, 1:rank, drop = FALSE]
}
s$d <- s$d / sqrt(max(1, nrow(x) - 1))
! dimnames(s$v) <-
! list(colnames(x), paste("PC", seq(len = ncol(s$v)), sep = ""))
! r <- list(sdev = s$d, rotation = s$v)
! if (retx) r$x <- x %*% s$v
class(r) <- "prcomp"
r
}
--- 1,21 ----
prcomp <- function (x, retx = TRUE, center = TRUE, scale. = FALSE,
! tol = NULL)
! {
x <- as.matrix(x)
x <- scale(x, center = center, scale = scale.)
! s <- La.svd(x, nu = 0)
if (!is.null(tol)) {
rank <- sum(s$d > (s$d[1] * tol))
if (rank < ncol(x))
! s$vt <- s$vt[, 1:rank, drop = FALSE]
}
s$d <- s$d/sqrt(max(1, nrow(x) - 1))
!
! dimnames(s$vt) <- list(paste("PC", seq(len = nrow(s$vt)), sep = ""),
! colnames(x), )
! r <- list(sdev = s$d, rotation = t(s$vt) )
! if (retx)
! r$x <- x %*% t(s$vt)
class(r) <- "prcomp"
r
}
LEGAL NOTICE
Unless expressly stated otherwise, this message is confidential and may be privileged. It is intended for the addressee(s) only. Access to this E-mail by anyone else is unauthorized. If you are not an addressee, any disclosure or copying of the contents of this E-mail or any action taken (or not taken) in reliance on it is unauthorized and may be unlawful. If you are not an addressee, please inform the sender immediately.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._