### R commands for section 6.2.3
### G. Givens & J. Hoeting, Computational Statistics, 2nd edition, Wiley, 2013.

## example: use squeezed rejection sampling to sample N(0,1)
# target function: f(x)=exp(-x^2/2)/sqrt(2*pi)
# envelope function: e(x)=exp(-|x|)*sqrt(e/(2*pi))
# squeezing function: s(x)= 1/sqrt(2*pi*e) if -1<x<1; 0 otherwise
fg <- function(x) {            # f(x)
  exp(-x^2/2)/sqrt(2*pi);
}
fe <- function(x) {            # e(x)
  exp(-abs(x))*sqrt(exp(1)/(2*pi));
}
fs <- function(x) {            # s(x)
  ifelse( (x>-1)&(x<1), 1/sqrt(2*pi*exp(1)), 0);
}
par(mfrow=c(1,1))              # verify the relationship among f(x), e(x), s(x)
plot(fe, -3, 3, lty=2, xlab="x", ylab="")
plot(fg, -3, 3, lty=1, add=T)
plot(fs, -3, 3, lty=3, add=T)
legend(-3, fe(0), legend=c("e(x)","f(x)","s(x)"), lty=c(2,1,3))
# squeezed rejection sampling
Ginv <- function(y) {          # function to sample e(x) (cf. example for section 6.2.2 in 6_2.R)
  ifelse( y<0.5, log(2*y), -log(2*(1-y)) );
}
fexample6231 <- function(n, LIMIT=Inf) { 
  # input: sample size "n", upper limit "LIMIT" (in case acceptance ratio is too low)
  # output: random sample "y", acceptance ratios "acratio.sx" (below s(x)) and "acratio" (below f(x)) 
  y <- rep(NA,n);
  i <- 0;                      # index of y
  j <- 0;                      # index of g
  i1 <- 0;                     # index of acratio.sx
  while((i<n)&(j<LIMIT)) {
    g <- Ginv(runif(1));
    u <- runif(1);
    if( (g>-1)&(g<1) ) {
      if( u < exp(abs(g)-1) ) { i<-i+1; y[i]<-g; i1<-i1+1; } else {  # that is, if u < s(x)/e(x)
        if( u < fg(g)/fe(g) ) { i<-i+1; y[i]<-g; } 
      }
    } else if( u < fg(g)/fe(g) ) { i<-i+1; y[i]<-g; } 
    j <- j + 1;
  }
  if(i<n) cat("\n No enough random numbers! ", round(i/n*100), "% complete.\n");
  list(y=y, acratio.sx=i1/j, acratio=i/j);
}
temp <- fexample6231(100)    
temp$acratio
temp$acratio.sx
temp <- fexample6231(100, LIMIT=200)    
temp$acratio                            
temp$acratio.sx
# plot sampling results
temp <- fexample6231(100000)    
temp$acratio                   # 0.2910564
temp$acratio.sx
# theoretical acceptance ratios
integrate(fg, lower=-Inf,upper=Inf)$value/integrate(fe, lower=-Inf,upper=Inf)$value # 0.7601735
integrate(fs, lower=-Inf,upper=Inf)$value/integrate(fe, lower=-Inf,upper=Inf)$value # 0.3678794
# plot histogram
par(mfrow=c(1,1))
hist(temp$y, xlim=c(-3,3), freq=F, nclass=50, xlab=expression(lambda), 
     main="Sampling Standard Normal Distribution")
curve(fg, -3, 3, add=T)
# plot Q-Q plot
qqnorm(temp$y)
abline(0,1)