HW21 solution

Bmat = array( c(1.2,1.4,-0.6,2.3,-0.6,0.8,-1,0,1, 2.9,-.7,-.7,1.1,.2,2,1,0,-1,
         -.3,.7,-.5,1.9,-.5,-1,0,1,0, -.4,2.8,.3,,1.5,2.4,-.3,0,-1,0), c(9,4))

> Bmat
      [,1] [,2] [,3] [,4]
 [1,]  1.2  2.9 -0.3 -0.4
 [2,]  1.4 -0.7  0.7  2.8
 [3,] -0.6 -0.7 -0.5  0.3
 [4,]  2.3  1.1  1.9  1.5
 [5,] -0.6  0.2 -0.5  2.4
 [6,]  0.8  2.0 -1.0 -0.3
 [7,] -1.0  1.0  0.0  0.0
 [8,]  0.0  0.0  1.0 -1.0
 [9,]  1.0 -1.0  0.0  0.0

### From the 7th and 9th constraints, we learn that abs(x[1]-x[2]) <=1 .
### From the 8th:   x[3] - x[4] >= -1
### Then using these and the 5th:  1.9 x[4] -.4 x[1] >= -1.7
### Further checking of the same sort implies the set D is bounded.

### Metropolis-Hastings: start with 

x = rep(0,4)

## and update with random-walk steps in random direction of fixed length
## so that the  q(x,y) proposal kernel is symmetric, and check only if y in D

> StepChk = function(x.in, lgth=0.3) {
    vec = rnorm(4)
    vec=vec/sqrt(sum(vec^2))
    xtmp = x.in + lgth*vec
    if(min( 1+c(Bmat%*%xtmp)) > 0) xtmp else x.in  }

> Xarr = array(0, c(1.5e4,4))
  for(i in 1:1.5e4) {
         x = StepChk(x)
         Xarr[i,] = x }

> apply(Xarr[5001:1.5e4,],2,range)
           [,1]       [,2]      [,3]       [,4]
[1,] -0.7718544 -0.5235437 -1.329350 -0.4468657
[2,]  1.5887388  1.4386245  2.074533  2.7776675

### compare hists of Xarr coords on 5001:1e4 and (1e4+1):1.5e4 blocks
### to check convergence

> par(mfrow=c(2,4))
> for(j in 1:2) for(i in 1:4) 
      hist(Xarr[5000*j+(1:5000),i], nclass=30, prob=T,
      xlab=paste("Coord",i), main = "")

### Not quite the same in two rows, so do several further blocks of 1e4.

> for(j in 1:10) 
       for(i in 1:1e4) {
         x = StepChk(x)
         Xarr[i,] = x }

  for(j in 1:2) for(i in 1:4) 
      hist(Xarr[5000*(j-1)+(1:5000),i], nclass=30, prob=T,
      xlab=paste("Coord",i), main = "")

## Maybe still not converged
> sum(diff(Xarr[,1])==0)
[1] 5059

### So maybe not enough of the new steps are accepted. Try smaller

> for(j in 1:10) 
       for(i in 1:1e4) {
         x = StepChk(x, lgth=.15)
         Xarr[i,] = x }

 for(j in 1:2) for(i in 1:4) 
      hist(Xarr[5000*(j-1)+(1:5000),i], nclass=30, prob=T,
      xlab=paste("Coord",i), main = "")

> sum(diff(Xarr[,1])==0)
[1] 3348

### Maybe we'll get a better result with a much larger array and small steps

Xarr2 = array(0, c(1e6,4))

> system.time({  for(j in 1:2) 
       for(i in 1:1e6) {
         x = StepChk(x, .1)
         Xarr2[i,] = x } })
   user  system elapsed 
    8.5     0.0     8.5 

par(mfrow=c(2,2))
for(j in 1:2) for(i in 1:2) 
      hist(Xarr2[5e5*(j-1)+10*(1:5e4),i], nclass=40, prob=T,
      xlab=paste("Coord",i), main = "")
for(j in 1:2) for(i in 3:4) 
      hist(Xarr2[5e5*(j-1)+10*(1:5e4),i], nclass=40, prob=T,
      xlab=paste("Coord",i), main = "")

### Seems roughly OK, still not perfectly the same; 
###  note that we are using every 10th coord in array of 1e6
### Now answer further questions:


