coverages

Author: Jelena Markovic

selectiveInference/R/RcppExports.R

@@ -0,0 +1,35 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+solve_QP <- function(Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, parameter_tol, max_active, kkt_stop, objective_stop, param_stop) {
+    .Call('selectiveInference_solve_QP', PACKAGE = 'selectiveInference', Sigma, bound, maxiter, theta, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol, parameter_tol, max_active, kkt_stop, objective_stop, param_stop)
+}
+
+solve_QP_wide <- function(X, bound, ridge_term, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, parameter_tol, max_active, kkt_stop, objective_stop, param_stop) {
+    .Call('selectiveInference_solve_QP_wide', PACKAGE = 'selectiveInference', X, bound, ridge_term, maxiter, theta, linear_func, gradient, X_theta, ever_active, nactive, kkt_tol, objective_tol, parameter_tol, max_active, kkt_stop, objective_stop, param_stop)
+}
+
+update1_ <- function(Q2, w, m, k) {
+    .Call('selectiveInference_update1_', PACKAGE = 'selectiveInference', Q2, w, m, k)
+}
+
+downdate1_ <- function(Q1, R, j0, m, n) {
+    .Call('selectiveInference_downdate1_', PACKAGE = 'selectiveInference', Q1, R, j0, m, n)
+}
+
+log_density_gaussian_ <- function(noise_scale, internal_linear, internal_state, optimization_linear, optimization_state, offset) {
+    .Call('selectiveInference_log_density_gaussian_', PACKAGE = 'selectiveInference', noise_scale, internal_linear, internal_state, optimization_linear, optimization_state, offset)
+}
+
+log_density_gaussian_conditional_ <- function(noise_scale, optimization_linear, optimization_state, offset) {
+    .Call('selectiveInference_log_density_gaussian_conditional_', PACKAGE = 'selectiveInference', noise_scale, optimization_linear, optimization_state, offset)
+}
+
+log_density_laplace_ <- function(noise_scale, internal_linear, internal_state, optimization_linear, optimization_state, offset) {
+    .Call('selectiveInference_log_density_laplace_', PACKAGE = 'selectiveInference', noise_scale, internal_linear, internal_state, optimization_linear, optimization_state, offset)
+}
+
+log_density_laplace_conditional_ <- function(noise_scale, optimization_linear, optimization_state, offset) {
+    .Call('selectiveInference_log_density_laplace_conditional_', PACKAGE = 'selectiveInference', noise_scale, optimization_linear, optimization_state, offset)
+}
+

selectiveInference/R/funs.randomized.R

@@ -258,7 +258,7 @@ conditional_density = function(noise_scale, lasso_soln) {
   return(lasso_soln)
 }
 
-randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term){
+randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term, level=0.9){
 
   n = nrow(X)
   p = ncol(X)
@@ -283,7 +283,7 @@ 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], 
@@ -300,20 +300,29 @@ randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term){
                                                      opt_transform,
                                                      target_transform,
                                                      observed_raw)
-      return(mean((target_sample<observed_target[i])*weights)/mean(weights))
+      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,34 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
 }
 
 
-collect_results = function(n,p,s, nsim=10){
-  rho=0.3
+collect_results = function(n,p,s, nsim=10, level=0.9){
+  rho=0.
   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
+    beta=data$beta
     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)
+    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
+      }
+      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()
 }