### R commands for section 6.2 (continued)
### G. Givens & J. Hoeting, Computational Statistics, Wiley, 2005.

## example for section 6.2.2
## g(x) = (e^x)/2 if x<0; (e^{-x})/2 otherwise
## G(x) = (e^x)/2 if x<0; 1-(e^{-x})/2 otherwise
## G^{-1}(y) = ln(2y) if 0<y<0.5; -ln[2(1-y)] otherwise 
g <- function(x) {             # function g(x)
  ifelse( x<0, exp(x)/2, exp(-x)/2 );
}
G <- function(x) {             # function G(x)
  ifelse( x<0, exp(x)/2, 1-exp(-x)/2 );
}
Ginv <- function(y) {          # function G^{-1}, 0<x<1
  ifelse( y<0.5, log(2*y), -log(2*(1-y)) );
}
par(mfrow=c(1,1))
plot(c(-2,2), c(-2,2), xlab="x", ylab="", main="Example for Section 6.2.2", type="n")
plot(G, -2, 2, lty=1, add=T)
plot(g, -2, 2, lty=2, lwd=2, add=T)
plot(Ginv, 0, 1, lty=3, lwd=2, add=T)
abline(0,1, lty=4, lwd=1)
legend(-2, 2, legend=c("G(x)", "g(x)", "G^{-1}", "y=x"), lty=c(1,2,3,4), lwd=c(1,2,2,1))
# function to sample from g(x)
gsample <- function(n) {       # input sample size n and output sample vector
  ans <- runif(n);
  Ginv(ans);
}
x <- gsample(100000)
hist(x, xlim=c(-3, 3), freq=F, nclass=200)
curve(g, add=T)

## Example 6.2 on page 150: Sampling a Bayesian posterior
x <- c(8,3,4,3,1,7,2,6,2,7)    # observations
mean(x)                        # sample mean estimate = 4.3
a <- sum(x)                    # sum of observations
b <- sum(log(factorial(x)))    # log(8!) + log(3!) + ... + log(7!)
nx <- length(x)                # number of observations
fLx <- function(lambda) {      # function L(lambda|x)
  exp(a*log(lambda) - b - nx*lambda);
}
flambda <- function(lambda) {  # function f(lambda)
  dlnorm(lambda,meanlog=log(4),sdlog=0.5);
}
# re-produce Figure 6.2 on page 150
fe <- function(lambda) {       # e(lambda)*10^11    
  flambda(lambda)*fLx(4.3)*10^11;
}
fq <- function(lambda) {       # q(lambda)*10^11
  flambda(lambda)*fLx(lambda)*10^11;
}
par(mfrow=c(1,1))              # for more math notations in plot, cf "plotmath"
plot(fe, 0, 20, xaxt="n", yaxt="n", xlab=expression(lambda), 
     ylab=expression("Unnormalized density  " %*% 10^11), 
     main="Fig.6.2: Target (Dotted) and Envelope (Solid) Functions")
plot(fq, 0, 20, lty=3, add=T)
axis(1, at=c(0,10,20), labels=c("0","10","20"))
axis(2, at=c(0,1,2,3), labels=c("0","1","2","3"))

# sampling function for example 6.2
L43 <- fLx(4.3)                # L(4.3|x)
fexample62 <- function(n, LIMIT=Inf) { 
  # input: sample size "n", upper limit "LIMIT" (in case acceptance ratio is too low)
  # output: random sample "y" and acceptance ratio "acratio" 
  y <- rep(NA,n);
  i <- 0;                      # index of y
  j <- 0;                      # index of lambda
  while((i<n)&(j<LIMIT)) {
    lambda <- rlnorm(1, meanlog = log(4), sdlog = 0.5);
    u <- runif(1);
    if(u < fLx(lambda)/L43) {
      i <- i + 1;
      y[i] <- lambda;
    };
    j <- j + 1;
  }
  if(i<n) cat("\n No enough random numbers! ", round(i/n*100), "% complete.\n");
  list(y=y, acratio=i/j);
}
temp <- fexample62(100)    
temp$acratio
temp <- fexample62(100, LIMIT=200)    
temp$acratio                            
# plot sampling results
temp <- fexample62(100000)    
temp$acratio                   # 0.2910564
par(mfrow=c(1,1))
hist(temp$y, xlim=c(0,10), freq=F, nclass=50, xlab=expression(lambda), 
     main="Example 6.2: Sampling a Bayesian Posterior")
cfq <- integrate(fq, 0, 20)$value
fq1 <- function(lambda) { fq(lambda)/cfq; }    # normalized target density function
curve(fq1, 0, 10, add=T)