HW 19 Solution ----------------       12/17/2015

(1) Here I did not specify whether the pairs are 
to be (approximately) independent, so I would accept 
a HW solution that did not necessarily do that.

In the following function, we do not worry about independence
if   indp=F  (the default); otherwise, we select only every 1 
simulated pair out of  "gap"  (default=100)  successive ones.


BivExpBeta = function(Nrep, indp=F, gap = 100, burnin=1000) {
    outarr = array(0, c(if(indp) burnin+1+(Nrep-1)*gap 
                        else burnin+Nrep,2) )
    N = nrow(outarr)
    outarr[,1] = rexp(N)
    outarr[1,1] = outarr[1,1]/2
    outarr[1,2] = rbeta(1,1*2*outarr[1,1],3+outarr[1,1]^2)
    for(i in 2:N) {
       outarr[i,1] = outarr[i,1]/(1+2*outarr[i-1,2])
       outarr[i,2] = rbeta(1,1+2*outarr[i,1],3+outarr[i,1]^2) }
    if(indp) outarr[burnin+1+gap*(0:(Nrep-1)),]
     else outarr[burnin+1:Nrep,]  }

> dim(tmpXY1)
[1] 10000     2
> apply(tmpXY1,2,mean)
[1] 0.6110157 0.3552449

### mean of X and Y found to be respectively around 0.611, 0.355

> hist(tmpXY1[,1], nclass=60, xlim=c(0,2), prob=T)
> lines(density(tmpXY1[,1]))

   ### There are other ways to fit/display a density curve, but 
### the fact is that the correct marginal density of X is not a 
## monotone decreasing function of its argument.

> hist(tmpXY1[,2], nclass=60, prob=T)
  curve(dbeta(x,1+2*0.61, 3+0.61^2), add=T, col="red")
### The actual marginal Y-density is shifted left with respect
### to this conditional density with the mean of X plugged in for X=x.

> par(mfrow=c(3,3))
  for(i in 2:10) 
       hist(tmpXY1[1000*(i-1)+(1:1000),2], prob=T, 
       nclass=40, xlab="Y",main=paste("Sampler Obs ",
       (i-1)*1000+1,"+",sep="")) 
  par(mfrow=c(1,1))
       
### These simultaneously plotted histograms all look essentially the same
### The overall estimate of the density can be found by 

> plot(density(tmpXY1[1:1e4,2]))
  plot(density(tmpXY1[1:1e4,2], bw=.04))
    ### second one slightly smoothed; differs from any standard density
  
par(mfrow=c(3,3))
  for(i in 2:10) 
       hist(tmpXY1[1000*(i-1)+(1:1000),1], prob=T, xlim=c(0,2.5),
       nclass=40, xlab="X",main=paste("Sampler Obs ",
       (i-1)*1000+1,"+",sep="")) 
  par(mfrow=c(1,1))

### A little hard to assess convergence, but looks OK
### Overall density:

> plot(density(tmpXY1[1:1e4,1]), xlim=c(0,3.5))
  
### Any kind of assessment of summary statistics across blocks of 
generated Gibbs-Sampler values contributes to testing convergence, e.g.

round(apply(array(tmpXY1,c(1000,10,2)),2, function(mat)
   cor(mat[,1],mat[,2])), 4)
 [1] 0.1995 0.1720 0.1605 0.2340 0.1942 0.1921 0.1714 0.1715 0.1891 0.1872

## These are the correlations between X and Y in blocks of 1000; 
## the overall correlation is 0.1864. One could also look at this 
## as a time series (with correlations taken in batches of 80)

plot(apply(array(tmpXY1,c(80,125,2)),2, function(mat)
   cor(mat[,1],mat[,2])))

### Behavior looks stationary enough, with some scatter


(2) 
set.seed(7793)
xv = runif(80)
yv = 1.5*xv + 0.8*rnorm(80)

Tsum0 = c(sum(yv^2), sum(xv*yv), sum(xv^2))
b0 = Tsum0[2]/Tsum0[3]

## For the posterior density, just do Gibbs sampler, using  
##     parameters  (b, tau) where tau=1/sig^2  and 

f(tau|b, xv,yv) ~ Gamma(41, .05 + 0.5*sum(Tsum * c(1,-2*b,b^2)) )

f( b |tau, xv,yv) ~ N( mu, 1/(tau*(.01+Tsum[3])) )

PostGibbs = function(Nrep,Tsum) {
   outarr = array(0,c(Nrep,2))
   outarr[,1] = rnorm(Nrep,Tsum[2]/Tsum[3])
   outarr[1,2] = rgamma(1,41, 0.5*sum(Tsum*c(1,
             -2*outarr[1,1],outarr[1,1]^2)))
   for(i in 2:Nrep) {
      tau=outarr[i-1,2]
      b = Tsum[2]/(.01/tau+Tsum[3]) + 
             sqrt(1/(.01+tau*Tsum[3]))*outarr[i,1]
      outarr[i,2] = rgamma(1,41,0.5*sum(Tsum*c(1,-2*b,b^2)))
      outarr[i,1]=b }
   outarr }


> tmp2 = PostGibbs(1e5,Tsum0)
  hist(tmp2[,1],nclass=50,prob=T)
  curve(dnorm(x, mean(tmp2[,1]),sd(tmp2[,1])), add=T, col="red")
  ### very accurate agreement

  round(rbind(ML=c(b=b0, sig2=sum(Tsum0*c(1,-2*b0,b0^2))/80),
      post.mean = c(mean(tmp2[,1]), mean(1/tmp2[,2]))),4)
               b   sig2
ML        1.5645 0.5461
post.mean 1.7772 0.5705

## NB true values are 1.5, 0.64

## This result suggests there are still small-sample effects

  hist(1/tmp2[,2],nclass=50,prob=T)
  curve(dnorm(x, mean(1/tmp2[,2]),sd(1/tmp2[,2])), add=T, col="red")
 ### The variance posterior density is a little skewed compared to normal

### Let's over-plot the histogram of the first 20,000 tau's (after burn-in) 
versus the last 20,000.

 hist(tmp2[1:2e4,2],nclass=50, prob=T, col="orange")
 lines(hist(tmp2[8e4+(1:2e4),2], nclass=50, plot=F), freq=F)
### Rather small differences between the two histograms