Skip to content

Commit d1ebca3

Browse files
jonathan-taylorjonathan-taylor
committed
BF: noise_scale was NULL
1 parent 46856e1 commit d1ebca3

File tree

4 files changed

+67
-36
lines changed

4 files changed

+67
-36
lines changed

‎C-software

‎selectiveInference/R/funs.randomized.R

Lines changed: 63 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -65,7 +65,6 @@ randomizedLasso = function(X,
6565
ever_active = rep(0, p)
6666
nactive = as.integer(0)
6767

68-
print('here')
6968
result = solve_QP_wide(X, # design matrix
7069
lam / n, # vector of Lagrange multipliers
7170
ridge_term / n, # ridge_term
@@ -86,7 +85,6 @@ randomizedLasso = function(X,
8685

8786
sign_soln = sign(result$soln)
8887

89-
print('now')
9088
unpenalized = lam == 0
9189
active = (!unpenalized) & (sign_soln != 0)
9290
inactive = (!unpenalized) & (sign_soln == 0)
@@ -156,16 +154,17 @@ randomizedLasso = function(X,
156154

157155
log_optimization_density = function(opt_state
158156
) {
159-
160-
161157
if ((sum(abs(opt_state[(inactive_start + 1):p]) > inactive_lam) > 0) ||
162158
(sum(opt_state[(active_start+1):inactive_start] < 0) > 0)) {
163159
return(-Inf)
164160
}
165-
D = log_density_gaussian_conditional_(noise_scale,
166-
opt_transform$linear_term,
167-
as.matrix(opt_state),
168-
observed_raw+opt_transform$offset_term)
161+
162+
A = opt_transform$linear_term %*% opt_state + observed_raw + opt_transform$offset_term
163+
D = -apply(A^2, 2, sum) / noise_scale^2
164+
# D = log_density_gaussian_conditional_(noise_scale,
165+
# opt_transform$linear_term,
166+
# as.matrix(opt_state),
167+
# observed_raw + opt_transform$offset_term)
169168
return(D)
170169
}
171170

@@ -177,7 +176,8 @@ randomizedLasso = function(X,
177176
internal_transform = internal_transform,
178177
log_optimization_density = log_optimization_density,
179178
observed_opt_state = observed_opt_state,
180-
observed_raw = observed_raw
179+
observed_raw = observed_raw,
180+
noise_scale = noise_scale
181181
))
182182

183183
}
@@ -222,17 +222,31 @@ importance_weight = function(noise_scale,
222222
target_transform,
223223
observed_raw) {
224224

225-
log_num = log_density_gaussian_(noise_scale,
226-
target_transform$linear_term,
227-
as.matrix(target_sample),
228-
opt_transform$linear_term,
229-
as.matrix(opt_sample),
230-
target_transform$offset_term + opt_transform$offset_term)
231-
232-
log_den = log_density_gaussian_conditional_(noise_scale,
233-
opt_transform$linear_term,
234-
as.matrix(opt_sample),
235-
observed_raw+opt_transform$offset_term)
225+
use_C_code = FALSE
226+
if (!use_C_code) {
227+
A = (opt_transform$linear_term %*% opt_sample +
228+
target_transform$linear_term %*% target_sample)
229+
A = apply(A, 2, function(x) {return(x + target_transform$offset_term + opt_transform$offset_term)})
230+
log_num = -apply(A^2, 2, sum) / noise_scale^2
231+
} else {
232+
log_num = log_density_gaussian_(noise_scale,
233+
target_transform$linear_term,
234+
as.matrix(target_sample),
235+
opt_transform$linear_term,
236+
as.matrix(opt_sample),
237+
target_transform$offset_term + opt_transform$offset_term)
238+
}
239+
240+
if (!use_C_code) {
241+
A = opt_transform$linear_term %*% opt_sample
242+
A = apply(A, 2, function(x) {return(x + observed_raw + opt_transform$offset_term)})
243+
log_den = -apply(A^2, 2, sum) / noise_scale^2
244+
} else {
245+
log_den = log_density_gaussian_conditional_(noise_scale,
246+
opt_transform$linear_term,
247+
as.matrix(opt_sample),
248+
observed_raw+opt_transform$offset_term)
249+
}
236250
W = log_num - log_den
237251
W = W - max(W)
238252
return(exp(W))
@@ -263,11 +277,19 @@ conditional_density = function(noise_scale, lasso_soln) {
263277
if (sum(opt_state < 0) > 0) {
264278
return(-Inf)
265279
}
266-
D = log_density_gaussian_conditional_(noise_scale,
267-
reduced_B,
268-
as.matrix(opt_state),
269-
reduced_beta_offset)
270-
return(D)
280+
281+
use_C_code = FALSE
282+
if (!use_C_code) {
283+
A = reduced_B %*% as.matrix(opt_state) + reduced_beta_offset
284+
A = apply(A, 2, function(x) {x + reduced_beta_offset})
285+
log_den = -apply(A^2, 2, sum) / noise_scale^2
286+
} else {
287+
log_den = log_density_gaussian_conditional_(noise_scale,
288+
reduced_B,
289+
as.matrix(opt_state),
290+
reduced_beta_offset)
291+
}
292+
return(log_den)
271293
}
272294
lasso_soln$log_optimization_density = log_condl_optimization_density
273295
lasso_soln$observed_opt_state = observed_opt_state[1:nactive]
@@ -282,7 +304,9 @@ randomizedLassoInf = function(X,
282304
noise_scale=NULL,
283305
ridge_term=NULL,
284306
condition_subgrad=TRUE,
285-
level=0.9) {
307+
level=0.9,
308+
nsample=10000,
309+
burnin=2000) {
286310

287311
n = nrow(X)
288312
p = ncol(X)
@@ -291,14 +315,20 @@ randomizedLassoInf = function(X,
291315
inactive_set = lasso_soln$inactive_set
292316
nactive = length(active_set)
293317

294-
if (condition_subgrad==TRUE){
295-
lasso_soln=conditional_density(noise_scale,lasso_soln)
296-
}
318+
noise_scale = lasso_soln$noise_scale # set to default value in randomizedLasso
319+
320+
if (condition_subgrad==TRUE){
321+
lasso_soln=conditional_density(noise_scale, lasso_soln)
322+
}
297323

298324
dim = length(lasso_soln$observed_opt_state)
299325
print(paste("chain dim", dim))
300-
S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=10000)
301-
opt_samples = S$samples[2001:10000,]
326+
327+
# print(lasso_soln)
328+
329+
330+
S = sample_opt_variables(lasso_soln, jump_scale=rep(1/sqrt(n), dim), nsample=nsample)
331+
opt_samples = S$samples[(burnin+1):nsample,]
302332
print(paste("dim opt samples", toString(dim(opt_samples))))
303333

304334
X_E = X[, active_set]
@@ -331,12 +361,12 @@ randomizedLassoInf = function(X,
331361

332362
pivot = function(candidate){
333363
weights = importance_weight(noise_scale,
334-
t(as.matrix(target_sample)) + candidate,
364+
t(as.matrix(target_sample) + candidate),
335365
t(opt_samples),
336366
opt_transform,
337367
target_transform,
338368
observed_raw)
339-
return(mean((target_sample+candidate<observed_target[i])*weights)/mean(weights))
369+
return(mean((target_sample + candidate < observed_target[i]) * weights)/mean(weights))
340370
}
341371
rootU = function(candidate){
342372
return (pivot(observed_target[i]+candidate)-(1-level)/2)

‎selectiveInference/src/Rcpp-randomized.cpp

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#include <Rcpp.h> // need to include the main Rcpp header file
22
#include <randomized_lasso.h> // where densities are defined
3-
3+
#include <stdio.h>
44
// [[Rcpp::export]]
55
Rcpp::NumericVector log_density_gaussian_(double noise_scale, // Scale of randomization
66
Rcpp::NumericMatrix internal_linear, // A_D -- linear part for data
@@ -53,6 +53,7 @@ Rcpp::NumericVector log_density_gaussian_conditional_(double noise_scale,
5353

5454
int ipt;
5555
for (ipt=0; ipt<npt; ipt++) {
56+
// fprintf(stderr, "here 5: %f %d %d %d %d\n", noise_scale, ndim, noptimization, ipt, npt);
5657
result[ipt] = log_density_gaussian_conditional(noise_scale,
5758
ndim,
5859
noptimization,

‎tests/randomized/test_instances.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ collect_results = function(n,p,s, nsim=100, level=0.9){
3636
X=data$X
3737
y=data$y
3838
beta=data$beta
39-
result = selectiveInference:::randomizedLassoInf(X, y, lam, level=level)
39+
result = selectiveInference:::randomizedLassoInf(X, y, lam, level=level, burnin=2000, nsample=4000)
4040
true_beta = beta[result$active_set]
4141
coverage = rep(0, nrow(result$ci))
4242
for (i in 1:nrow(result$ci)){

0 commit comments

Comments
 (0)