Merge branch 'jelena-markovic-randomized_jelena'

Author: jonathan-taylor

selectiveInference/R/funs.randomized.R

@@ -33,7 +33,7 @@ randomizedLasso = function(X,
     if (is.null(noise_scale)) {
         noise_scale = 0.5 * sd(y) * sqrt(mean_diag)
     }
-
+    
     noise_type = match.arg(noise_type)
 
     if (noise_scale > 0) {
@@ -65,8 +65,8 @@ randomizedLasso = function(X,
     nactive = as.integer(0)
 
     result = solve_QP_wide(X,                  # design matrix
-    	                   lam / n,            # vector of Lagrange multipliers
-		           ridge_term / n,     # ridge_term 
+    	                     lam / n,            # vector of Lagrange multipliers
+		                       ridge_term / n,     # ridge_term 
                            max_iter, 
                            soln, 
                            linear_func, 
@@ -76,12 +76,12 @@ randomizedLasso = function(X,
                            nactive, 
                            kkt_tol, 
                            objective_tol, 
-			   parameter_tol,
+			                     parameter_tol,
                            p,
-		           objective_stop,     # objective_stop
-			   kkt_stop,           # kkt_stop
-			   parameter_stop)         # param_stop
-
+		                       objective_stop,     # objective_stop
+			                     kkt_stop,           # kkt_stop
+			                     parameter_stop)         # param_stop
+    
     sign_soln = sign(result$soln)
 
     unpenalized = lam == 0
@@ -96,7 +96,11 @@ randomizedLasso = function(X,
 
     observed_scalings = abs(result$soln)[active]
     observed_unpen = result$soln[unpenalized]
-    observed_subgrad = result$gradient[inactive]
+    observed_subgrad = -n*result$gradient[inactive]
+    
+    if (length(which(abs(observed_subgrad)>lam[1]))){
+      print("subgradient eq not satisfied")
+    }
 
     observed_opt_state = c(observed_unpen, observed_scalings, observed_subgrad)
 
@@ -111,14 +115,15 @@ randomizedLasso = function(X,
     coef_term = L_E
 
     signs_ = c(rep(1, sum(unpenalized)), sign_soln[active])
+    
+    coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active)))  # ridge term
+  
     if (length(signs_) == 1) {
-        coef_term = coef_term * signs_
+      coef_term = coef_term * signs_
     } else {
-        coef_term = coef_term %*% diag(signs_)  # scaligns are non-negative
+      coef_term = coef_term %*% diag(signs_)  # scaligns are non-negative
     }
-
-    coef_term[active,] = coef_term[active,] + ridge_term * diag(rep(1, sum(active)))  # ridge term
-
+    
     subgrad_term = matrix(0, p, sum(inactive)) # for subgrad
     for (i in 1:sum(inactive)) {
         subgrad_term[inactive_set[i], i] = 1
@@ -155,7 +160,8 @@ randomizedLasso = function(X,
     inactive_lam = lam[inactive_set]
     inactive_start = sum(unpenalized) + sum(active)
     active_start = sum(unpenalized)
-
+    
+    
     # XXX only for Gaussian so far
 
     log_optimization_density = function(opt_state) {
@@ -185,9 +191,11 @@ randomizedLasso = function(X,
                 optimization_transform = opt_transform,
                 internal_transform = internal_transform,
                 log_optimization_density = log_optimization_density,
-		observed_opt_state = observed_opt_state,
+		            observed_opt_state = observed_opt_state,
                 observed_raw = observed_raw,
-		noise_scale = noise_scale
+		            noise_scale = noise_scale,
+		            soln = result$soln,
+		            perturb = perturb_
                 ))
 
 }
@@ -239,7 +247,6 @@ importance_weight = function(noise_scale,
         A = apply(A, 2, function(x) {return(x + target_transform$offset_term + opt_transform$offset_term)})
         log_num = -apply(A^2, 2, sum) / noise_scale^2
     } else {
-
         log_num = log_density_gaussian_(noise_scale,
                                         target_transform$linear_term,
                                         as.matrix(target_sample),
@@ -262,6 +269,15 @@ importance_weight = function(noise_scale,
     W = W - max(W)
     return(exp(W))
 }
+
+get_mean_cov = function(noise_scale, linear_term, offset_term){
+    temp = solve(t(linear_term) %*% linear_term)
+    cov = noise_scale^2*temp
+    mean = temp %*% t(linear_term) %*% offset_term
+    return(list(mean=mean, cov=cov))
+}
+
+
 
 conditional_density = function(noise_scale, lasso_soln) {
 
@@ -306,7 +322,9 @@ conditional_density = function(noise_scale, lasso_soln) {
   lasso_soln$log_optimization_density = log_condl_optimization_density
   lasso_soln$observed_opt_state = observed_opt_state[1:nactive]
   lasso_soln$optimization_transform = opt_transform
-  return(lasso_soln)
+  reduced_opt_transform =list(linear_term = reduced_B, offset_term = reduced_beta_offset)
+  return(list(lasso_soln=lasso_soln, 
+              reduced_opt_transform = reduced_opt_transform))
 }
 
 randomizedLassoInf = function(X, 
@@ -317,8 +335,9 @@ randomizedLassoInf = function(X,
                               ridge_term=NULL, 
                               condition_subgrad=TRUE, 
                               level=0.9,
-			      nsample=10000,
-			      burnin=2000,
+                              sampler=c("norejection", "adaptMCMC"),
+                              nsample=10000,
+                              burnin=2000,
                               max_iter=100,        # how many iterations for each optimization problem
                               kkt_tol=1.e-4,       # tolerance for the KKT conditions
                               parameter_tol=1.e-8, # tolerance for relative convergence of parameter
@@ -345,22 +364,47 @@ randomizedLassoInf = function(X,
                                parameter_stop=parameter_stop)
 
   active_set = lasso_soln$active_set
-  if (length(active_set)==0){
+  nactive = length(active_set)
+
+  if (nactive==0){
     return (list(active_set=active_set, pvalues=c(), ci=c()))
   }
   inactive_set = lasso_soln$inactive_set
-  nactive = length(active_set)
-
+  
   noise_scale = lasso_soln$noise_scale # set to default value in randomizedLasso
-
- if (condition_subgrad==TRUE){
-   lasso_soln=conditional_density(noise_scale, lasso_soln)
- } 
+  
+  constraints = matrix(0,nactive,2)
+  constraints[,2] = Inf
+  if (condition_subgrad==TRUE){
+   condl_lasso=conditional_density(noise_scale, lasso_soln)
+   lasso_soln = condl_lasso$lasso_soln
+   cur_opt_transform = condl_lasso$reduced_opt_transform
+   } else{
+   if (nactive<p){
+     subgrad_constraints = matrix(-lam, p-nactive, 2)
+     subgrad_constraints[,2]=lam
+     constraints = rbind(constraints, subgrad_constraints)
+   }
+    cur_opt_transform = list(linear_term = lasso_soln$optimization_transform$linear_term,
+                             offset_term = lasso_soln$optimization_transform$offset_term+lasso_soln$observed_raw)
+  }
 
   ndim = length(lasso_soln$observed_opt_state)
-
-  S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), ndim), nsample=nsample)
-  opt_samples = as.matrix(S$samples[(burnin+1):nsample,,drop=FALSE])
+  
+  sampler = match.arg(sampler)
+
+  if (sampler == "adaptMCMC"){
+    S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), ndim), nsample=nsample)
+    opt_samples = as.matrix(S$samples[(burnin+1):nsample,,drop=FALSE])
+  } else if (sampler == "norejection") {
+    opt_samples = gaussian_sampler(noise_scale, 
+                                 lasso_soln$observed_opt_state, 
+                                 cur_opt_transform$linear_term,
+                                 cur_opt_transform$offset_term,
+                                 constraints,
+                                 nsamples=nsample)
+    opt_sample = opt_samples[(burnin+1):nsample,]
+  }
 
   X_E = X[, active_set]
   X_minusE = X[, inactive_set]
@@ -404,7 +448,6 @@ randomizedLassoInf = function(X,
         cur_linear = reduced_target_opt_linear[,2:ncol(reduced_target_opt_linear)]
         cur_offset = temp %*% opt_transform$offset_term
         cur_transform = list(linear_term = as.matrix(cur_linear), offset_term = cur_offset)
-
         raw = target_transform$linear_term * observed_target[i] + target_transform$offset_term
     } else {
         cur_transform = opt_transform
@@ -450,3 +493,6 @@ randomizedLassoInf = function(X,
   }
   return(list(active_set=active_set, pvalues=pvalues, ci=ci))
 }
+
+
+

selectiveInference/R/sampler.R

@@ -0,0 +1,98 @@
+
+log_concave_sampler = function(negative_log_density, 
+                               grad_negative_log_density, 
+                               constraints,
+                               observed,
+                               nsamples){
+  #print(constraints)
+  constraints = as.matrix(constraints)
+  dim = nrow(constraints)
+  
+  get_poisson_process = function(state){
+    pos = as.matrix(state$pos)
+    velocity = as.matrix(state$velocity)
+    neg_velocity = velocity<0
+    pos_velocity = velocity>0
+    tau_min = 0
+    tau_max = 10
+    if (sum(neg_velocity)>0){
+      R = (-constraints[neg_velocity,1]+pos[neg_velocity])/(-velocity[neg_velocity])
+      tau_max = min(tau_max, min(R))
+      L = (-constraints[neg_velocity,2]+pos[neg_velocity])/(-velocity[neg_velocity])
+      tau_min = max(tau_min, max(L))
+    }
+    if (sum(pos_velocity)>0){
+      R = (constraints[pos_velocity,2]-pos[pos_velocity])/velocity[pos_velocity]
+      tau_max = min(tau_max, min(R))
+      L = (constraints[pos_velocity,1]-pos[pos_velocity])/velocity[pos_velocity]
+      tau_min = max(tau_min, max(L))
+    }
+    
+    f=function(t){as.numeric(t(velocity) %*% grad_negative_log_density(pos+velocity*t))}
+    tau_star = tau_max
+    if (f(tau_min)*f(tau_max)<0){
+      tau_star = uniroot(f, c(tau_min, tau_max))$root
+    } else{
+      if (negative_log_density(pos+velocity*tau_min)<negative_log_density(pos+velocity*tau_max)){
+        tau_star = tau_min
+      }
+    }
+    
+    tau_min = max(tau_min, tau_star)
+    
+    RHS = negative_log_density(pos+velocity*tau_star)+rexp(1)
+    g = function(t){negative_log_density(pos+velocity*t)-RHS}
+    if (g(tau_min)*g(tau_max)<0){
+      tau = uniroot(g, c(tau_min, tau_max))$root
+    } else{
+      tau = tau_max
+    }
+    return (tau)
+  }
+  
+  update_velocity = function(){
+    Z=rnorm(dim)
+    return(Z/sqrt(t(Z)%*%Z))
+  }
+  
+  compute_next = function(state){
+    bounce_time = get_poisson_process(state)/2
+    #print(paste("bounce time", bounce_time))
+    next_pos = state$pos+state$velocity*bounce_time
+    next_velocity=update_velocity()
+    return(list(pos=next_pos, velocity=next_velocity))
+  }
+  
+  state = list(pos=observed, velocity = update_velocity())
+  samples = matrix(0, nrow = nsamples, ncol = dim)
+  for (i in 1:nsamples){
+    #print(paste("pos", toString(state$pos)))
+    #print(paste("velocity", toString(state$velocity)))
+    samples[i,]=state$pos
+    state = compute_next(state)
+  }
+  return (samples)
+}
+
+gaussian_sampler = function(noise_scale, 
+                            observed, 
+                            linear_term, 
+                            offset_term, 
+                            constraints,
+                            nsamples){
+  
+  negative_log_density = function(x) {
+    recon = linear_term %*% x+offset_term
+    return(as.numeric(t(recon)%*%recon/(2*noise_scale^2)))
+  }
+  grad_negative_log_density=function(x){
+    recon = linear_term %*% x+offset_term
+    return(t(linear_term)%*% recon/(noise_scale^2))
+  }
+  
+  return(log_concave_sampler(negative_log_density, 
+                             grad_negative_log_density,
+                             constraints,
+                             observed,
+                             nsamples))
+}

selectiveInference/man/randomizedLassoInf.Rd

@@ -19,6 +19,7 @@ randomizedLassoInf(X,
                    ridge_term=NULL, 
                    condition_subgrad=TRUE, 
                    level=0.9,
+                   sampler=c("norejection", "adaptMCMC"),
                    nsample=10000,
                    burnin=2000,
                    max_iter=100,       
@@ -71,6 +72,10 @@ Default is TRUE.
 \item{level}{
 Level for confidence intervals.
 }
+\item{sampler}{
+Which sampler to use -- default is a no-rejection sampler. Otherwise
+use MCMC from the adaptMCMC package.
+}
 \item{nsample}{
 Number of samples of optimization variables to sample.
 }