Skip to content

Commit fd0070a

Browse files
jonathan-taylorjonathan-taylor
committed
using TG base for all poly calcs
1 parent 9d41749 commit fd0070a

File tree

4 files changed

+66
-68
lines changed

4 files changed

+66
-68
lines changed

‎selectiveInference/R/funs.fixed.R

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -122,14 +122,19 @@ else{
122122
vj = vj / mj # Standardize (divide by norm of vj)
123123
sign[j] = sign(sum(vj*y))
124124
vj = sign[j] * vj
125-
a = poly.pval(y,G,u,vj,sigma,bits)
125+
126+
limits.info = TG.limits(y, -G, -u, vj, Sigma=diag(rep(sigma^2, n)))
127+
a = TG.pvalue.base(limits.info, bits=bits)
126128
pv[j] = a$pv
127129
vlo[j] = a$vlo * mj # Unstandardize (mult by norm of vj)
128130
vup[j] = a$vup * mj # Unstandardize (mult by norm of vj)
129131
vmat[j,] = vj * mj * sign[j] # Unstandardize (mult by norm of vj)
130132

131-
a = poly.int(y,G,u,vj,sigma,alpha,gridrange=gridrange,
132-
flip=(sign[j]==-1),bits=bits)
133+
a = TG.interval.base(limits.info,
134+
alpha=alpha,
135+
gridrange=gridrange,
136+
flip=(sign[j]==-1),
137+
bits=bits)
133138
ci[j,] = a$int * mj # Unstandardize (mult by norm of vj)
134139
tailarea[j,] = a$tailarea
135140
}

‎selectiveInference/R/funs.fs.R

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -299,15 +299,21 @@ fsInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic
299299
vj = vreg[j,]
300300
mj = sqrt(sum(vj^2))
301301
vj = vj / mj # Standardize (divide by norm of vj)
302-
a = poly.pval(y,Gj,uj,vj,sigma,bits)
302+
303+
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
304+
a = TG.pvalue.base(limits.info, bits=bits)
305+
303306
pv[j] = a$pv
304307
sxj = sx[vars[j]]
305308
vlo[j] = a$vlo * mj / sxj # Unstandardize (mult by norm of vj / sxj)
306309
vup[j] = a$vup * mj / sxj # Unstandardize (mult by norm of vj / sxj)
307310
vmat[j,] = vj * mj / sxj # Unstandardize (mult by norm of vj / sxj)
308311

309-
a = poly.int(y,Gj,uj,vj,sigma,alpha,gridrange=gridrange,
310-
flip=(sign[j]==-1),bits=bits)
312+
a = TG.interval.base(limits.info,
313+
alpha=alpha,
314+
gridrange=gridrange,
315+
flip=(sign[j]==-1),
316+
bits=bits)
311317
ci[j,] = a$int * mj / sxj # Unstandardize (mult by norm of vj / sxj)
312318
tailarea[j,] = a$tailarea
313319
}
@@ -349,15 +355,19 @@ fsInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","aic
349355
Gj = rbind(G,vj)
350356
uj = c(u,0)
351357

352-
a = poly.pval(y,Gj,uj,vj,sigma,bits)
358+
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
359+
a = TG.pvalue.base(limits.info, bits=bits)
353360
pv[j] = a$pv
354361
sxj = sx[vars[j]]
355362
vlo[j] = a$vlo * mj / sxj # Unstandardize (mult by norm of vj / sxj)
356363
vup[j] = a$vup * mj / sxj # Unstandardize (mult by norm of vj / sxj)
357364
vmat[j,] = vj * mj / sxj # Unstandardize (mult by norm of vj / sxj)
358365

359-
a = poly.int(y,Gj,uj,vj,sigma,alpha,gridrange=gridrange,
360-
flip=(sign[j]==-1),bits=bits)
366+
a = TG.interval.base(limits.info,
367+
alpha=alpha,
368+
gridrange=gridrange,
369+
flip=(sign[j]==-1),
370+
bits=bits)
361371
ci[j,] = a$int * mj / sxj # Unstandardize (mult by norm of vj / sxj)
362372
tailarea[j,] = a$tailarea
363373
}

‎selectiveInference/R/funs.inf.R

