Skip to content

Commit a1ad181

Browse files
author
Jelena Markovic
authored
Merge pull request #1 from jonathan-taylor/master
BF: segfault should be gone now
2 parents 8f99570 + 95b98c2 commit a1ad181

File tree

5 files changed

+105
-151
lines changed

5 files changed

+105
-151
lines changed

‎C-software

‎selectiveInference/R/RcppExports.R

Lines changed: 0 additions & 35 deletions
This file was deleted.

‎selectiveInference/R/funs.randomized.R

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -221,3 +221,99 @@ importance_weight = function(noise_scale,
221221
return(exp(W))
222222
}
223223

224+
conditional_density = function(noise_scale, lasso_soln) {
225+
226+
active_set = lasso_soln$active_set
227+
observed_raw = lasso_soln$observed_raw
228+
opt_linear = lasso_soln$optimization_transform$linear_term
229+
opt_offset = lasso_soln$optimization_transform$offset_term
230+
observed_opt_state = lasso_soln$observed_opt_state
231+
232+
nactive = length(active_set)
233+
B = opt_linear[,1:nactive]
234+
beta_offset = opt_offset
235+
p=length(observed_opt_state)
236+
if (nactive<p){
237+
beta_offset = beta_offset+(opt_linear[,(nactive+1):p] %*% observed_opt_state[(nactive+1):p])
238+
}
239+
opt_transform = list(linear_term=B,
240+
offset_term = beta_offset)
241+
reduced_B = chol(t(B) %*% B)
242+
beta_offset = beta_offset + observed_raw
243+
reduced_beta_offset = solve(t(reduced_B)) %*% (t(B) %*% beta_offset)
244+
245+
log_condl_optimization_density = function(opt_state) {
246+
if (sum(opt_state < 0) > 0) {
247+
return(-Inf)
248+
}
249+
D = selectiveInference:::log_density_gaussian_conditional_(noise_scale,
250+
reduced_B,
251+
as.matrix(opt_state),
252+
reduced_beta_offset)
253+
return(D)
254+
}
255+
lasso_soln$log_optimization_density = log_condl_optimization_density
256+
lasso_soln$observed_opt_state = observed_opt_state[1:nactive]
257+
lasso_soln$optimization_transform = opt_transform
258+
return(lasso_soln)
259+
}
260+
261+
randomized_inference = function(X, y, sigma, lam, noise_scale, ridge_term){
262+
263+
n = nrow(X)
264+
p = ncol(X)
265+
lasso_soln = selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
266+
active_set = lasso_soln$active_set
267+
inactive_set = lasso_soln$inactive_set
268+
nactive = length(active_set)
269+
270+
dim = length(lasso_soln$observed_opt_state)
271+
print(paste("chain dim", dim))
272+
S = selectiveInference:::sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=10000)
273+
opt_samples = S$samples[2001:10000,]
274+
print(paste("dim opt samples", toString(dim(opt_samples))))
275+
276+
X_E = X[, active_set]
277+
X_minusE = X[, inactive_set]
278+
target_cov = solve(t(X_E) %*% X_E)*sigma^2
279+
cov_target_internal = rbind(target_cov, matrix(0, nrow=p-nactive, ncol=nactive))
280+
observed_target = solve(t(X_E) %*% X_E) %*% t(X_E) %*% y
281+
observed_internal = c(observed_target, t(X_minusE) %*% (y-X_E%*% observed_target))
282+
internal_transform = lasso_soln$internal_transform
283+
opt_transform = lasso_soln$optimization_transform
284+
observed_raw = lasso_soln$observed_raw
285+
286+
pvalus = rep(0, nactive)
287+
ci = matrix(0, nactive, 2)
288+
for (i in 1:nactive){
289+
target_transform = selectiveInference:::linear_decomposition(observed_target[i],
290+
observed_internal,
291+
target_cov[i,i],
292+
cov_target_internal[,i],
293+
internal_transform)
294+
target_sample = rnorm(nrow(opt_samples)) * sqrt(target_cov[i,i])
295+
296+
pivot = function(candidate){
297+
weights = selectiveInference:::importance_weight(noise_scale,
298+
t(as.matrix(target_sample)) + candidate,
299+
t(opt_samples),
300+
opt_transform,
301+
target_transform,
302+
observed_raw)
303+
return(mean((target_sample<observed_target[i])*weights)/mean(weights))
304+
}
305+
level = 0.9
306+
rootU = function(candidate){
307+
return (pivot(observed_target[i]+candidate)-(1-level)/2)
308+
}
309+
rootL = function(candidate){
310+
return (pivot(observed_target[i]+candidate)-(1+level)/2)
311+
}
312+
pvalues[i] = pivot(0)
313+
line_min = -10*sd(target_sample)
314+
line_max = 10*sd(target_sample)
315+
ci[i,1] = uniroot(rootU, c(line_min, line_max))$root+observed_target[i]
316+
ci[i,2] = uniroot(rootL, c(line_min, line_max))$root+observed_target[i]
317+
}
318+
return(list(pvalues=pvalues, ci=ci))
319+
}

