P <- matrix(c(1/3,0,1/3,0,1/3,0,0,1/3,1/3,0,1/3,0,1/3,0,1/3,0,0,1/3,
              1/3,0,0,1/3,0,1/3,0,1/3,0,1/3,1/3,0,0,1/3,0,1/3,0,1/3),nrow=6)
P
                      [,1]             [,2]            [,3]            [,4]            [,5]            [,6]
[1,] 0.3333333 0.0000000 0.3333333 0.3333333 0.0000000 0.0000000
[2,] 0.0000000 0.3333333 0.0000000 0.0000000 0.3333333 0.3333333
[3,] 0.3333333 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000
[4,] 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333
[5,] 0.3333333 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000
[6,] 0.0000000 0.0000000 0.3333333 0.3333333 0.0000000 0.3333333

MC.multiply <- function(P.in,N)  {
       S <- c(1,0,0,0,0,0)%*%P.in
       for (i in 2:N)  {
           S <- S%*%P.in
           print(S)
       }
}

MC.multiply(P,15)


B <- 5; k <- 15; m <- 500; x <- NULL; y <- NULL
while (length(x) < m)  {
    x.val <- c(runif(1,0,B),rep((B+1),length=k))
    y.val <- c(runif(1,0,B),rep((B+1),length=k))
    for (j in 2:(k+1))  {
        while(x.val[j] > B) x.val[j] <- rexp(1,y.val[j-1])
        while(y.val[j] > B) y.val[j] <- rexp(1,x.val[j])
    }
    x <- c(x,x.val[(k+1)])
    y <- c(y,y.val[(k+1)])
}
\end{VM}


coal.mining.disasters <- c(4,5,4,0,1,4,3,4,0,6,3,3,4,0,2,6,
                           3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,
                           2,2,3,4,2,1,3,2,2,1,1,1,1,3,0,0,
                           1,0,1,1,0,0,3,1,0,3,2,2,0,1,1,1,
                           0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,
                           3,3,1,1,2,1,1,1,1,2,4,2,0,0,1,4,
                           0,0,0,1,0,0,0,0,0,1,0,0,1,0,1)

num.reps <- 2000
coal.mat <- matrix(NA,ncol=3,nrow=num.reps)
coal.mat[1,] <- c(1,1,100)
alpha <- 4; beta <- 1; gamma <- 1; delta <- 2
coal.mat <- bcp(coal.mat,coal.mining.disasters,alpha,beta,gamma,delta)
summary(coal.mat[1000:2000,])

                bcp <- function(theta.matrix,y,a,b,g,d)  {
                    n <- length(y)
                    k.prob <- rep(0,length=n)
                    for (i in 2:nrow(theta.matrix))  {
                        lambda <- rgamma(1,a+sum(y[1:theta.matrix[(i-1),3]]), 
                                           b+theta.matrix[(i-1),3])
                        phi    <- rgamma(1,g+sum(y[theta.matrix[(i-1),3]:n]), 
                                           d+length(y)-theta.matrix[(i-1),3])
                        for (j in 1:n)  k.prob[j] <- exp(j*(phi-lambda))*
                                                     (lambda/phi)^sum(y[1:j])
                        k.prob <- k.prob/sum(k.prob)
                        k      <- sample(1:n,size=1,prob=k.prob)
                        theta.matrix[i,] <- c(lambda,phi,k)
                    }
                    return(theta.matrix)
                }      


metropolis <- function(num.iter,I.mat)  {
      theta.matrix <- rbind(c(3,-2),matrix(NA,nrow=num.iter,ncol=2))
      for (i in 2:nrow(theta.matrix))  {
          theta.prime <- mvrnorm(1,theta.matrix[(i-1),],I.mat)/(sqrt(rchisq(2,5)/5))
          a <-dmultinorm(theta.prime[1],theta.prime[2],c(0,0),I.mat)/
              dmultinorm(theta.matrix[(i-1),1],theta.matrix[(i-1),2], c(0,0),I.mat)
          if (a > runif(1)) theta.matrix[i,] <- theta.prime
          else theta.matrix[i,] <- theta.matrix[(i-1),]
      }
      theta.matrix
}

mcmc.mat <- metropolis(300,matrix(c(1,0.9,0.9,1),ncol=2))
par(mfrow=c(1,1),mar=c(3,2,2,2),oma=c(3,3,1,1),col.axis="white",col.lab="white",
    col.sub="white",col="white",bg="black")
plot(mcmc.mat,pch=1,cex=0.3,xlim=c(-3,3),ylim=c(-3,3))
lines(mcmc.mat,lwd=0.5,col="red")
axis(side=1,at=-3:3,NULL,col="white"); axis(side=2,at=-3:3,NULL,col="white")

dmultinorm <- function(xval,yval,mu.vector,sigma.matrix)  {
        normalizer <- (2*pi*sigma.matrix[1,1]*sigma.matrix[2,2]*sqrt(1-sigma.matrix[1,2]^2))^(-1)
        like <- exp(-(1/(2*(1-sigma.matrix[1,2]^2)))* (
             ((xval-mu.vector[1])/sigma.matrix[1,1])^2
             -2*sigma.matrix[1,2]*(((xval-mu.vector[1])/sigma.matrix[1,1])*((yval-mu.vector[2])/sigma.matrix[2,2]))
             + ((yval-mu.vector[2])/sigma.matrix[2,2])^2 ) )
        pdf.out <- normalizer * like
        return(pdf.out)
}

hit.run <- function(num.iter,I.mat)  {
    theta.mat <- rbind(c(0,0),matrix(NA,nrow=num.iter,ncol=2))
    for (i in 2:nrow(theta.mat))  {
        rad.draw <- runif(1,0,2*pi)
        if ((rad.draw > 0 & rad.draw < pi/2) | (rad.draw > pi & rad.draw < 3*pi/2)) 
                distance <- rgamma(1,1,3)
        else    distance <- rgamma(1,1,6)
        xy.theta <- c(distance*cos(rad.draw),distance*sin(rad.draw)) + theta.mat[(i-1),]
        a <- dmultinorm(xy.theta[1], xy.theta[2],              c(0,0),I.mat)/
             dmultinorm(theta.mat[(i-1),1], theta.mat[(i-1),2],c(0,0),I.mat)
        u <- runif(1)
        if (runif(1) < a) { theta.mat[i,] <- xy.theta }
        else theta.mat[i,] <- theta.mat[(i-1),]
    }
    theta.mat
}
 
mcmc.mat <- hit.run(10000,matrix(c(1,0.95,0.95,1),ncol=2))
postscript("Class.Multilevel/hit.run.contours.ps")
par(mfrow=c(1,1),mar=c(3,2,2,2),oma=c(3,3,1,1),col.axis="white",col.lab="white",
    col.sub="white",col="white",bg="black")
plot(mcmc.mat[(nrow(mcmc.mat)/2+1):nrow(mcmc.mat),],pch=1,cex=0.1,
    xlim=c(-4,4),ylim=c(-4,4))

alpha.grid <- seq(-4,4,length=300)
beta.grid <- seq(-4,4,length=300)
density <- outer(alpha.grid,beta.grid,dmultinorm,mu.vector=c(0,0),
    sigma.matrix=matrix(c(1,0.95,0.95,1),ncol=2))
contour(alpha.grid,beta.grid,density,levels=c(.2,.1,.05),col="red",add=TRUE)