Merge branch 'master' of https://github.com/selective-inference/R-software

Author: ryantibs

forLater/josh/selectiveInference/R/cv.R

@@ -1,11 +1,18 @@
 # ------------------------------------------------
 # Cross-validation, preliminary
 
-cvMakeFolds <- function(x, nfolds = 10) {
-    #inds <- sample(1:nrow(x), replace=FALSE)
-    inds <- 1:nrow(x)
+cvMakeFolds <- function(x, nfolds = 5) {
+    inds <- sample(1:nrow(x), replace=FALSE)
+    #inds <- 1:nrow(x)
     foldsize <- floor(nrow(x)/nfolds)
-    lapply(1:nfolds, function(f) return(inds[1:foldsize+(f-1)*foldsize]))
+    folds <- lapply(1:nfolds, function(f) return(inds[1:foldsize+(f-1)*foldsize]))
+    if (nfolds*foldsize < nrow(x)) {
+      # remainder observations added to first several folds
+      for (i in 1:(nrow(x) - nfolds*foldsize)) {
+        folds[[i]] <- c(folds[[i]], inds[nfolds*foldsize + i])
+      }
+    }
+    return(folds)
 }
 
 ############################################
@@ -87,6 +94,7 @@ cvfs <- function(x, y, index = 1:ncol(x), maxsteps, sigma = NULL, intercept = TR
         fold <- folds[[f]]
         fit <- groupfs(X[-fold,], Y[-fold], index=index, maxsteps=maxsteps, sigma=sigma, intercept=FALSE, center=FALSE, normalize=FALSE)
         fit$fold <- fold
+        # Why is this commented out?
         ## projections[[f]] <- lapply(fit$projections, function(step.projs) {
         ##     lapply(step.projs, function(proj) {
         ##         # Reduce from n by n matrix to svdu_thresh
@@ -122,3 +130,48 @@ cvfs <- function(x, y, index = 1:ncol(x), maxsteps, sigma = NULL, intercept = TR
     invisible(fit)
 }
 
+
+cvlar <- function(x, y) { # other args
+    folds <- cvMakeFolds(x)
+    models <- lapply(folds, function(fold) {
+        x.train <- X
+        y.train <- Y
+        x.train[fold,] <- 0
+        y.train[fold] <- 0
+        x.test <- X[fold,]
+        y.test <- Y[fold]
+          larpath.train <- lar(x.train, y.train, maxsteps = maxsteps, intercept = F, normalize = F)
+        return(lff)
+    })
+
+    active.sets <- lapply(models, function(model) model$action)
+    lambdas <- lapply(models, function(model) model$lambda)
+    lmin <- min(unlist(lambdas))
+
+# Interpolate lambda grid or parametrize by steps?
+# interpolation probably requires re-writing cvRSSquads for
+# penalized fits in order to make sense
+
+# do steps for now just to have something that works?
+
+    RSSquads <- list()
+    for (s in 1:maxsteps) {
+        initial.active <- lapply(active.sets, function(a) a[1:s])
+        RSSquads[[s]] <- cvRSSquad(X, folds, initial.active)
+    }
+
+    RSSs <- lapply(RSSquads, function(Q) t(Y) %*% Q %*% Y)
+    sstar <- which.min(RSSs)
+    quadstar <- RSSquads[sstar][[1]]
+
+    RSSquads <- lapply(RSSquads, function(quad) quad - quadstar)
+    RSSquads[[sstar]] <- NULL # remove the all zeroes case
+
+    fit <- lar(X, Y, maxsteps=sstar, intercept = F, normalize = F)
+
+# Very tall Gamma encoding all cv-model paths
+    Gamma <- do.call(rbind, lapply(models, function(model) return(model$Gamma)))
+
+# more to do here    
+}
+

forLater/josh/sim.cvlar.R

@@ -0,0 +1,58 @@
+# Choices
+
+# RSS: least-squares or penalized beta?
+# depends on final model. Go with least-squares for now
+
+# fixed vs lar? (lar, apparently)
+# fixed probably slower, but advantage of same lambda grid?
+# is same lambda grid necessary? -- doesn't lar algorithm give all possible models anyway?
+# i.e. for non-knot lambda just find where it is in lar path, take corresponding model
+
+# groups? later
+
+# TODO
+
+# copy larInf or groupfsInf?
+# larInf: add CV quadratic constraints* & break/fix p-value computation
+# -------- *but can we even use the ydecomp we use for quadratic?
+# groupfsInf: some ugly rewriting, no cumprojs etc, but straightforward
+# -------- downside: need to implement larInf basically
+
+# larInf
+# [ ] is.null(sigma) don't estimate it
+
+# plan:
+# expand Gamma for [-fold] indices?
+# stack all the Gammas? or iterate through them?
+# work backward from poly.pval <- larInf
+
+
+# big picture / long term
+# what OOP kind of design would lend itself to easily implementing more cv things?
+
+# Gamma: something x n
+# Gamma %*% y >= 0
+
+# pass 0-padded x[-fold] and y[-fold] to lar?
+
+library(selectiveInference)
+setwd("/Users/joftius/Dropbox/work/R-software/forLater/josh")
+source("selectiveInference/R/cv.R")
+
+set.seed(1)
+n <- 100
+p <- 50
+maxsteps <- 10
+sparsity <- 3
+snr <- 2
+rho <- 0.1
+nfolds <- 5
+
+x <- matrix(rnorm(n*p), nrow=n)
+y <- rnorm(n)
+beta <- rep(0, p)
+beta[1:sparsity] <- 2* sqrt(2*log(p)/n) * sample(c(-1,1), sparsity, replace=T)
+y <- y + x %*% beta
+my <- mean(y)
+y <- y - my
+

forLater/josh/sim.datasplit.R

@@ -4,30 +4,31 @@ source("../../selectiveInference/R/funs.groupfs.R")
 source("../../selectiveInference/R/funs.quadratic.R")
 source("../../selectiveInference/R/funs.common.R")
 
-set.seed(1)
-niters <- 400
+set.seed(19)
+niters <- 500
 known <- FALSE
-n <- 100
+n <- 50
 p <- 100
-maxsteps <- 10
+maxsteps <- 8
 sparsity <- 5
-snr <- 1
+snr <- 2
 rho <- 0.1
-ratio <- 0.7
-ratio2 <- 0.85
+ratio <- 0.6
+ratio2 <- 0.8
 train <- 1:(ratio*n)
 test <- setdiff(1:n, train)
 train2 <- 1:(ratio2*n)
-test <- setdiff(1:n, train2)
+test2 <- setdiff(1:n, train2)
 index <- 1:p
 
-instance <- function(n, p, sparsity, snr, maxsteps, rho) {
-
-    x <- matrix(rnorm(n*p), nrow=n)
+x <- matrix(rnorm(n*p), nrow=n)
     if (rho != 0) {
         z <- matrix(rep(t(rnorm(n)), p), nrow = n)
         x <- sqrt(1-rho)*x + sqrt(rho)*z
     }
+
+instance <- function(n, p, sparsity, snr, maxsteps, rho) {
+
     y <- rnorm(n)
 
     if (sparsity > 0) {
@@ -47,23 +48,24 @@ instance <- function(n, p, sparsity, snr, maxsteps, rho) {
     xte2 <- x[test2, ]
 
     if (known) {
-        trfit <- groupfs(xtr, ytr, index, maxsteps=maxsteps, sigma=1, aicstop=1, k = 2*log(p))
-        fit <- groupfs(xtr2, ytr2, index, maxsteps=maxsteps, sigma=1, aicstop=1, k = 2*log(p))
+        trfit <- groupfs(xtr, ytr, index, maxsteps=maxsteps, sigma=1, aicstop=1, k = log(length(train)))
+        fit <- groupfs(xtr2, ytr2, index, maxsteps=maxsteps, sigma=1, aicstop=1, k = log(length(train2)))
     } else {
-        trfit <- groupfs(xtr, ytr, index, maxsteps=maxsteps, aicstop=1, k = log(length(train)))
-        fit <- groupfs(xtr2, ytr2, index, maxsteps=maxsteps, aicstop=1, k = log(length(train2)))
+        trfit <- groupfs(xtr, ytr, index, maxsteps=maxsteps, aicstop=1, k = 2*log(p))
+        fit <- groupfs(xtr2, ytr2, index, maxsteps=maxsteps, aicstop=1, k = 2*log(p))
     }
 
     trcols <- which(1:p %in% trfit$action)
     tr2cols <- which(1:p %in% fit$action)
     tepv <- summary(lm(yte~xte[, trcols]-1))$coefficients[,4]
     tepv2 <- summary(lm(yte2~xte2[, tr2cols]-1))$coefficients[,4]
     names(tepv) <- as.character(sort(trfit$action))
-    names(tepv2) <- as.character(sort(trfit$action))
+    names(tepv2) <- as.character(sort(fit$action))
     pv <- groupfsInf(fit)
     trpv <- groupfsInf(trfit)
     return(list(vars = fit$action, pvals = pv$pv,
                 splitvars = sort(trfit$action), splitpvals = tepv,
+                splitvars2 = sort(fit$action), splitpvals2 = tepv2,                
                 trpvals = trpv$pv))
 }
 
@@ -75,9 +77,12 @@ vars <- do.call(c, list(output[1,]))
 pvals <- do.call(c, list(output[2,]))
 splitvars <- do.call(c, list(output[3,]))
 splitpvals <- do.call(c, list(output[4,]))
-trpvals <- do.call(c, list(output[5,]))
+splitvars2 <- do.call(c, list(output[5,]))
+splitpvals2 <- do.call(c, list(output[6,]))
+trpvals <- do.call(c, list(output[7,]))
 
-save(vars, pvals, splitvars, splitpvals, trpvals,
+save(vars, pvals, splitvars, splitpvals,
+     splitvars2, splitpvals2, trpvals,
      file = paste0("results/datasplit",
          "_", ifelse(known, "TC", "TF"),
          "_n", n,

forLater/josh/sim.selectedmodel.R

@@ -0,0 +1,60 @@
+library(selectiveInference)
+library(intervals)
+setwd("~/Dropbox/work/R-software/forLater/josh")
+source("selectiveInference/R/cv.R")
+source("../../selectiveInference/R/funs.groupfs.R")
+source("../../selectiveInference/R/funs.quadratic.R")
+source("../../selectiveInference/R/funs.common.R")
+source("../../selectiveInference/R/funs.fs.R")
+source("../../selectiveInference/R/funs.lar.R")
+source("../../selectiveInference/R/funs.inf.R")
+library(MASS)
+pinv = ginv
+
+set.seed(19)
+niters <- 500
+known <- TRUE
+n <- 50
+p <- 100
+maxsteps <- 8
+sparsity <- 5
+snr <- 2
+index <- 1:p
+
+x <- matrix(rnorm(n*p), nrow=n)
+
+instance <- function(n, p, sparsity, snr, maxsteps) {
+    y <- rnorm(n)
+    if (sparsity > 0) {
+      beta <- rep(0, p)
+      beta[1:sparsity] <- snr * sample(c(-1,1), sparsity, replace=T)
+      y <- y + x %*% beta
+    }
+    y <- y - mean(y)
+    fit <- groupfs(x, y, index, maxsteps=maxsteps, sigma=1, intercept=F, center=F, normalize=F)
+    fitfs <- fs(x, y, maxsteps=maxsteps, intercept=F, normalize=F)
+    if (any(fit$action != fitfs$action)) stop("Model paths did not agree")
+    pvfs <- fsInf(fitfs, sigma=1, k = maxsteps, type = "all")
+    pv <- groupfsInf(fit)
+    return(list(vars = fit$action, pvals = pv$pv, selpvals = pvfs$pv))
+}
+
+time <- system.time({
+          output <- replicate(niters, instance(n, p, sparsity, snr, maxsteps))
+})
+
+vars <- do.call(c, list(output[1,]))
+pvals <- do.call(c, list(output[2,]))
+selpvals <- do.call(c, list(output[3,]))
+
+save(vars, pvals, selpvals,
+     file = paste0("results/selected",
+         "_", ifelse(known, "TC", "TF"),
+         "_n", n,
+         "_p", p,
+         "_sparsity", sparsity,
+         "_snr", as.character(snr),
+         ".RData"))
+
+print(time)
+