### Simple Antithetic & Control-variates example for class:
===========================================================

Let Y = X1 + X2 + X3 + X4
where X1 ~ Expon(1.25), X2 ~ Unif[2,4], X3 ~ Logistic, 
X4 ~ Gamma(1.5,1)  are indep r.v.'s

Problem: to find   P(Y > 5.5) .

### Method 1. Naive:
===================

tmp = .8*rexp(1e5) + runif(1e5,2,4) + rlogis(1e5) + rgamma(1e5,1.5)
round(c(est = mean(tmp>5.5), SE = sd(tmp>5.5)/sqrt(1e5)),5)
    est      SE 
0.44667 0.00157 

### Method 2. Antithetic
========================

uvec1 = runif(1e5)
uvec2 = runif(1e5)
lgs3 = rlogis(1e5)
gam4 = rgamma(1e5,1.5)
gam4A = qgamma(1-pgamma(gam4,1.5),1.5)
tmp = -0.8*log(uvec1) + 2 + 2*uvec2 + lgs3 + gam4
tmp2 = -0.8*log(1-uvec1) + 4 - 2*uvec2 - lgs3 + gam4A

> cor(tmp > 5.5, tmp2 > 5.5)
[1]  -0.6840151

>  tmp3 = c(tmp, .8*rexp(1e5) + runif(1e5,2,4) + rlogis(1e5) + rgamma(1e5,1.5))
   round(c(est = mean(tmp3>5.5), SE = sd(tmp3>5.5)/sqrt(2e5)),5)
   est      SE
0.44517 0.00111

round(c(est2 = 0.5 *(mean(tmp>5.5)+mean(tmp2>5.5)),
   SE2 = sqrt(var(tmp>5.5)+var(tmp2>5.5)-2*0.6848248*sqrt(
                  var(tmp>5.5)*var(tmp2>5.5)))/sqrt(2e5)),5)
   est2     SE2 
0.44400 0.00088       ### OK, nice improvement in SE !!

### Method 3. Control Variates:
==============================

> cor(tmp, tmp>5.5)
[1] 0.780859            ### suggests to use Y as control variate
               ### since we know E(Y) = 0.8 + 3 + 0 + 1.5 = 5.3

   ### according to theory in Sec3NotF09.pdf, the standard error
### should be decreased by a factor of sqrt(1-.781^2) = .624
### over naive estimator (whose SE we saw to be .00157

> round(summary(lm(I(tmp>5.5) ~ I(tmp-5.3)))$coef[1,1:2],5)
  Estimate Std. Error 
   0.44426    0.00098   ### OK, note that THIS SE is based on 1e5 
          ### variates, while previous one used 2e5


### Method 4. Exact numerical-integral calculation
==================================================

> fcn1 = function(x) {
    outval = numeric(length(x))
    for(i in 1:length(x)) 
       outval[i] = integrate(function(y) exp(-1.25*y)*
           dgamma(x[i]-y,1.5), 0, x[i])$val
    1.25*outval }
> fcn1(c(1,3,5))
[1] 0.25093626 0.14690617 0.03617701
> fcn2 = function(x) {
    outval = numeric(length(x))
    for(i in 1:length(x)) 
       outval[i] = integrate(function(y) (1-plogis(x[i]-y))*0.5,
           2,4)$val
    outval }

> integrate(function(x) fcn1(x)*fcn2(5.5-x),0,Inf)$val
[1] 0.4438655          ### Exact value, to high accuracy