Merge branch 'jelena-markovic-randomized_jelena'

Author: jonathan-taylor

selectiveInference/R/funs.randomized.R

@@ -3,11 +3,11 @@
 #
 # min 1/2 || y - \beta_0 - X \beta ||_2^2 + \lambda || \beta ||_1 - \omega^T\beta + \frac{\epsilon}{2} \|\beta\|^2_2
 
-randomizedLASSO = function(X, 
+randomizedLasso = function(X, 
                            y, 
                            lam, 
-                           noise_scale, 
-                           ridge_term, 
+                           noise_scale=NULL, 
+                           ridge_term=NULL, 
                            noise_type=c('gaussian', 'laplace'),
                            max_iter=100,        # how many iterations for each optimization problem
                            kkt_tol=1.e-4,       # tolerance for the KKT conditions
@@ -20,6 +20,21 @@ randomizedLASSO = function(X,
 
     n = nrow(X); p = ncol(X)
 
+    mean_diag = mean(apply(X^2, 2, sum))
+
+    # default ridge term
+
+    if (is.null(ridge_term)) {
+        ridge_term = sqrt(mean_diag) * sd(y) / sqrt(n)
+    }
+
+    # default noise level
+
+    if (is.null(noise_scale)) {
+        noise_scale = 0.5 * sd(y) * sqrt(mean_diag)
+    }
+
+    print(c(noise_scale, ridge_term))
     noise_type = match.arg(noise_type)
 
     if (noise_scale > 0) {
@@ -246,10 +261,10 @@ conditional_density = function(noise_scale, lasso_soln) {
     if  (sum(opt_state < 0) > 0) {
       return(-Inf)
     }
-    D = selectiveInference:::log_density_gaussian_conditional_(noise_scale,
-                                                               reduced_B,
-                                                               as.matrix(opt_state),
-                                                               reduced_beta_offset)
+    D = log_density_gaussian_conditional_(noise_scale,
+                                          reduced_B,
+                                          as.matrix(opt_state),
+                                          reduced_beta_offset)
     return(D)
   }
   lasso_soln$log_optimization_density = log_condl_optimization_density
@@ -258,23 +273,42 @@ conditional_density = function(noise_scale, lasso_soln) {
   return(lasso_soln)
 }
 
-randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term){
+randomizedLassoInf = function(X, 
+                              y, 
+                              lam, 
+                              sigma=NULL, 
+                              noise_scale=NULL, 
+                              ridge_term=NULL, 
+                              condition_subgrad=TRUE, 
+                              level=0.9) {
 
   n = nrow(X)
   p = ncol(X)
-  lasso_soln = selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
+  lasso_soln = randomizedLasso(X, y, lam, noise_scale, ridge_term)
   active_set = lasso_soln$active_set
   inactive_set = lasso_soln$inactive_set
   nactive = length(active_set)
-
+  
+  if (condition_subgrad==TRUE){
+    lasso_soln=conditional_density(noise_scale,lasso_soln)
+  } 
+    
   dim = length(lasso_soln$observed_opt_state)
   print(paste("chain dim", dim))
-  S = selectiveInference:::sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=10000)
+  S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=10000)
   opt_samples = S$samples[2001:10000,]
   print(paste("dim opt samples", toString(dim(opt_samples))))
 
   X_E = X[, active_set]
   X_minusE = X[, inactive_set]
+
+  # if no sigma given, use OLS estimate
+
+  if (is.null(sigma)) {
+        lm_y = lm(y ~ X[,E] - 1)
+        sigma = sum(resid(lm_y)^2 / lm_y$df.resid)
+  }        
+  print(c(sigma, 'sigma'))
   target_cov = solve(t(X_E) %*% X_E)*sigma^2
   cov_target_internal = rbind(target_cov, matrix(0, nrow=p-nactive, ncol=nactive))
   observed_target = solve(t(X_E) %*% X_E) %*% t(X_E) %*% y
@@ -283,37 +317,46 @@ randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term){
   opt_transform = lasso_soln$optimization_transform
   observed_raw = lasso_soln$observed_raw
 
-  pvalus = rep(0, nactive)
+  pvalues = rep(0, nactive)
   ci = matrix(0, nactive, 2)
   for (i in 1:nactive){
-    target_transform = selectiveInference:::linear_decomposition(observed_target[i], 
-                                                  observed_internal, 
-                                                  target_cov[i,i], 
-                                                  cov_target_internal[,i],
-                                                  internal_transform)
+    target_transform = linear_decomposition(observed_target[i], 
+                                            observed_internal, 
+                                            target_cov[i,i], 
+                                            cov_target_internal[,i],
+                                            internal_transform)
     target_sample = rnorm(nrow(opt_samples)) * sqrt(target_cov[i,i])
 
     pivot = function(candidate){
-      weights = selectiveInference:::importance_weight(noise_scale,
-                                                     t(as.matrix(target_sample)) + candidate,
-                                                     t(opt_samples),
-                                                     opt_transform,
-                                                     target_transform,
-                                                     observed_raw)
-      return(mean((target_sample<observed_target[i])*weights)/mean(weights))
+      weights = importance_weight(noise_scale,
+                                  t(as.matrix(target_sample)) + candidate,
+                                  t(opt_samples),
+                                  opt_transform,
+                                  target_transform,
+                                  observed_raw)
+      return(mean((target_sample+candidate<observed_target[i])*weights)/mean(weights))
     }
-    level = 0.9
     rootU = function(candidate){
       return (pivot(observed_target[i]+candidate)-(1-level)/2)
     }
     rootL = function(candidate){
       return (pivot(observed_target[i]+candidate)-(1+level)/2)
     }
     pvalues[i] = pivot(0)
-    line_min = -10*sd(target_sample)
-    line_max = 10*sd(target_sample)
-    ci[i,1] = uniroot(rootU, c(line_min, line_max))$root+observed_target[i]
-    ci[i,2] = uniroot(rootL, c(line_min, line_max))$root+observed_target[i]
+    line_min = -20*sd(target_sample)
+    line_max = 20*sd(target_sample)
+    if (rootU(line_min)*rootU(line_max)<0){
+      ci[i,2] = uniroot(rootU, c(line_min, line_max))$root+observed_target[i]
+    } else{
+      print("non inv u")
+      ci[i,2]=line_max
+    }
+    if (rootL(line_min)*rootL(line_max)<0){
+      ci[i,1] = uniroot(rootL, c(line_min, line_max))$root+observed_target[i]
+    } else{
+      print("non inv u")
+      ci[i,1] = line_min
+    }
   }
-  return(list(pvalues=pvalues, ci=ci))
+  return(list(active_set=active_set, pvalues=pvalues, ci=ci))
 }

tests/randomized/test_instances.R

@@ -25,28 +25,32 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
 }
 
 
-collect_results = function(n,p,s, nsim=10){
+collect_results = function(n,p,s, nsim=100, level=0.9){
   rho=0.3
   lam=1.
   sigma=1
-  sample_pivots = c()
+  sample_pvalues = c()
+  sample_coverage = c()
   for (i in 1:nsim){
     data = gaussian_instance(n=n,p=p,s=s, rho=rho, sigma=sigma)
     X=data$X
     y=data$y
-    ridge_term=sd(y)/sqrt(n)
-    noise_scale = sd(y)/2
-    #X = matrix(rnorm(n * p), n, p)
-    #y = rnorm(n)
-    #lam = 20 / sqrt(n)
-    #noise_scale = 0.01 * sqrt(n)
-    #ridge_term = .1 / sqrt(n)
-    result = selectiveInference:::randomized_inference(X,y,sigma,lam,noise_scale,ridge_term)
-    sample_pivots = c(sample_pivots, result$pivots)
+    beta=data$beta
+    result = selectiveInference:::randomizedLassoInf(X, y, sigma, lam, level=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
+      }
+      print(paste("ci", toString(result$ci[i,])))
+    }
+    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_pivots), xlim=c(0,1),  main="Empirical CDF of null p-values", xlab="p-values", ylab="ecdf")
+  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()
 }

tests/randomized/test_randomized.R

@@ -7,7 +7,7 @@ smoke_test = function() {
     lam = 20 / sqrt(n)
     noise_scale = 0.01 * sqrt(n)
     ridge_term = .1 / sqrt(n)
-    selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
+    selectiveInference:::randomizedLasso(X, y, lam, noise_scale, ridge_term)
 }
 
 A = smoke_test()
@@ -20,7 +20,7 @@ sampler_test = function() {
     lam = 20 / sqrt(n)
     noise_scale = 0.01 * sqrt(n)
     ridge_term = .1 / sqrt(n)
-    obj = selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
+    obj = selectiveInference:::randomizedLasso(X, y, lam, noise_scale, ridge_term)
     S = selectiveInference:::sample_opt_variables(obj, jump_scale=rep(1/sqrt(n), p), nsample=10000)
     return(S$samples[2001:10000,])
 }