library(testthat)
context("Helpers")
test_that("compute_P", {
set.seed(135435)
dat <- lfmm_sampler(n = 100, p = 1000, K = 3,
outlier.prop = 0.1,
cs = c(0.8),
sigma = 0.2,
B.sd = 1.0,
U.sd = 1.0,
V.sd = 1.0)
res <- compute_P(X = dat$X,
lambda = 1e-5)
expect_lte(mean(abs(diag(1, nrow(dat$X), nrow(dat$X))) - res$sqrt.P %*% res$sqrt.P.inv), 1e-16)
## mean(abs(res$P - res$sqrt.P %*% res$sqrt.P))
})
test_that("compute_B_ridge", {
dat <- lfmm_sampler(n = 100, p = 1000, K = 3,
outlier.prop = 0.1,
cs = c(0.8),
sigma = 0.2,
B.sd = 1.0,
U.sd = 1.0,
V.sd = 1.0)
X <- cbind(matrix(1,100,1), dat$X)
B <- compute_B_ridge(A = dat$Y,
X = X,
lambda = 1e-3)
expect_equal(dim(B), c(1000, 2))
## hist(B[,1]) ## RMK when lambda -> infinity B -> 0 OK
})
test_that("compute_B_lasso", {
dat <- lfmm_sampler(n = 100, p = 1000, K = 3,
outlier.prop = 0.1,
cs = c(0.8),
sigma = 0.2,
B.sd = 1.0,
U.sd = 1.0,
V.sd = 1.0)
B_lasso <- function(A, X, lambda) {
B_hat <- solve((crossprod(X,X) ), crossprod(X, A))
sign(B_hat) * ((abs(B_hat) - lambda) %>% purrr::map_dbl(~ max(.x, 0)))
}
lambda <- 1.5e0
X <- cbind(matrix(1,100,1), dat$X)
X <- svd(X)$u ## to have othogonal value
B.pkg <- compute_B_lasso(A = dat$Y,
X = X,
lambda = lambda)
expect_equal(dim(B.pkg), c(1000, 2))
## hist(B[,2])
## mean(B[,2] == 0)
## with R
B.r <- B_lasso(dat$Y, X, lambda)
## comp
expect_lt(mean(abs(t(B.r) - B.pkg)), 1e-15)
expect_equal(mean(B.r != 0), mean(B.pkg != 0))
})
test_that("X svd with and ridge estimator", {
dat <- lfmm_sampler(n = 100, p = 1000, K = 3,
outlier.prop = 0.1,
cs = c(0.8),
sigma = 0.2,
B.sd = 1.0,
U.sd = 1.0,
V.sd = 1.0)
dat$X <- cbind(dat$X, rnorm(100))
lambda <- 1
n <- nrow(dat$X)
d <- ncol(dat$X)
res <- compute_eigen_svd(X = dat$X)
B <- compute_B_ridge(dat$Y, dat$X, lambda)
## compute B as on my cahier 6/07/2017
D2 <- diag(res$sigma, d, n)
diag(D2) <- res$sigma / (res$sigma ^ 2 + lambda)
B.hat <- t(res$R %*% D2 %*% t(res$Q) %*% dat$Y)
expect_lt(mean(abs(B - B.hat)), 1e-15)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.