Lines changed: 26 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -1,48 +1,3 @@
1-
# Main p-value function
2-
3-
poly.pval <- function(y, G, u, v, sigma, bits=NULL) {
4-
z = sum(v*y)
5-
vv = sum(v^2)
6-
sd = sigma*sqrt(vv)
7-
8-
rho = G %*% v / vv
9-
vec = (u - G %*% y + rho*z) / rho
10-
vlo = suppressWarnings(max(vec[rho>0]))
11-
vup = suppressWarnings(min(vec[rho<0]))
12-
13-
pv = tnorm.surv(z,0,sd,vlo,vup,bits)
14-
return(list(pv=pv,vlo=vlo,vup=vup))
15-
}
16-
17-
# Main confidence interval function
18-
19-
poly.int <- function(y, G, u, v, sigma, alpha, gridrange=c(-100,100),
20-
gridpts=100, griddepth=2, flip=FALSE, bits=NULL) {
21-
22-
z = sum(v*y)
23-
vv = sum(v^2)
24-
sd = sigma*sqrt(vv)
25-
26-
rho = G %*% v / vv
27-
vec = (u - G %*% y + rho*z) / rho
28-
vlo = suppressWarnings(max(vec[rho>0]))
29-
vup = suppressWarnings(min(vec[rho<0]))
30-
31-
xg = seq(gridrange[1]*sd,gridrange[2]*sd,length=gridpts)
32-
fun = function(x) { tnorm.surv(z,x,sd,vlo,vup,bits) }
33-
34-
int = grid.search(xg,fun,alpha/2,1-alpha/2,gridpts,griddepth)
35-
tailarea = c(fun(int[1]),1-fun(int[2]))
36-
37-
if (flip) {
38-
int = -int[2:1]
39-
tailarea = tailarea[2:1]
40-
}
41-
42-
return(list(int=int,tailarea=tailarea))
43-
}
44-
45-
##############################
461

