|
| 1 | +library(testthat) |
| 2 | +library(biglasso) |
| 3 | +library(ncvreg) |
| 4 | + |
| 5 | +context("Testing linear regression:") |
| 6 | + |
| 7 | +test_that("Test against OLS:", { |
| 8 | + ## test against OLS |
| 9 | + set.seed(1234) |
| 10 | + n <- 100 |
| 11 | + p <- 10 |
| 12 | + eps <- 1e-10 |
| 13 | + tolerance <- 1e-6 |
| 14 | + X <- matrix(rnorm(n*p), n, p) |
| 15 | + b <- rnorm(p) |
| 16 | + y <- rnorm(n, X %*% b) |
| 17 | + fit.ols <- lm(y ~ X) |
| 18 | + beta <- fit.ols$coefficients |
| 19 | + |
| 20 | + X.bm <- as.big.matrix(X) |
| 21 | + fit <- biglasso(X.bm, y, screen = 'None', lambda = 0, eps = eps) |
| 22 | + fit.edpp <- biglasso(X.bm, y, screen = 'SEDPP', lambda = 0, eps = eps) |
| 23 | + fit.edpp.no.active <- biglasso(X.bm, y, screen = 'SEDPP-No-Active', lambda = 0, eps = eps) |
| 24 | + fit.ssr <- biglasso(X.bm, y, screen = 'SSR', eps = eps, lambda = 0) |
| 25 | + fit.ssr.dome <- biglasso(X.bm, y, screen = 'SSR-Dome', eps = eps, lambda = 0) |
| 26 | + fit.ssr.edpp <- biglasso(X.bm, y, screen = 'SSR-BEDPP', eps = eps, lambda = 0) |
| 27 | + |
| 28 | + expect_equal(as.numeric(beta), as.numeric(fit$beta), tolerance = tolerance) |
| 29 | + expect_equal(as.numeric(beta), as.numeric(fit.edpp$beta), tolerance = tolerance) |
| 30 | + expect_equal(as.numeric(beta), as.numeric(fit.edpp.no.active$beta), tolerance = tolerance) |
| 31 | + expect_equal(as.numeric(beta), as.numeric(fit.ssr$beta), tolerance = tolerance) |
| 32 | + expect_equal(as.numeric(beta), as.numeric(fit.ssr.dome$beta), tolerance = tolerance) |
| 33 | + expect_equal(as.numeric(beta), as.numeric(fit.ssr.edpp$beta), tolerance = tolerance) |
| 34 | + expect_equal(as.numeric(beta), as.numeric(fit.edpp.no.active$beta), tolerance = tolerance) |
| 35 | + |
| 36 | +}) |
| 37 | + |
| 38 | +set.seed(1234) |
| 39 | +n <- 100 |
| 40 | +p <- 200 |
| 41 | +X <- matrix(rnorm(n*p), n, p) |
| 42 | +b <- c(rnorm(50), rep(0, p-50)) |
| 43 | +y <- rnorm(n, X %*% b) |
| 44 | +eps <- 1e-8 |
| 45 | +tolerance <- 1e-3 |
| 46 | +lambda.min <- 0.05 |
| 47 | + |
| 48 | +fit.ncv <- ncvreg(X, y, penalty = 'lasso', eps = sqrt(eps), lambda.min = lambda.min) |
| 49 | +cvfit.ncv <- cv.ncvreg(X, y, penalty = 'lasso', eps = sqrt(eps), |
| 50 | + lambda.min = lambda.min, seed = 1234, nfolds = 5) |
| 51 | + |
| 52 | +X.bm <- as.big.matrix(X) |
| 53 | +# fit <- biglasso(X.bm, y, screen = 'None', eps = eps) |
| 54 | +fit.edpp <- biglasso(X.bm, y, screen = 'SEDPP', eps = eps) |
| 55 | +# fit.edpp.no.active <- biglasso(X.bm, y, screen = 'SEDPP-No-Active', eps = eps) |
| 56 | +fit.ssr <- biglasso(X.bm, y, screen = 'SSR', eps = eps) |
| 57 | +fit.ssr.dome <- biglasso(X.bm, y, screen = 'SSR-Dome', eps = eps) |
| 58 | +fit.ssr.edpp <- biglasso(X.bm, y, screen = 'SSR-BEDPP', eps = eps) |
| 59 | + |
| 60 | +# cvfit <- cv.biglasso(X.bm, y, screen = 'None', eps = eps, |
| 61 | +# ncores = 1, nfolds = 5, seed = 1234) |
| 62 | +cvfit.edpp <- cv.biglasso(X.bm, y, screen = 'SEDPP', eps = eps, |
| 63 | + ncores = 1, nfolds = 5, seed = 1234) |
| 64 | +# cvfit.edpp.no.active <- cv.biglasso(X.bm, y, screen = 'SEDPP-No-Active', eps = eps, |
| 65 | +# ncores = 1, nfolds = 5, seed = 1234) |
| 66 | +cvfit.ssr <- cv.biglasso(X.bm, y, screen = 'SSR', eps = eps, |
| 67 | + ncores = 1, nfolds = 5, seed = 1234) |
| 68 | +cvfit.ssr.dome <- cv.biglasso(X.bm, y, screen = 'SSR-Dome', eps = eps, |
| 69 | + ncores = 1, nfolds = 5, seed = 1234) |
| 70 | +cvfit.ssr.edpp <- cv.biglasso(X.bm, y, screen = 'SSR-BEDPP', eps = eps, |
| 71 | + ncores = 1, nfolds = 5, seed = 1234) |
| 72 | + |
| 73 | +## parallel computing |
| 74 | +fit.edpp2 <- biglasso(X.bm, y, screen = 'SEDPP', eps = eps, ncores = 2) |
| 75 | +# fit.edpp.no.active2 <- biglasso(X.bm, y, screen = 'SEDPP-No-Active', eps = eps, ncores = 2) |
| 76 | +fit.ssr2 <- biglasso(X.bm, y, screen = 'SSR', eps = eps, ncores = 2) |
| 77 | +fit.ssr.dome2 <- biglasso(X.bm, y, screen = 'SSR-Dome', eps = eps, ncores = 2) |
| 78 | +fit.ssr.edpp2 <- biglasso(X.bm, y, screen = 'SSR-BEDPP', eps = eps, ncores = 2) |
| 79 | + |
| 80 | +test_that("Test against ncvreg for entire path:", { |
| 81 | + # expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit$beta), tolerance = tolerance) |
| 82 | + expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit.edpp$beta), tolerance = tolerance) |
| 83 | + # expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit.edpp.no.active$beta), tolerance = tolerance) |
| 84 | + expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit.ssr$beta), tolerance = tolerance) |
| 85 | + expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit.ssr.dome$beta), tolerance = tolerance) |
| 86 | + expect_equal(as.numeric(fit.ncv$beta), as.numeric(fit.ssr.edpp$beta), tolerance = tolerance) |
| 87 | +}) |
| 88 | + |
| 89 | +test_that("Test parallel computing: ",{ |
| 90 | + expect_identical(fit.edpp, fit.edpp2) |
| 91 | + # expect_identical(fit.edpp.no.active, fit.edpp.no.active2) |
| 92 | + expect_identical(fit.ssr, fit.ssr2) |
| 93 | + expect_identical(fit.ssr.dome, fit.ssr.dome2) |
| 94 | + expect_identical(fit.ssr.edpp, fit.ssr.edpp2) |
| 95 | +}) |
| 96 | + |
| 97 | +test_that("Test cross validation: ",{ |
| 98 | + # expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit$cve), tolerance = tolerance) |
| 99 | + expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit.edpp$cve), tolerance = tolerance) |
| 100 | + # expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit.edpp.no.active$cve), tolerance = tolerance) |
| 101 | + expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit.ssr$cve), tolerance = tolerance) |
| 102 | + expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit.ssr.dome$cve), tolerance = tolerance) |
| 103 | + expect_equal(as.numeric(cvfit.ncv$cve), as.numeric(cvfit.ssr.edpp$cve), tolerance = tolerance) |
| 104 | + |
| 105 | + expect_equal(as.numeric(cvfit.ncv$cvse), as.numeric(cvfit.edpp$cvse), tolerance = tolerance) |
| 106 | + expect_equal(as.numeric(cvfit.ncv$cvse), as.numeric(cvfit.ssr$cvse), tolerance = tolerance) |
| 107 | + expect_equal(as.numeric(cvfit.ncv$cvse), as.numeric(cvfit.ssr.dome$cvse), tolerance = tolerance) |
| 108 | + expect_equal(as.numeric(cvfit.ncv$cvse), as.numeric(cvfit.ssr.edpp$cvse), tolerance = tolerance) |
| 109 | + |
| 110 | + expect_equal(as.numeric(cvfit.ncv$lambda.min), as.numeric(cvfit.edpp$lambda.min), tolerance = tolerance) |
| 111 | + expect_equal(as.numeric(cvfit.ncv$lambda.min), as.numeric(cvfit.ssr$lambda.min), tolerance = tolerance) |
| 112 | + expect_equal(as.numeric(cvfit.ncv$lambda.min), as.numeric(cvfit.ssr.dome$lambda.min), tolerance = tolerance) |
| 113 | + expect_equal(as.numeric(cvfit.ncv$lambda.min), as.numeric(cvfit.ssr.edpp$lambda.min), tolerance = tolerance) |
| 114 | + |
| 115 | +}) |
0 commit comments