Homework Problems 5, STAT 705, Fall 2015.

Assigned 9/21/2015,  Due Friday 10/2/2015

(1) Exercise 11.1.(b), p.431  in J. Gentle book:

Find the integral of  x^2*cos(x*y)  over  (x,y) in [0,1] x [-2,2]. You should 
find the integral using Monte Carlo Integration over the 2-dimensional domain, 
for which you should provide the R code, including an error estimate (Confidence Interval).

-------------- Solution by a couple of different methods: ----------

The integrand + 1 is strictly positive and bounded above by 2, so the Monte-Carlo 
(pure simulation, accept-reject) method integral is found as 

> tmp = cbind(runif(1e7), 4*runif(1e7)-2, 2*runif(1e7))
  aux = (tmp[,3] < (1+tmp[,1]^2*cos(tmp[,1]*tmp[,2])))
  range(1+tmp[,1]^2*cos(tmp[,1]*tmp[,2])  ### = .5855, 2.0

>  c(integral= 8*mean(aux)-4, num=num, CIwidBd=2*qnorm(.995)*0.5*8/sqrt(num))
    integral          num      CIwidBd 
8.696288e-01 6.087036e+06 8.352262e-03 

### Factor of 8 is the reciprocal of the density values for (X,Y,U)
--------------------------------------------------------------------

## Another Monte Carlo Method based on treatinng X,Y as indep uniforms

Xvec = runif(1e7); Yvec = runif(1e7,-2,2)

> 4*c(Mean = mean(Xvec^2*cos(Xvec*Yvec)), 
     HalfCIwidth = (1.96/sqrt(1e7))*sd(Xvec^2*cos(Xvec*Yvec)))
        Mean  HalfCIwidth 
0.8707709649 0.0005535841 

-----------------------------------------------------------------------------
## Analytical method simplified to \int_0^1 2*x*sin(2x) dx (analytical result = 
##   (-x*cos(2*x) + 0.5*sin(2*x))_0^1 = 0.5*sin(2) - cos(2) = 0.8707955

----------------------------------------------------------------------------
## univariate numerical integral

> integrate(function(x) 2*x*sin(2*x), 0, 1)
0.8707955 with absolute error < 9.7e-15
----------------------------------------------------------------------------

## TWO MORE METHODs TO COME: bivariate trapezoid rule via linear algebra, 
## and  "integrate" using inner integral computed via trapezoid rule

## Bivariate Trapezoid Rule

wts = c(0.5,rep(1,99),0.5)
mat = outer((0:100)/100, seq(-2,2,by=4/100), function(x,y) x^2*cos(x*y))
sum(wts * (mat %*% wts)) * 4.e-4
[1] 0.8707336

wts = c(0.5,rep(1,999),0.5)
mat = outer((0:1000)/1000, seq(-2,2,by=4/1000), function(x,y) x^2*cos(x*y))
sum(wts * (mat %*% wts)) * 4.e-6
[1] 0.8707949

### Nested "integrate"

> finner = function(xvec, npt=1000) {
            m = length(xvec)  
            out = numeric(m)          
            wts = c(0.5,rep(1,npt-1),0.5)
            for(i in 1:m) out[i] = 
              xvec[i]^2*sum(cos(xvec[i]*seq(-2,2, by=4/npt))*wts)
            out*4/npt }
               
> integrate(finner,0,1)
0.8707949 with absolute error < 9.7e-15

### Nested "integrate" used to define integrand for outer "integrate"

> finner2 = function(xvec) {
            m = length(xvec)  
            out = numeric(m)          
            for(i in 1:m) out[i] = integrate(function(y,.x)  .x^2*cos(.x*y), 
                  -2,2, .x=xvec[i])$value
            out
            }

> integrate(finner2,0,1)
0.8707955 with absolute error < 9.7e-15 
      ### but this error bound is not correct because
     ### it ignores the "integrate" errors  in finner2

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

(2) Do exercises Sim.1 and Sim.3 on pp.3-4 of the handout

http://www.math.umd.edu/~evs/s400/RNG2.pdf

after reading both that handout and http://www.math.umd.edu/~evs/s400/RNG1.pdf

In Sim.3, find both a theoretical and computational justification for your answers.
(This means, a theoretically based computed answer and also one based on simulation.

#-------- Solution of (2): -----------------------

Sim.1: Since the distribution function is  1 - exp(-2*w^5), the simplest algorithm 
is based on the transformation method W_i = (-log(1-U_i)/2)^0.2, leading to

> ftmp = function(u) (-log(1-u)/2)^0.2
> round(ftmp(c(0.7432, 0.6545, 0.2110)),4)
[1] 0.9257 0.8812 0.6527

Sim.3: (a) N ~ Binom(1000, P(U>=0.6)) = Binom(1000,0.4)
       (b) Law of Large Numbers says N/1000 is close to 0.4 with high probability.

> summary(apply(array(runif(1e6)>0.6, c(1000,1000)),1,mean))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.357   0.390   0.400   0.400   0.410   0.444 

       (c) Law of Large Numbers says with high probability, A/1000 is close to 
E(exp(-U)) = \int_0^1 e^{-u} du = 1-exp(-1) = 0.6321.

> summary(apply(array(exp(-runif(1e6)), c(1000,1000)),1,mean))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6164  0.6284  0.6322  0.6323  0.6360  0.6506 

       (d) By transformation formula with g(u) = exp(-u), g'(g^{-1}(r)) = -1/r, 
the density transformation formula yields density  1/r  on [exp(-1),1]. To check 
this via simulation:

>  aux = hist(exp(-runif(1e6)), nclass=80, prob=T, plot=F)
>  curve(1/x,exp(-1),1,add=T) ### Note single leftmost histogram bar of `wrong' height