# 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)
            }