### R commands for section 6.3.2 -- Antithetic sampling
###            and section 6.3.3 -- Control variates
### G. Givens & J. Hoeting, Computational Statistics, Wiley, 2005.

## Example 6.6 (Normal expectation) on page 171
# need to estimate mu=E[h(X)], h(x)=x/(2^x-1), X~N(0,1)
fh <- function(x) {        # h(x)=x/(2^x-1)
  ifelse(x!=0, x/(2^x-1), 1/log(2));
}
plot(fh, -2, 2)
temp <- function(x) {      # function to calculate theoretical value
  fh(x)*dnorm(x);
}
integrate(temp, lower=-Inf, upper=Inf)$value  # theoretical value: 1.499143

n <- 100000                # simple Monte Carlo sample size
x <- rnorm(n)              # simple Monte Carlo sample
y <- fh(x)
mean(y)                    # MC estimate: 1.501694
sqrt(var(y)/n)             # estimated std of MC estimator: 0.001610868

m <- 50000                 # antithetic sample size
cor(y[1:m], fh(-x[1:m]))   # estimate of cor(h(X), h(-X)):  -0.9517503
y1 <- (y[1:m] + fh(-x[1:m]))/2 
mean(y1)                   # antithetic estimate: 1.499378
sqrt(var(y1)/m)            # estiamted std of antithetic estimator: 0.0003523504


## Example 6.8 (A control variate for importance sampling) on page 173
# need to estimate E[sqrt(x)], X~exponential(1)
temp <- function(x) { sqrt(x)*dexp(x); }
integrate(temp, lower=0, upper=Inf)$value   # theoretical value
#[1] 0.8862265

n <- 10000

## Envelope: abs(Normal(0,1))
set.seed(1)
x <- abs(rnorm(n))
wts <- dexp(x)/(dnorm(x)/.5)
tx <- sqrt(x)*wts             # t(X)=h(x)*f(X)/g(X)
mu.is <- mean(tx)             # Importance sampling estimate (weights are not standardized)
mu.is
#[1] 0.8535775
sqrt(var(tx)/n)               # std of IS estimate with non-standardized weights, eq. (6.32)
#[1] 0.01500608

mu.wts <- mean(wts)           # bar{c} in eq. (6.61)
mu.wts
#[1] 0.9862718
var.theta <- sum((wts-mu.wts)^2)/(n*(n-1))            # eq. (6.61)
cov.mutheta <- sum((tx-mu.is)*(wts-mu.wts))/(n*(n-1)) # eq. (6.62)
lambda <- -cov.mutheta/var.theta                      # eq. (6.60)   
lambda
#[1] -1.880021

mu.is + lambda*(mu.wts-1)     # eq. (6.63)
#[1] 0.8793868
sqrt(var(tx)/n - cov.mutheta^2/var.theta)     # std of Mu_ISCV in eq. (6.63)
#[1] 0.004444861