dealing with zero length active set

Author: Jelena Markovic

selectiveInference/R/funs.randomized.R

@@ -264,6 +264,9 @@ randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term, level
   p = ncol(X)
   lasso_soln = selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
   active_set = lasso_soln$active_set
+  if (length(active_set)==0){
+    return (list(active_set=active_set, pvalues=c(), ci=c()))
+  }
   inactive_set = lasso_soln$inactive_set
   nactive = length(active_set)
 
@@ -293,8 +296,9 @@ randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term, level
                                                   target_cov[i,i], 
                                                   cov_target_internal[,i],
                                                   internal_transform)
+  
     target_sample = rnorm(nrow(as.matrix(opt_samples))) * sqrt(target_cov[i,i])
-    
+    print(length(target_sample))
     pivot = function(candidate){
       weights = selectiveInference:::importance_weight(noise_scale,
                                                      t(as.matrix(target_sample)) + candidate,

tests/randomized/test_instances.R

@@ -39,22 +39,26 @@ collect_results = function(n,p,s, nsim=1, level=0.9){
     ridge_term=sd(y)/sqrt(n)
     noise_scale = sd(y)/2
     result = selectiveInference:::randomized_inference(X,y,sigma,lam,noise_scale,ridge_term, level)
-    true_beta = beta[result$active_set]
-    coverage = rep(0, nrow(result$ci))
-    for (i in 1:nrow(result$ci)){
-      if (result$ci[i,1]<true_beta[i] & result$ci[i,2]>true_beta[i]){
-        coverage[i]=1
+    if (length(result$active_set)>0){
+      true_beta = beta[result$active_set]
+      coverage = rep(0, nrow(result$ci))
+      for (i in 1:nrow(result$ci)){
+        if (result$ci[i,1]<true_beta[i] & result$ci[i,2]>true_beta[i]){
+          coverage[i]=1
+        }
+        print(paste("ci", toString(result$ci[i,])))
       }
-      print(paste("ci", toString(result$ci[i,])))
+      sample_pvalues = c(sample_pvalues, result$pvalues)
+      sample_coverage = c(sample_coverage, coverage)
     }
-    sample_pvalues = c(sample_pvalues, result$pvalues)
-    sample_coverage = c(sample_coverage, coverage)
   }
-  print(paste("coverage", mean(sample_coverage)))
-  jpeg('pivots.jpg')
-  plot(ecdf(sample_pvalues), xlim=c(0,1),  main="Empirical CDF of null p-values", xlab="p-values", ylab="ecdf")
-  abline(0, 1, lty=2)
-  dev.off()
+  if (length(sample_coverage)>0){
+    print(paste("coverage", mean(sample_coverage)))
+    jpeg('pivots.jpg')
+    plot(ecdf(sample_pvalues), xlim=c(0,1),  main="Empirical CDF of null p-values", xlab="p-values", ylab="ecdf")
+    abline(0, 1, lty=2)
+    dev.off()
+  }
 }
 
 set.seed(1)