gibbs.binorm <- function( N = 1000, means = c(0,0), corr = 0.5, start = means )
  {
    s <- sign( corr )
    ## generate correlated bivariate normal
    n <- 1
    x1 <- rnorm( n )
    x2 <- rnorm( n )
    z <- rnorm( n, 0, corr / (s - corr ) )
    y <- c( means[1] + x1 + z, means[2] + s * ( x2 + z))
    y <- means
    sd <- sqrt( 1 - corr^2 )

    ## Gibbs sampler
    gibbs <- matrix(NA,N+1,2)
    gibbs[1,] <- start
    for( i in seq(N)) {
      j <- sample(1:2,1)
      jj <- 3 - j
      gibbs[i+1,j] <- rnorm(1,y[j] + corr * ( gibbs[i,jj] - y[jj] ), sd )
      gibbs[i+1,jj] <- rnorm(1,y[jj] + corr * ( gibbs[i+1,j] - y[j] ), sd )
    }
    gibbs
  }
ellipse <- function(x=0,y=0,sx=1,sy=1,rho=0,fine=.005,center=0,show=T,add=T,
                       xlab="",ylab="",bty="n",xaxt="n",yaxt="n",type="n",lty=1,pty="s",err=-1,
                       ...) {
     if(min(sx) < 0 | min(sy) < 0 | max(abs(rho)) > 1)
             stop("sx,sy or rho out of bounds")
        if(fine <= 0 | fine > 1)
                stop("fine out of bounds")
        n <- max(length(x),length(y),length(sx),length(sy),length(rho))
        x <- array(x,n)
        y <- array(y,n)
        lty <- array(lty,n)
        sx <- array(sx,n)
        sy <- array(sy,n)
        rho <- array(rho,n)
        center <- array(center,n)
        prin1 <- sqrt((1+rho)/2)
        prin2 <- sqrt((1-rho)/2)
        xarc <- sin(pi * seq(0,.5,by=4*fine))
        yarc <- c(rev(xarc),-xarc)
        xarc <- c(xarc,yarc,-rev(xarc))
        yarc <- c(yarc,rev(yarc))
        arc <- list()
        for (i in 1:n)
                arc[[i]] <- data.frame(x=x[i]+sx[i]*(prin1[i]*xarc+prin2[i]*yarc),
                                                y=y[i]+sy[i]*(prin1[i]*xarc-prin2[i]*yarc))

        if(show) {
                if(!add)
                  {
                              tmp <- options( warn = -1 )
                                        plot(c(min(x-sx),max(x+sx)),c(min(y-sy),max(y+sy)),
                                                            bty=bty,xaxt=xaxt,yaxt=yaxt,type=type,xlab=xlab,ylab=ylab,pty=pty,
                                                            ...)
                                        options( tmp )
                            }
                      for (i in 1:n) {
                                 lines(arc[[i]]$x,arc[[i]]$y,lty=lty[i],err=err,...)
                                          if(center[i] > 0) {
                                                        lines(x[i]+center[i]*c(-1,1)*sx[i]*prin1[i],
                                                                             y[i]+center[i]*c(-1,1)*sy[i]*prin1[i],err=err)
                                                                    lines(x[i]+center[i]*c(-1,1)*sx[i]*prin2[i],
                                                                                         y[i]+center[i]*c(1,-1)*sy[i]*prin2[i],err=err)
                                                      }
                               }
              }
        invisible(arc)
   }


gibbs.plot <- function(corr=.75,N=200,cex=1,means=c(0,0),...)
  {
    par(mfcol=c(2,2),mar=c(3.1,3.1,0.1,0.1))
    tmp <- gibbs.binorm(corr=corr,N=N,means=means,...)
    for(i in 1:2) {
      plot(tmp[,i],type="b",xlab="",ylab="",cex=cex,lwd=2)
      abline( h = means[i], col = "blue", lty = 2, lwd = 2 )
      mtext("Markov chain index",1,2)
      mtext(paste("Gibbs: mean",i),2,2)
    }
    plot(tmp,type="b",xlab="",ylab="",cex=cex,lwd=2)
    abline( 0, corr, lwd = 2, col = "blue", lty = 2)
    for(i in 1:2)
      mtext(paste("Gibbs: mean",i),i,2)
    plot(tmp,xlab="",ylab="",cex=cex,lwd=2)
    for(i in 1:2)
      mtext(paste("Gibbs: mean",i),i,2)
    for(i in 1:2)
      ellipse(rho=corr,add=T,center=0,col="blue",x=means[1],y=means[2],sx=i,sy=i,lwd=3)
    invisible(tmp)
  }

#Metropolis-Hastings: local
mhlocal <- function( l = 10000,bw=5,jump=.001,plotit=(l<100000),rejectit=TRUE, cex = 0.8,
	x = seq(0,10,len=201),
	y = .75*dnorm(x,4,.625)+.25*dnorm(x,7,.625))
{
	z <- sample( seq(-bw,bw), l, replace = TRUE )
	znow <- sample(seq(201),1)
	z[1] <- x[znow]
	ynow <- y[znow]

	if( plotit ) {
		plot( range(x), c(0,l),type="n",xlab="",ylab="")
		mtext(expression(theta),1,2)
		mtext("mcmc sequence",2,2)
#		lines( x, y * l / max( y ), col = "blue" )
		points(z[1],1,cex=cex)
	}
	jumpit <- FALSE
	rejects <- 0
	for( i in 2:l ) {
		znew <- znow + z[i]
		if( znew < 1 )
			znew <- 2 - znew
		if( znew > 201 )
			znew <- 402 - znew
		if( jump > 0 ) {
			jumpit <- runif(1) < jump 
			if( jumpit )
				znew <- sample( seq( 201 ), 1 )
		}

		ynew <- y[ znew] 
		if( ynew >= ynow )
			znow <- znew
		else {
			if( runif(1) <= ynew / ynow )
				znow <- znew
		}
		z[i] <- x[znow]
		ynow <- y[znow]
		if( z[i] > 10 | z[i] < 0 ) browser()
		rejects <- rejects + ( znow != znew )
		if( plotit ) {
			if( rejectit & znow != znew )
				points(x[znew],i,col="red",cex=cex)
			if( jumpit ) {
				points(z[i],i,col="blue",cex=cex)
				abline( h = i, col = "black" )
			}
			else
				points(z[i],i,col="blue",cex=cex)
		}
		jumpit <- FALSE
	}
	cat( "rejection rate =", round( 100 * rejects / l ), "\n" )
	hist( z, nc = 201, prob = TRUE, main="",xlab="",ylab="")
	mtext( expression( theta ), 1, 2 )
	mtext( expression(paste("pr(",theta,"|",italic(Y),")")), 2, 2 )
	lines(x,y,lwd=4,col="blue")
	invisible(z)
}