@@ -562,13 +562,17 @@ print.ROSI <- function(x, ...) {
562
562
563
563
564
564
compute_coverage = function (ci , beta ){
565
+ print(ci )
566
+ print(beta )
565
567
nactive = length(beta )
566
- coverage_vector = rep(0 , nactive )
568
+ coverage_vector = rep(NA , nactive )
567
569
for (i in 1 : nactive ){
568
- if (beta [i ]> = ci [1 ,i ] && beta [i ]< = ci [2 ,i ]){
570
+ if (! is.na(ci [i ,1 ]) & ! is.na(ci [i ,2 ])) {
571
+ if (beta [i ]> = ci [i ,1 ] && beta [i ]< = ci [i ,2 ]){
569
572
coverage_vector [i ]= 1
570
573
}
571
574
}
575
+ }
572
576
return (coverage_vector )
573
577
}
574
578
@@ -590,91 +594,3 @@ family_label = function(loss){
590
594
}
591
595
592
596
593
- # Unused
594
-
595
- selective_CI = function (observed ,
596
- variance ,
597
- sigma_est ,
598
- center ,
599
- radius ,
600
- alpha = 0.1 ,
601
- gridrange = c(- 20 ,20 ),
602
- gridpts = 10000 ,
603
- griddepth = 2 ) {
604
-
605
- pivot = function (param ){
606
- return (test_TG(param ,
607
- observed ,
608
- variance ,
609
- sigma_est ,
610
- center ,
611
- radius ,
612
- alt = " upper" ))
613
- }
614
- st.error = sqrt(variance )
615
- param_grid = seq(observed + gridrange [1 ]* st.error ,
616
- observed + gridrange [2 ]* st.error ,
617
- length = gridpts )
618
- interval = grid.search(param_grid ,
619
- pivot ,
620
- alpha / 2 ,
621
- 1 - alpha / 2 ,
622
- gridpts ,
623
- griddepth )
624
- return (interval )
625
- }
626
-
627
- # the pvalue if prob(Z>obs given |sigma_est^2/variance * Z+center|>radius)
628
- # where Z~N(param, variance)
629
- test_TG = function (param ,
630
- observed ,
631
- variance ,
632
- sigma_est ,
633
- center ,
634
- radius ,
635
- alt ) {
636
- st.error = sqrt(variance )
637
- lower = variance * (- center - radius )/ (sigma_est ^ 2 )
638
- upper = variance * (- center + radius )/ (sigma_est ^ 2 )
639
- if (observed < = lower ){
640
- case = 1
641
- num = (pnorm(upper ,
642
- mean = param ,
643
- sd = st.error ,
644
- lower.tail = FALSE ) +
645
- pnorm(lower ,
646
- mean = param ,
647
- sd = st.error ) -
648
- pnorm(observed , mean = param , sd = st.error ))
649
- } else if (observed > = upper ){
650
- case = 2
651
- num = pnorm(observed ,
652
- mean = param ,
653
- sd = st.error ,
654
- lower.tail = FALSE )
655
- } else {
656
- case = 3
657
- num = pnorm(upper ,
658
- mean = param ,
659
- sd = st.error ,
660
- lower.tail = FALSE )
661
- return (NULL )
662
- }
663
- den = (pnorm(upper ,
664
- mean = param ,
665
- sd = st.error ,
666
- lower.tail = FALSE ) +
667
- pnorm(lower ,
668
- mean = param ,
669
- sd = st.error ))
670
- pivot = num / den
671
-
672
- if (alt == " two-sided" ){
673
- return (2 * pmin(pivot , 1 - pivot ))
674
- } else if (alt == " upper" ) {
675
- return (pivot )
676
- } else if (alt == " lower" ){
677
- return (1 - pivot )
678
- }
679
- }
680
-
0 commit comments