# Begin constants # k00<- (1.0-exp(-0.5))/3.0 k0<- 2.0*(1-k00) k1<- 0.5*(exp(-0.5)+exp(-2.0)) k11<- 2.0*exp(-0.5)-exp(-2.0) k12<- (exp(-2.0)-exp(-0.5))/2.0 kg1<- -k11/(2.0*k12) k2<- exp(-2.0)/2.0 cp0<- 0.5*(3*k00)^1.5 ym0<- 2.0/sqrt(3.0*k00) # End constants# # Begin functions for integral constants # c00<- function(a,b){ (b-a)-k00*(b^3-a^3) } c11<- function(a,b){ k11*(b-a)+k12*(b^2-a^2) } c22<- function(a,b){ exp(-a^2/2.0)*(1.0-exp((a^2-b^2)/2))/2.0 } # End functions for integral constants # # Begin within-region generators # genregion0<- function(a,b,c){ repeat{ q<- a^3-(a+runif(1)*c)/k00 C<- cp0*q x<- ym0*cos((acos(C)+4*pi)/3.0) if(runif(1)*(1.0-3.0*k00*x^2) < exp(-x^2/2)) break } x } genregion1<- function(a,b,c){ repeat{ x<- kg1+sqrt((k11+2.0*k12*a)^2+4.0*k12*runif(1)*c)/(2.0*k12); if(runif(1)*(k11+2.0*k12*x) < exp(-x^2/2)) break } x } genregion2<- function(a,b,c){ repeat{x<- sqrt(a^2-2*log(1-runif(1)*(1-exp(-(b^2-a^2)/2)))); if(runif(1)<= a/x) break} x } # End within-region generators # # A function to test regional generators # testr<- function(a,b,c,region,n){ x<- 1:n; i<- 1; while(i<=n){x[i]<-generator[[region+1]](a,b,c); i<- i+1} hist(x,prob=TRUE) ulim=b; if(b==Inf) ulim<- max(x); xx<- seq(from=a,to=ulim,length.out=100) lines(xx,dnorm(xx)/(pnorm(b)-pnorm(a))) } # Begin case list and case functions # case<- list() case[[11]]<- genregion2 case[[10]]<- function(a,b){ p1<- c11(a,2); p2<- c22(2,b); if(runif(1)*(p1+p2) < p1) sample<- genregion1(a,2,p1) else sample<- genregion2(2,b,p2) sample } case[[9]]<- function(a,b){ p1<- c00(a,1); p2<- k1; p3<- c22(2,b); pickregion<- runif(1)*(p1+p2+p3); if(pickregion < p1) sample<- genregion0(a,1,p1) else if(pickregion < p1+p2) sample<- genregion1(1,2,p2) else sample<- genregion2(2,b,p3) sample } case[[8]]<- function(a,b){ p1<- c11(1,-a); p2<- k0; p3<- c11(1,2); p4<- c22(2,b); pickregion<- runif(1)*(p1+p2+p3+p4); if(pickregion < p1) sample<- -genregion1(1,-a,p1) else if(pickregion < p1+p2) sample<- genregion0(-1,1,p2) else if(pickregion < p1+p2+p3) sample<- genregion1(1,2,p3) else sample<- genregion2(2,b,p4) sample } case[[7]]<- function(a,b){ p1<- c22(2,-a); p2<- k1; p3<- k0; p4<- k1; p5<- c22(2,b); pickregion<- runif(1)*(p1+p2+p3+p4+p5); if(pickregion < p1) sample<- -genregion2(2,-a,p1) else if(pickregion < p1+p2) sample<- -genregion1(1,2,p2) else if(pickregion < p1+p2+p3) sample<- genregion0(-1,1,p3) else if(pickregion < p1+p2+p3+p4)sample<- genregion1(1,2,p4) else sample<- genregion2(2,b,p5) sample } case[[6]]<- function(a,b){ genregion1(a,b,c11(a,b)) } case[[5]]<- function(a,b){ p1<- c00(a,1); p2<- c11(1,b); if(runif(1)*(p1+p2) < p1) sample<- genregion0(a,1,p1) else sample<- genregion1(1,b,p2) sample } case[[4]]<- function(a,b){ p1<- c11(1,-a); p2<- k0; p3<- c11(1,b); pickregion<- runif(1)*(p1+p2+p3); if(pickregion < p1) sample<- -genregion1(1,-a,p1) else if(pickregion < p1+p2) sample<- genregion0(-1,1,p2) else sample<- genregion1(1,b,p3) sample } case[[1]]<- function(a,b){ genregion0(a,b,c00(a,b)) } # End case list and case functions # # Main function to generate a truncated normal # rnormtrunc<- function(a,b){ flip<- 1; indexa<- sign(a)*min(2,floor(abs(a))); indexb<- sign(b)*min(2,floor(abs(b))); if(indexa+indexb < 0){t<- indexa; indexa<- -indexb; indexb<- -t; t<- a; a<- -b; b<- -t; flip<- -1} flip*case[[4*indexb+indexa+1]](a,b) } # A function to compute truncated normal density in extreme tails # dnormtrunc<- function(x,a,b){ exp(dnorm(x,log=TRUE)-pnorm(a,lower.tail=FALSE,log.p=TRUE) -log(1-exp(pnorm(b,lower.tail=FALSE,log.p=TRUE) -pnorm(a,lower.tail=FALSE,log.p=TRUE)))) } # A function to test the truncated normal generator # testt<- function(a,b,n){ x<- 1:n; i<- 1; while(i<=n){x[i]<- rnormtrunc(a,b); i<- i+1} hist(x,prob=TRUE) ulim=b; if(b==Inf) ulim<- max(x); llim=a; if(a== -Inf) llim<- min(x); xx<- seq(from=llim,to=ulim,length.out=100) lines(xx,dnormtrunc(xx,a,b)) } # A competing function using the simple inverse cdf method # rinversePhi<- function(a,b){ qnorm(pnorm(a)+runif(1)*(pnorm(b)-pnorm(a))) } # A function to test the simple inverse cdf generator # testc<- function(a,b,n){ x<- 1:n; i<- 1; while(i<=n){x[i]<- rinversePhi(a,b); i<- i+1} hist(x,prob=TRUE) ulim=b; if(b==Inf) ulim<- max(x); llim=a; if(a== -Inf) llim<- min(x); xx<- seq(from=llim,to=ulim,length.out=100) lines(xx,dnormtrunc(xx,a,b)) } # A function for full (untrucated) normal generation via rejection sampling # rnormfull<- function(n){ ptail<- exp(-2.0); pmid<- 2.0*(1.0-(1.0-exp(-0.5))/3.0); p1to2<- exp(-0.5)+exp(-2.0); phalf<- p1to2/2; sample<- 1:n; i<- 1; while(i<= n){ pickregion<- runif(1)*(pmid+p1to2+ptail); if(pickregion < pmid) sample[i]<- genmid(pmid) else if(pickregion < pmid+p1to2)sample[i]<- sign(runif(1)-0.5)*gen1to2(phalf) else sample[i]<- sign(runif(1)-0.5)*gentail() i<- i+1} sample } # Begin three generator functions for rnormfull # genmid<- function(c){ repeat{ q<- -1-(-1+runif(1)*c)/k00 C<- cp0*q x<- ym0*cos((acos(C)+4*pi)/3.0) if(runif(1)*(1.0-3.0*k00*x^2) < exp(-x^2/2)) break } x } gen1to2<- function(c){ repeat{ x<- kg1+sqrt((k11+2.0*k12)^2+4.0*k12*runif(1)*c)/(2.0*k12); if(runif(1)*(k11+2.0*k12*x) < exp(-x^2/2)) break } x } gentail<- function(){ repeat{x<- sqrt(4-2*log(runif(1))); if(runif(1)<= 2/x) break} x } # End of generators for rnormfull # # A function to test the full normal generator # testf<- function(n){ x<- rnormfull(n) hist(x,prob=TRUE) ulim<- max(x); llim<- min(x); xx<- seq(from=llim,to=ulim,length.out=100) lines(xx,dnorm(xx)) } # A function to generate truncated normal using inverse cdf robust in tails # rnorminvcdftrunc<- function(a,b){ lphibara<- pnorm(a,lower.tail=FALSE,log.p=TRUE) qnorm(lphibara+ log(1-runif(1)* (1-exp(pnorm(b,lower.tail=FALSE,log.p=TRUE)-lphibara)) ),lower.tail=FALSE,log.p=TRUE ) } # A function to test inverse cdf truncated normal generator # testtinvcdf<- function(a,b,n){ x<- 1:n; i<- 1; while(i<=n){x[i]<- rnorminvcdftrunc(a,b); i<- i+1} hist(x,prob=TRUE) ulim=b; if(b==Inf) ulim<- max(x); llim=a; if(a== -Inf) llim<- min(x); xx<- seq(from=llim,to=ulim,length.out=100) lines(xx,dnormtrunc(xx,a,b)) } # Begin functions for time tests # # for full normal via rejection sampling # timernormfull<- function(n){ ptail<- exp(-2.0); pmid<- 2.0*(1.0-(1.0-exp(-0.5))/3.0); p1to2<- exp(-0.5)+exp(-2.0); phalf<- p1to2/2; i<- 1; while(i<= n){ pickregion<- runif(1)*(pmid+p1to2+ptail); if(pickregion < pmid) sample<- genmid(pmid) else if(pickregion < pmid+p1to2)sample<- sign(runif(1)-0.5)*gen1to2(phalf) else sample<- sign(runif(1)-0.5)*gentail() i<- i+1} sample } # for basic R function rnorm # timernorm<- function(n){ i<- 1; while(i<= n){ sample<- rnorm(1) i<- i+1} sample } # for timing Weibull tail generator # timernormtail<- function(a,b,n){ i<- 1; while(i<= n){ repeat{x<- sqrt(a^2-2*log(1-runif(1)*(1-exp(-(b^2-a^2)/2)))); if(runif(1)<= a/x) break} i<- i+1} x } # for robust inverse cdf truncated normal generator # timernorminvcdf<- function(a,b,n){ i<- 1; while(i<= n){ sample<- rnorminvcdftrunc(a,b) i<- i+1} sample } # for general truncated normal via rejection sampling # timernormtrunc<- function(a,b,n){ i<- 1; while(i<= n){ sample<- rnormtrunc(a,b) i<- i+1} sample } # for general truncated normal via rejection sampling; fixed a,b # timernormtrunc0<- function(a,b,n){ flip<- 1; indexa<- sign(a)*min(2,floor(abs(a))); indexb<- sign(b)*min(2,floor(abs(b))); if(indexa+indexb < 0){t<- indexa; indexa<- -indexb; indexb<- -t; t<- a; a<- -b; b<- -t; flip<- -1} caseindex<- 4*indexb+indexa+1 i<- 1; while(i<= n){ sample<- flip*case[[caseindex]](a,b) i<- i+1} sample } # for truncated normals in region0 of rejection sampler # timeregion0<- function(a,b,n){ i<- 1; while(i<= n){ sample<- genregion0(a,b,c00(a,b)) i<- i+1} sample } # for truncated normals in region1 of rejection sampler # timeregion1<- function(a,b,n){ i<- 1; while(i<= n){ sample<- genregion1(a,b,c11(a,b)) i<- i+1} sample } # End of functions for time tests # # A function to compute percent MC error (see notes of 12/31/09) # pcmcerrbd<- function(a){ exp(log(100)-a^2/4 -0.5*(log(a)+0.5*(log(2)+log(pi)) +pnorm(a,lower.tail=FALSE,log.p=TRUE) -log(1-exp(pnorm(a,lower.tail=FALSE,log.p=TRUE) +log(a)+0.5*(log(2)+log(pi)+a^2)))) ) } # A function to compute normal tail probability # pnormtail<- function(a,mcsize){s<- 0; ssq<- 0; for(i in 1:mcsize){ y<- sqrt(a^2-2*log(runif(1))) s<- s+1/y; ssq<- ssq+1/y^2; } phibar<- s/mcsize; mcerror<- sqrt((ssq/mcsize-phibar^2)/mcsize); phibar<- exp(-0.5*a^2-0.5*log(2*pi)+log(phibar)); mcerror<- exp(-0.5*a^2-0.5*log(2*pi)+log(mcerror)); c(phibar,mcerror) } # > out<- pnormtail(cut,100000) # > array(out,dim=c(6,2)) # [,1] [,2] # [1,] 1.585300e-01 1.360114e-04 # [2,] 2.274030e-02 9.542134e-06 # [3,] 1.349548e-03 3.281790e-07 # [4,] 2.866177e-07 3.081847e-11 # [5,] 7.619595e-24 2.308249e-28 # [6,] 3.056679e-138 1.541814e-143 # A function to compute log of normal tail probability # pnormtaillog<- function(a,mcsize){s<- 0; ssq<- 0; for(i in 1:mcsize){ y<- sqrt(a^2-2*log(runif(1))) s<- s+1/y; ssq<- ssq+1/y^2; } phibar<- s/mcsize; mcerror<- sqrt((ssq/mcsize-phibar^2)/mcsize); mcerror<- mcerror/phibar; phibar<- -0.5*a^2-0.5*log(2*pi)+log(phibar); c(phibar,mcerror) }