472
# Assuming that grid is in sorted order from smallest to largest,
483
# and vals are monotonically increasing function values over the
@@ -251,6 +206,8 @@ aicStop <- function(x, y, action, df, sigma, mult=2, ntimes=2) {
251206

252207
TG.limits = function(Z, A, b, eta, Sigma=NULL) {
253208

209+
target_estimate = sum(as.numeric(eta) * as.numeric(Z))
210+
254211
if (is.null(Sigma)) {
255212
Sigma = diag(rep(1, n))
256213
}
@@ -271,16 +228,14 @@ TG.limits = function(Z, A, b, eta, Sigma=NULL) {
271228
vup = suppressWarnings(min(vec[rho > 0]))
272229

273230
sd = sqrt(var_estimate)
274-
return(list(vlo=vlo, vup=vup, sd=sd))
231+
return(list(vlo=vlo, vup=vup, sd=sd, estimate=target_estimate))
275232
}
276233

277234
TG.pvalue = function(Z, A, b, eta, Sigma=NULL, null_value=0, bits=NULL) {
278235

279236
limits.info = TG.limits(Z, A, b, eta, Sigma)
280-
target_estimate = sum(as.numeric(eta) * as.numeric(Z))
281-
pv = tnorm.surv(target_estimate, null_value, limits.info$sd, limits.info$vlo, limits.info$vup, bits)
282237

283-
return(list(pv=pv, vlo=limits.info$vlo, vup=limits.info$vup, sd=limits.info$sd))
238+
return(TG.pvalue.base(limits.info, null_value=null_value, bits=bits))
284239
}
285240

286241
TG.interval = function(Z, A, b, eta, Sigma=NULL, alpha=0.1,
@@ -290,15 +245,29 @@ TG.interval = function(Z, A, b, eta, Sigma=NULL, alpha=0.1,
290245
flip=FALSE,
291246
bits=NULL) {
292247

248+
limits.info = TG.limits(Z, A, b, eta, Sigma)
249+
250+
return(TG.interval.base(limits.info,
251+
alpha=alpha,
252+
gridrange=gridrange,
253+
griddepth=griddepth,
254+
flip=flip,
255+
bits=bits))
256+
}
257+
258+
TG.interval.base = function(limits.info, alpha=0.1,
259+
gridrange=c(-100,100),
260+
gridpts=100,
261+
griddepth=2,
262+
flip=FALSE,
263+
bits=NULL) {
264+
293265
# compute sel intervals from poly lemmma, full version from Lee et al for full matrix Sigma
294266

295-
limits.info = TG.limits(Z, A, b, eta, Sigma)
296-
target_estimate = sum(as.numeric(eta) * as.numeric(Z))
297-
298267
param_grid = seq(gridrange[1] * limits.info$sd, gridrange[2] * limits.info$sd, length=gridpts)
299268

300269
pivot = function(param) {
301-
tnorm.surv(target_estimate, param, limits.info$sd, limits.info$vlo, limits.info$vup, bits)
270+
tnorm.surv(limits.info$estimate, param, limits.info$sd, limits.info$vlo, limits.info$vup, bits)
302271
}
303272

304273
interval = grid.search(param_grid, pivot, alpha/2, 1-alpha/2, gridpts, griddepth)
@@ -314,6 +283,10 @@ TG.interval = function(Z, A, b, eta, Sigma=NULL, alpha=0.1,
314283
tailarea=tailarea))
315284
}
316285

286+
TG.pvalue.base = function(limits.info, null_value=0, bits=NULL) {
287+
pv = tnorm.surv(limits.info$estimate, null_value, limits.info$sd, limits.info$vlo, limits.info$vup, bits)
288+
return(list(pv=pv, vlo=limits.info$vlo, vup=limits.info$vup, sd=limits.info$sd))
289+
}
317290

318291

319292
mydiag=function(x){

‎selectiveInference/R/funs.lar.R

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -379,15 +379,20 @@ larInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","ai
379379
vj = vreg[j,]
380380
mj = sqrt(sum(vj^2))
381381
vj = vj / mj # Standardize (divide by norm of vj)
382-
a = poly.pval(y,Gj,uj,vj,sigma,bits)
382+
383+
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
384+
a = TG.pvalue.lowlevel(limits.info, bits=bits)
383385
pv[j] = a$pv
384386
sxj = sx[vars[j]]
385387
vlo[j] = a$vlo * mj / sxj # Unstandardize (mult by norm of vj / sxj)
386388
vup[j] = a$vup * mj / sxj # Unstandardize (mult by norm of vj)
387389
vmat[j,] = vj * mj / sxj # Unstandardize (mult by norm of vj / sxj)
388390

389-
a = poly.int(y,Gj,uj,vj,sigma,alpha,gridrange=gridrange,
390-
flip=(sign[j]==-1),bits=bits)
391+
a = TG.interval.lowlevel(limits.info,
392+
alpha=alpha,
393+
gridrange=gridrange,
394+
flip=(sign[j]==-1),
395+
bits=bits)
391396
ci[j,] = a$int * mj / sxj # Unstandardize (mult by norm of vj / sxj)
392397
tailarea[j,] = a$tailarea
393398

@@ -433,15 +438,20 @@ larInf <- function(obj, sigma=NULL, alpha=0.1, k=NULL, type=c("active","all","ai
433438
Gj = rbind(G,vj)
434439
uj = c(u,0)
435440

436-
a = poly.pval(y,Gj,uj,vj,sigma,bits)
441+
limits.info = TG.limits(y, -Gj, -uj, vj, Sigma=diag(rep(sigma^2, n)))
442+
a = TG.pvalue.lowlevel(limits.info, bits=bits)
443+
437444
pv[j] = a$pv
438445
sxj = sx[vars[j]]
439446
vlo[j] = a$vlo * mj / sxj # Unstandardize (mult by norm of vj / sxj)
440447
vup[j] = a$vup * mj / sxj # Unstandardize (mult by norm of vj / sxj)
441448
vmat[j,] = vj * mj / sxj # Unstandardize (mult by norm of vj / sxj)
442449

443-
a = poly.int(y,Gj,uj,vj,sigma,alpha,gridrange=gridrange,
444-
flip=(sign[j]==-1),bits=bits)
450+
a = TG.interval.lowlevel(limits.info,
451+
alpha=alpha,
452+
gridrange=gridrange,
453+
flip=(sign[j]==-1),
454+
bits=bits)
445455
ci[j,] = a$int * mj / sxj # Unstandardize (mult by norm of vj / sxj)
446456
tailarea[j,] = a$tailarea
447457
}

0 commit comments

Comments
 (0)