Class Log to Show MCMC Steps
=============================     12/4/17

> BlockExmp = function(x,y, N) {
     XYout = array(0, c(N,2))
	 tmp = c(x,y)
	 for (j in 1:N) {
	    XYout[j,1] = rexp(1,1+4*tmp[2])
            XYout[j,2] = rexp(1,1+4*XYout[j,1])
                                  ### these last lines were previously [in class] written
                                ### incorrectly as XYout[j,] = XYout[,j]/(1+4*tmp[2:1])
         tmp = XYout[j,] }
	 XYout  }
	 
### The previous way this was written produces a Markov Chain iteration that is different
###  from the definition of Gibbs Sampler. But the output (ie the equilibrium distribution) 
###  looks the same either way.

> XY = BlockExmp(1,1,5000)
> XY = BlockExmp(XY[5000,1],XY[5000,2], 5000)

> const = 1/integrate(function(x) exp(-x)/(1+4*x))$val
> const
[1] 2.983103

> curve(const*exp(-x)/(1+4*x), 0, 2, col="blue")
  lines(density(XY[1:1000,1]), lty=2)
  lines(density(XY[1001:2000,1]), lty=4)
  lines(density(XY[2001:3000,1]), lty=6)
  lines(density(XY[3001:4000,1]), lty=8)
  lines(density(XY[4001:5000,1]), lty=9)
  
  
> curve(log(const*exp(-x)/(1+4*x)), 0, 2, col="blue")
  aux = density(XY[1:1000,1])
  lines(aux$x,log(aux$y), lty=2)
  aux = density(XY[1001:2000,1])
  lines(aux$x,log(aux$y), lty=4)
  aux = density(XY[2001:3000,1])
  lines(aux$x,log(aux$y), lty=6)
  aux = density(XY[3001:4000,1])
  lines(aux$x,log(aux$y), lty=8)
  aux = density(XY[4001:5000,1])
  lines(aux$x,log(aux$y), lty=9)

                                        ### the lack of agreement at the left side is an artifact 
                                  ### of the density estimator, not the Gibbs Sampler.
> hist(XY[1:5000,1][XY[1:5000,1]<=2], breaks=seq(0,2,length=41), prob=T)
  
## Recall 
> kerfcn = function(xnew, bw, pts,  kern=function(x) dnorm(x)) {
       n = length(pts)
       c(outer(xnew, pts, function(x,y) kern((x-y)/bw)/bw) %*% rep(1/n,n)) }
  
  curve(kerfcn(x, 0.1, XY[1:5000,1]), add=T, lty=3, col="red")
  curve(kerfcn(x, 0.031, XY[1:5000,1]), add=T, lty=3, col="blue")
  curve(kerfcn(x, 0.005, XY[1:5000,1]), add=T, lty=3, col="orange")
  
  ### YOU CAN SEE THAT THE SMALLER BANDWIDTHS REMOVE THE PROBLEM
  ### THAT IS WHY KERNEL-DENSITY ESTIMATES ON A (LEFT- OR RIGHT-) BOUNDED
  ### INTERVAL ARE OFTEN GIVEN ONE- OR TWO-SIDED MODIFICATIONS NEAR THE ENDS.
 
####################################################### 

NucPmp = cbind(Fail=c(5,1,5,14,3,19,1,1,4,22), 
   Tim=c(94.3, 15.7, 62.9, 125.8, 5.2, 31.4, 1.1, 1.0, 2.1, 10.5))
   
> t(NucPmp)
     [,1] [,2] [,3]  [,4] [,5] [,6] [,7] [,8] [,9] [,10]
Fail  5.0  1.0  5.0  14.0  3.0 19.0  1.0    1  4.0  22.0
Tim  94.3 15.7 62.9 125.8  5.2 31.4  1.1    1  2.1  10.5

> p0vec = c(rep(.1,10), .05)

> BlockEx2 = function(lamvec, bet, N) {
     parvec = c(lamvec,bet)
     Parray = array(0, c(N,11))
	 for(j in 1:N) {
	 parvec[1:10] = rgamma(10, NucPmp[,1]+1.8, NucPmp[,2]+parvec[11])
	 parvec[11] = rgamma(1,18.01, 1+sum(parvec[1:10]))
	 Parray[j,] = parvec }
	 Parray }
	 
> tmp = BlockEx2(p0vec[1:10],p0vec[11],5000)
  tmp = BlockEx2(tmp[5000,1:10],tmp[5000,11],5000)
  plot(density(tmp[,11])
> tmp = BlockEx2(tmp[5000,1:10],tmp[5000,11],5000)
  tmp = BlockEx2(tmp[5000,1:10],tmp[5000,11],5000)
  lines(density(tmp[,11]), col="blue")