‎selectiveInference/src/Rcpp-debias.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -160,8 +160,8 @@ Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X,
160160
(double *) X.begin(),
161161
(double *) linear_func.begin(),
162162
(int *) need_update.begin(),
163-
nfeature,
164163
ncase,
164+
nfeature,
165165
(double *) bound.begin(),
166166
ridge_term,
167167
kkt_tol);
@@ -175,8 +175,8 @@ Rcpp::List solve_QP_wide(Rcpp::NumericMatrix X,
175175
(double *) X.begin(),
176176
(double *) linear_func.begin(),
177177
(int *) need_update.begin(),
178-
nfeature,
179-
ncase);
178+
ncase,
179+
nfeature);
180180

181181
return(Rcpp::List::create(Rcpp::Named("soln") = theta,
182182
Rcpp::Named("gradient") = gradient,

‎tests/randomized/test_instances.R

Lines changed: 5 additions & 112 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,4 @@
1-
#library(devtools)
2-
#devtools::install_github('jonathan-taylor/R-selective/selectiveInference')
3-
library(selectiveInference, lib.loc='/Users/Jelena/anaconda/lib/R/library')
4-
1+
library(selectiveInference)
52

63
gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
74
random_signs=TRUE, scale=TRUE, center=TRUE, seed=NA){
@@ -10,7 +7,7 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
107
}
118

129
if (is.na(X)){
13-
X = sqrt(1-rho)*matrix(rnorm(n*p),n) + sqrt(rho)*matrix(rep(rnorm(n), p), nrow = n)
10+
X = sqrt(1-rho)*matrix(rnorm(n*p),n, p) + sqrt(rho)*matrix(rep(rnorm(n), p), nrow = n)
1411
X = scale(X)/sqrt(n)
1512
}
1613
beta = rep(0, p)
@@ -27,119 +24,15 @@ gaussian_instance = function(n, p, s, sigma=1, rho=0, signal=6, X=NA,
2724
return(result)
2825
}
2926

30-
conditional_density = function(noise_scale, lasso_soln){
31-
32-
active_set = lasso_soln$active_set
33-
observed_raw = lasso_soln$observed_raw
34-
opt_linear = lasso_soln$optimization_transform$linear_term
35-
opt_offset = lasso_soln$optimization_transform$offset_term
36-
observed_opt_state = lasso_soln$observed_opt_state
37-
38-
nactive = length(active_set)
39-
B = opt_linear[,1:nactive]
40-
beta_offset = opt_offset
41-
p=length(observed_opt_state)
42-
if (nactive<p){
43-
beta_offset = beta_offset+(opt_linear[,(nactive+1):p] %*% observed_opt_state[(nactive+1):p])
44-
}
45-
opt_transform = list(linear_term=B,
46-
offset_term = beta_offset)
47-
reduced_B = chol(t(B) %*% B)
48-
beta_offset = beta_offset+observed_raw
49-
reduced_beta_offset = solve(t(reduced_B)) %*% (t(B) %*% beta_offset)
50-
51-
log_condl_optimization_density = function(opt_state) {
52-
if (sum(opt_state < 0) > 0) {
53-
return(-Inf)
54-
}
55-
D = selectiveInference:::log_density_gaussian_conditional_(noise_scale,
56-
reduced_B,
57-
as.matrix(opt_state),
58-
reduced_beta_offset)
59-
return(D)
60-
}
61-
lasso_soln$log_optimization_density = log_condl_optimization_density
62-
lasso_soln$observed_opt_state = observed_opt_state[1:nactive]
63-
lasso_soln$optimization_transform = opt_transform
64-
return(lasso_soln)
65-
}
66-
67-
68-
randomized_inference = function(X,y,sigma, lam, noise_scale, ridge_term){
69-
n=nrow(X)
70-
p=ncol(X)
71-
lasso_soln=selectiveInference:::randomizedLASSO(X, y, lam, noise_scale, ridge_term)
72-
active_set = lasso_soln$active_set
73-
inactive_set = lasso_soln$inactive_set
74-
nactive = length(active_set)
75-
print(paste("nactive", nactive))
76-
77-
#lasso_soln = conditional_density(noise_scale, lasso_soln)
78-
79-
dim=length(lasso_soln$observed_opt_state)
80-
print(paste("chain dim", dim))
81-
S = selectiveInference:::sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=10000)
82-
opt_samples = S$samples[2001:10000,]
83-
print(paste("dim opt samples", toString(dim(opt_samples))))
84-
85-
X_E=X[, active_set]
86-
X_minusE=X[, inactive_set]
87-
target_cov = solve(t(X_E) %*% X_E)*sigma^2
88-
cov_target_internal = rbind(target_cov, matrix(0, nrow=p-nactive, ncol=nactive))
89-
observed_target = solve(t(X_E) %*% X_E) %*% t(X_E) %*% y
90-
observed_internal = c(observed_target, t(X_minusE) %*% (y-X_E%*% observed_target))
91-
internal_transform = lasso_soln$internal_transform
92-
opt_transform = lasso_soln$optimization_transform
93-
observed_raw = lasso_soln$observed_raw
94-
95-
pivots = rep(0, nactive)
96-
ci = matrix(0, nactive, 2)
97-
for (i in 1:nactive){
98-
target_transform = selectiveInference:::linear_decomposition(observed_target[i],
99-
observed_internal,
100-
target_cov[i,i],
101-
cov_target_internal[,i],
102-
internal_transform)
103-
target_sample = rnorm(nrow(opt_samples)) * sqrt(target_cov[i,i])
104-
105-
pivot = function(candidate){
106-
weights = selectiveInference:::importance_weight(noise_scale,
107-
t(as.matrix(target_sample))+candidate,
108-
t(opt_samples),
109-
opt_transform,
110-
target_transform,
111-
observed_raw)
112-
return(mean((target_sample<observed_target[i])*weights)/mean(weights))
113-
}
114-
level = 0.9
115-
rootU = function(candidate){
116-
return (pivot(observed_target[i]+candidate)-(1-level)/2)
117-
}
118-
rootL = function(candidate){
119-
return (pivot(observed_target[i]+candidate)-(1+level)/2)
120-
}
121-
pivots[i] = pivot(0)
122-
line_min = -10*sd(target_sample)
123-
line_max = 10*sd(target_sample)
124-
ci[i,1] = uniroot(rootU, c(line_min, line_max))$root+observed_target[i]
125-
ci[i,2] = uniroot(rootL,c(line_min, line_max))$root+observed_target[i]
126-
}
127-
print(paste("pivots", toString(pivots)))
128-
for (i in 1:nactive){
129-
print(paste("CIs", toString(ci[i,])))
130-
}
131-
return(list(pivots=pivots, ci=ci))
132-
}
13327

134-
collect_results = function(n,p,s, nsim=1){
28+
collect_results = function(n,p,s, nsim=10){
13529
rho=0.3
13630
lam=1.
13731
sigma=1
138-
sample_pivots = NULL
32+
sample_pivots = c()
13933
for (i in 1:nsim){
14034
data = gaussian_instance(n=n,p=p,s=s, rho=rho, sigma=sigma)
14135
X=data$X
142-
print(dim(X))
14336
y=data$y
14437
ridge_term=sd(y)/sqrt(n)
14538
noise_scale = sd(y)/2
@@ -148,7 +41,7 @@ collect_results = function(n,p,s, nsim=1){
14841
#lam = 20 / sqrt(n)
14942
#noise_scale = 0.01 * sqrt(n)
15043
#ridge_term = .1 / sqrt(n)
151-
result = randomized_inference(X,y,sigma,lam,noise_scale,ridge_term)
44+
result = selectiveInference:::randomized_inference(X,y,sigma,lam,noise_scale,ridge_term)
15245
sample_pivots = c(sample_pivots, result$pivots)
15346
}
15447

0 commit comments

Comments
 (0)