Homework Set 11 Solutions Monday October 24, 2016.
------------------------------------------------------

IMPORTANCE SAMPLING SIMULATION PROBLEM.

This problem was adapted from the Robert and Casella (1999) book 
on MonteCarlo, and was originally taken from a journal paper.

Problem:     Estimate P(W_j > 0 for j=1,...6)
where    W = c(  (diag(0:5) + 1) %*% Z )    and 
Z is a vector of 6 iid N(0,1) random variables.

## Note the correlations:

> V = A %*% A 
  varA = diag(V)
  diag(1/sqrt(varA)) %*% V %*% diag(1/sqrt(varA))
> CorMat = diag(1/sqrt(varA)) %*% V %*% diag(1/sqrt(varA))
  round(CorMat,3)
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
[1,] 1.000 0.953 0.873 0.802 0.745 0.701
[2,] 0.953 1.000 0.802 0.727 0.669 0.625
[3,] 0.873 0.802 1.000 0.642 0.586 0.543
[4,] 0.802 0.727 0.642 1.000 0.518 0.477
[5,] 0.745 0.669 0.586 0.518 1.000 0.428
[6,] 0.701 0.625 0.543 0.477 0.428 1.000

## So the correlations are initially large, progressively smaller

Straightforward approach:

A = diag(0:5)+1
W = array(rnorm(6e6), c(1e6,6)) %*% A
p = 1 - mean( c((W<0) %*% rep(1,6)) > 0 )
[1] 0.200849
### So this is the answer, with confidence interval +-
> 1.96*sqrt(p*(1-p))/1e3
[1] 0.0007852453

### As some students remarked, this is not a rare event and importance
###    sampling is by no means necessary to estimate its probability.
##-----------------------------------------------------------------

### It is also not so easy to find a density completely supported on
## the positive orthant in R^6 from which it is easy to simulate

Let's ease up on the restriction of positive-orthant support
and use g as density of (|W_1|,W_2,..,W_6)

Wst = cbind(abs(W[1:5e5,1]),W[1:5e5,-1])
Ainv = solve(A)

X = W[1:5e5,] %*% Ainv
Xst = Wst %*% Ainv
> 2*mean(  (c((Wst[,-1]<0) %*% rep(1,5)) == 0 )/
     (1+exp(0.5*(X^2-Xst^2) %*% rep(1,6))) )
[1] 0.2028472

### confidence interval width
> 1.96*2/sqrt(5e5)*sd((c((Wst[,-1]<0) %*% rep(1,5)) == 0 )/
     (1+exp(0.5*(X^2-Xst^2) %*% rep(1,6))) )
[1] 0.0011146
### just about exactly the same as
> 1.96*sqrt(p*(1-p)/5e5)
[1] 0.0011105

#----------------------------------------------------------------
### Now for a 2nd method of importance sampling, let the density g(w) 
###   for sampling be the density of coordinatewise  abs(W)
### But simulate only 2e5 batches of 6.

W = array(rnorm(2e5), c(1e6,6)) %*% A

Ainv = solve(A)

## Now use 6-digit binary to index the coordinatewise multiplications by -1

B0 = array(rep(diag(6),64), c(6,6,64))

M = cbind(rep(c(0,1),32), rep(c(0,0,1,1),16),
     rep(c(0,0,0,0,1,1,1,1),8), rep(c(rep(0,8),rep(1,8)),4),
     rep(c(rep(0,16),rep(1,16)),2), rep(c(0,1),each=32))

for(i in 1:6) B0[i,i,M[,i]==1]= -B0[i,i,M[,i]==1]

> table(c(B0))
  -1    0    1 
 192 1920  192 

### array B0 contains 64 6x6 matrices for applying all possible sign 
###   switches to 6 positive vector entries

> Baux = t(array(Ainv %*% array(B0, c(6,6*64)), c(36,64)))
  aux0 = t(array(array(Baux, c(64*6,6)) %*% abs(t(W[1:2e5,])), c(64,6*2e5)))
  ## 614Mb, 6 x 2e5 x 64 array
  aux0 = array(exp(-0.5*(rep(1,6) %*% array(aux0^2, c(6,2e5*64)))), c(2e5,64))
  
> mean(1/c(1 + (aux0[,-1]%*%rep(1,63))/aux0[,1]))
[1] 0.2016865
> (1.96/sqrt(2e5))*sd(1/c(1 + (aux0[,-1]%*%rep(1,63))/aux0[,1]))
[1] 0.0007978531

### This method achieved the same SE as the straightforward method, 
## with only 2e5 simulated vectors. But is is not an efficient method
## in this example.

##%%%-----------------------------------------------------------

Now try the same importance-sampling idea in the case where A is replaced 
by Ainv, with negative correlations and a MUCH smaller orthant probability!

> W2 =  array(rnorm(6e6), c(1e6,6)) %*% Ainv
> p2 = 1 - mean( c((W2<0) %*% rep(1,6)) > 0 )
> p2
[1] 9e-06    ### on another run I got 1.7e-05

> Baux2 = t(array(A %*% array(B0, c(6,6*64)), c(36,64)))
  aux2 = t(array(array(Baux2, c(64*6,6)) %*% abs(t(W2[1:2e5,])), c(64,6*2e5)))
  ## 614Mb, 6 x 2e5 x 64 array
  aux2 = array(exp(-0.5*(rep(1,6) %*% array(aux2^2, c(6,2e5*64)))), c(2e5,64))
  
> mean(1/c(1 + (aux2[,-1]%*%rep(1,63))/aux2[,1]))
[1] 1.679379e-05
> (1.96/sqrt(2e5))*sd(1/c(1 + (aux2[,-1]%*%rep(1,63))/aux2[,1]))
[1] 1.066795e-06

#----------------------------------------------------------------