|
| 1 | +#library(devtools) |
| 2 | +#devtools::install_github('jonathan-taylor/R-selective/selectiveInference') |
| 3 | +library(selectiveInference, lib.loc='/Users/Jelena/anaconda/lib/R/library') |
| 4 | + |
| 5 | + |
| 6 | +gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA, |
| 7 | + random_signs=TRUE, scale=TRUE, center=TRUE, seed=NA){ |
| 8 | + if (!is.na(seed)){ |
| 9 | + set.seed(seed) |
| 10 | + } |
| 11 | + |
| 12 | + if (is.na(X)){ |
| 13 | + X = sqrt(1-rho)*matrix(rnorm(n*p),n) + sqrt(rho)*matrix(rep(rnorm(n), p), nrow = n) |
| 14 | + X = scale(X)/sqrt(n) |
| 15 | + } |
| 16 | + beta = rep(0, p) |
| 17 | + if (s>0){ |
| 18 | + beta[1:s] = seq(3, 6, length.out=s) |
| 19 | + } |
| 20 | + beta = sample(beta) |
| 21 | + if (random_signs==TRUE & s>0){ |
| 22 | + signs = sample(c(-1,1), s, replace = TRUE) |
| 23 | + beta = beta * signs |
| 24 | + } |
| 25 | + y = X %*% beta + rnorm(n)*sigma |
| 26 | + result = list(X=X,y=y,beta=beta) |
| 27 | + return(result) |
| 28 | +} |
| 29 | + |
| 30 | + |
| 31 | +run_instance = function(n, p, s){ |
| 32 | + rho=0.3 |
| 33 | + lam=1.3 |
| 34 | + sigma=1 |
| 35 | + data = gaussian_instance(n=n,p=p,s=s, rho=rho, sigma=sigma) |
| 36 | + X=data$X |
| 37 | + print(dim(X)) |
| 38 | + y=data$y |
| 39 | + ridge_term=sd(y)/sqrt(n) |
| 40 | + noise_scale = sd(y)/2 |
| 41 | + lasso_soln=selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term) |
| 42 | + |
| 43 | + active_set = lasso_soln$active_set |
| 44 | + inactive_set = lasso_soln$inactive_set |
| 45 | + observed_raw = lasso_soln$observed_raw |
| 46 | + opt_linear = lasso_soln$optimization_transform$linear_term |
| 47 | + opt_offset = lasso_soln$optimization_transform$offset_term |
| 48 | + observed_opt_state = lasso_soln$observed_opt_state |
| 49 | + |
| 50 | + nactive = length(active_set) |
| 51 | + print(paste("nactive",nactive)) |
| 52 | + B = opt_linear[,1:nactive] |
| 53 | + beta_offset = observed_raw+opt_offset |
| 54 | + if (nactive<p){ |
| 55 | + U=opt_linear[,(nactive+1):p] |
| 56 | + beta_offset =+U %*% observed_opt_state[(nactive+1):p] |
| 57 | + } |
| 58 | + opt_transform = list(linear_term=B, offset_term = beta_offset) |
| 59 | + reduced_B = chol(t(B) %*% B) |
| 60 | + reduced_beta_offset = solve(t(reduced_B)) %*% (t(B) %*% beta_offset) |
| 61 | + |
| 62 | + log_condl_optimization_density = function(opt_state) { |
| 63 | + |
| 64 | + if (sum(opt_state < 0) > 0) { |
| 65 | + return(-Inf) |
| 66 | + } |
| 67 | + D = selectiveInference:::log_density_gaussian_conditional_(noise_scale, |
| 68 | + reduced_B, |
| 69 | + as.matrix(observed_opt_state[1:nactive]), |
| 70 | + reduced_beta_offset) |
| 71 | + return(D) |
| 72 | + } |
| 73 | + lasso_soln$log_optimization_density = log_condl_optimization_density |
| 74 | + lasso_soln$observed_opt_state = observed_opt_state[1:nactive] |
| 75 | + S = selectiveInference:::sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), nactive), nsample=10000) |
| 76 | + beta_samples = S$samples[2001:10000,] |
| 77 | + print(paste("dim beta samples", dim(beta_samples))) |
| 78 | + |
| 79 | + X_E=X[, active_set] |
| 80 | + X_minusE=X[, inactive_set] |
| 81 | + target_cov = solve(t(X_E)%*%X_E)*sigma^2 |
| 82 | + cov_target_internal = rbind(target_cov, matrix(0, nrow=p-nactive, ncol=nactive)) * sigma^2 |
| 83 | + observed_target = solve(t(X_E) %*% X_E) %*% t(X_E) %*% y |
| 84 | + observed_internal = c(observed_target, t(X_minusE) %*% (y-X_E%*% observed_target)) |
| 85 | + internal_transform = lasso_soln$internal_transform |
| 86 | + |
| 87 | + pivots = rep(0, nactive) |
| 88 | + for (i in 1:nactive){ |
| 89 | + target_transform = selectiveInference:::linear_decomposition(observed_target[i], |
| 90 | + observed_internal, |
| 91 | + target_cov[i,i], |
| 92 | + cov_target_internal[,i], |
| 93 | + internal_transform) |
| 94 | + target_sample = rnorm(nrow(beta_samples)) * sqrt(target_cov[i,i]) |
| 95 | + |
| 96 | + weights = selectiveInference:::importance_weight(noise_scale, |
| 97 | + t(as.matrix(target_sample)), |
| 98 | + t(beta_samples), |
| 99 | + opt_transform, |
| 100 | + target_transform, |
| 101 | + observed_raw) |
| 102 | + |
| 103 | + pivots[i] = mean((target_sample<observed_target[i])*weights)/mean(weights) |
| 104 | + print(pivots[i]) |
| 105 | + } |
| 106 | + |
| 107 | + return(pivots) |
| 108 | +} |
| 109 | + |
| 110 | +collect_instances = function(n,p,s, nsim=1){ |
| 111 | + |
| 112 | + for (i in 1:nsim){ |
| 113 | + result = run_instance(n,p,s) |
| 114 | + } |
| 115 | +} |
| 116 | + |
| 117 | + |
| 118 | +collect_instances(n=100, p=20, s=0) |
| 119 | + |
| 120 | + |
| 121 | + |
0 commit comments