require("irlba") # prcomp convenience function x <- matrix(rnorm(200), nrow=20) p1 <- prcomp_irlba(x, n=3) p2 <- prcomp(x, tol=0.7) if (!isTRUE(all.equal(p1$sdev[1:2], p2$sdev[1:2]))) { stop("Failed basic prcomp test") } s <- summary(p1) # scaling bug identified in issue #21 normalize_signs <- function(X, Y) { for (i in 1:ncol(X)) { if (sign(X[1, i]) != sign(Y[1, i])) { Y[, i] <- -Y[, i] } } return(Y) } all.equal_pca <- function(X, Y) { Y <- normalize_signs(X, Y) return(all.equal(X, Y, check.attributes=F, tolerance=1e-4)) } set.seed(1) X <- matrix(rnorm(2000), ncol=40) M <- 5 # number of PCA components centers <- colMeans(X) sds <- apply(X, 2, sd) rms <- apply(X, 2, function(x) sqrt(sum(x^2) / (length(x) - 1))) Xc <- sweep(X, 2, centers, `-`) Xs <- sweep(X, 2, sds, `/`) Xcs <- sweep(Xc, 2, sds, `/`) Xrms <- sweep(X, 2, rms, `/`) # unscaled scaled <- FALSE centered <- FALSE pca <- prcomp(X, center=centered, scale.=scaled) sv <- svd(X) svir <- irlba(X, nv=M, nu=M) pcair <- prcomp_irlba(X, n=M, center=centered, scale.=scaled) Xpca <- predict(pca)[, 1:M] Xsvl <- sv$u[, 1:M] %*% diag(sv$d[1:M]) Xsvr <- X %*% sv$v[, 1:M] Xsvirl <- svir$u %*% diag(svir$d) Xsvirr <- X %*% svir$v Xpcair <- predict(pcair) Xpcair2 <- X %*% pcair$rotation if (! isTRUE(all.equal_pca(Xsvl, Xsvr)) && isTRUE(all.equal_pca(Xpca, Xsvl)) && isTRUE(all.equal_pca(Xsvirl, Xsvirr)) && isTRUE(all.equal_pca(Xpca, Xsvirl)) && isTRUE(all.equal_pca(Xpcair, Xpcair2)) && isTRUE(all.equal_pca(Xpca, Xpcair)) && isTRUE(all.equal_pca(Xpcair, Xsvirl))) { stop("failed unscaled, uncentered prcomp") } # scaled, uncentered scaled <- TRUE centered <- FALSE pca <- prcomp(X, center=centered, scale.=scaled) sv <- svd(Xrms) svir <- irlba(X, nv=M, nu=M, scale=rms) pcair <- prcomp_irlba(X, n=M, center=centered, scale.=scaled) Xpca <- predict(pca)[, 1:M] Xsvl <- sv$u[, 1:M] %*% diag(sv$d[1:M]) Xsvr <- Xrms %*% sv$v[, 1:M] Xsvirl <- svir$u %*% diag(svir$d) Xsvirr <- Xrms %*% svir$v Xpcair <- predict(pcair) Xpcair2 <- Xrms %*% pcair$rotation if (! isTRUE(all.equal_pca(Xsvl, Xsvr)) && isTRUE(all.equal_pca(Xpca, Xsvl)) && isTRUE(all.equal_pca(Xsvirl, Xsvirr)) && isTRUE(all.equal_pca(Xpca, Xsvirl)) && isTRUE(all.equal_pca(Xpcair, Xpcair2)) && isTRUE(all.equal_pca(Xpca, Xpcair)) && isTRUE(all.equal_pca(Xpcair, Xsvirl))) { stop("failed scaled, uncentered prcomp") } # issue #25 prcomp_irlba regression (error in scale. handling) set.seed(1) x <- matrix(rnorm(100), 10) p <- prcomp_irlba(x, 3, scale.=TRUE, fastpath=FALSE) p <- prcomp_irlba(x, 3, scale.=TRUE, fastpath=TRUE) # issue #32 (and also issue #25 redux) more checks for proper # variance proportion computation library(irlba) set.seed(1) x <- matrix(rnorm(200), nrow=20) n <- 3 s1 <- summary(prcomp_irlba(x, n=n, center=TRUE, scale.=FALSE)) s2 <- summary(prcomp(x, tol=0.7, center=TRUE, scale.=FALSE)) if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) && all.equal(s1$importance, s2$importance[, 1:n]))) { stop("center=TRUE scale.=FALSE prcomp variance computation") } s1 <- summary(prcomp_irlba(x, n=3, center=TRUE, scale.=TRUE)) s2 <- summary(prcomp(x, tol=0.8, center=TRUE, scale.=TRUE)) if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) && all.equal(s1$importance, s2$importance[, 1:n]))) { stop("center=TRUE scale.=TRUE prcomp variance computation") } s1 <- summary(prcomp_irlba(x, n=3, center=FALSE, scale.=TRUE)) s2 <- summary(prcomp(x, tol=0.8, center=FALSE, scale.=TRUE)) if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) && all.equal(s1$importance, s2$importance[, 1:n]))) { stop("center=FALSE scale.=TRUE prcomp variance computation") } s1 <- summary(prcomp_irlba(x, n=3, center=FALSE, scale.=FALSE)) s2 <- summary(prcomp(x, tol=0.7, center=FALSE, scale.=FALSE)) if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) && all.equal(s1$importance, s2$importance[, 1:n]))) { stop("center=FALSE, scale.=FALSE prcomp variance computation") }
Generated by dwww version 1.15 on Sat May 18 08:01:13 CEST 2024.