## (ii)  Accept-Reject algorithm.
## Start by checking a box containing all the points
> apply(Xarr2,2,range)
           [,1]       [,2]      [,3]       [,4]
[1,] -0.8856825 -0.5744754 -1.443349 -0.4591163
[2,]  1.6940125  1.5486182  2.175326  3.1061695

## Use box [-1,1.8] x [-.7,1.7] x [-1.55,2.25] x [-.6,3.2]

> Xarr = array(0, c(1e5,4))
  ctr=0; j=0
  system.time({  while( j < 1e7 & ctr < 1e5) {
      j=j+1
      xtmp = runif(4, c(-1,-.7,-1.55,-.6), c(1.8,1.7,2.25,3.2))
      if( min(1+c(Bmat %*% xtmp)) > 0) {
             Xarr[ctr,] = xtmp
             ctr=ctr+1 } } })

   user  system elapsed 
   6.26    0.02    6.29 

### This makes it seem as though the two methods are comparable

> for(i in 1:4) {
          densA = density(Xarr2[10*(1:1e5),i], bw=.05)
          densB = density(Xarr[1:1e5,i], bw=.05)
      plot(densA, main=paste0("Coord ",i,", two densities"), 
          xlab=paste("Coord",i))
      lines(densB, col="red", lty=2)   }
 
###   P(X[1] > 0 and X[2] > 0 | X[3] > 0 and X[4] > 0)

> inds = (Xarr[,3]>0 & Xarr[,4]>0)
  mean( Xarr[inds,1]>0 & Xarr[inds,2] > 0 )
[1] 0.6406837

## numerator is ~ N(.40258 , .40258*(1-.40258)/1e5)
## while denominator is ~ N(0.62836 , 0.62836*(1-0.62836)/1e5)

> mean(Xarr[,3]>0 & Xarr[,1]>0 & Xarr[,2]>0  & Xarr[,4]>0 )
[1] 0.40258
> mean(Xarr[,3]>0 & Xarr[,4]>0 )
[1] 0.62836

### So the Delta Method gives the ratio ~ N( 0.6406837, 
###   (0.6406837^2/1e5) *(.40258*(1-.40258)/.40258^2 + 0.62836*(1-0.62836)/.62836^2 ) )

## SE = 
>  (0.6406837/sqrt(1e5))*sqrt((1-.40258)/.40258 + (1-0.62836)/.62836)
[1] 0.00291875
### and CI = + or - 1.96*0.00291875 =
> 0.6406837 + c(-1,1)*1.96*0.00291875
[1] 0.6349629 0.6464045

###==========================================================================
### SUPPLEMENT ON OTHER METHODS OF METROPOLIS-HASTINGS GENERATION OF X arrays
##---------------------------------------------------------------------------

## Try to generate uniform in largest sphere within D centered at xtmp

Bnrm = sqrt(apply(Bmat^2,1,sum))

qker = function(x,y) {
       rmax = min((1+c(Bmat %*% x))/Bnrm)
       rnew = sqrt(sum((x-y)^2))
       if(rnew<rmax) 4*rnew^3/rmax^4 else 0 }

MHstep = function(xin) {
       rmax = min((1+c(Bmat %*% xin))/Bnrm)
       utmp = rnorm(4)
       xout = xin + (rmax*(runif(1)^(.25))/sqrt(sum(utmp^2)))*utmp
       rho = qker(xout,xin)/qker(xin,xout)
       if(runif(1)<rho) xout else xin }

Xarr3 = array(0,c(1e5,4))
xtmp = rep(0,4)
for(i in 1:1e5) {
       xtmp = MHstep(xtmp)
       Xarr3[i,] = xtmp  }

system.time({ for(j in 1:5) for(i in 1:1e5) {
       xtmp = MHstep(xtmp)
       Xarr3[i,] = xtmp  } })
   user  system elapsed 
   5.98    0.00    6.01 

### Might already have converged! No, not yet

> for(i in 1:4) {
          densA = density(Xarr2[10*(1:1e5),i], bw=.05)
          densB = density(Xarr[1:1e5,i], bw=.05)
          densC = density(Xarr3[1:1e5,i], bw=.05)
      plot(densA, main=paste0("Coord ",i,", three densities"), 
          xlab=paste("Coord",i))
      lines(densB, col="red", lty=2)
      lines(densC, col="blue", lty=4)   }