# Metroplosis using Log-Normal to estimate postieor for mean
#  of normal (sigma known)
# edpsy 590BAY
#  c.j.anderson
#
#  Use as source file for metroLogNorm function
#
# Input:
#   niter   = number of iterations
#   y       = data
#   current = starting value for mean
#   tau     = standarid deviation of mean
#   jump    = the standard deviation of proposal distribution
#
# Output:
#   o x = estimates of posterior distribution
#            column 1:  current value
#            column 2:  acceptance of proposal where 1=accepted, 0=rejected
#   o plot of iteration by x (to see random walk)
#   histrogram of x
#
#
######################################################################################
# Start Function
######################################################################################
# If want to pause before next iteration active following line
#
#par(ask=F)

metroLogNorm <- function(niter, y, current, tau, jump) {

               x <- matrix(c(current,0),nrow=1,ncol=2)       # save the results of random walk 
               iteration <- matrix(1,nrow=1,ncol=1)            # to plot accepted values on y-axis }
               n <- length(y)
               
               mu <- mean(y)

               for (i in 2:niter) {
                
                propose <- rnorm(1,mean=x[i-1,1],sd=jump)
				 
                ln.ratio <-  dnorm(propose, mu, tau, log=TRUE) - dnorm(current, mu, tau, log=TRUE) 

      	        u <- log( runif(1,0,1) )
                if (u <= ln.ratio) {
                     current <- propose
                     accept <- 1 
                } else {
                     accept <- 0
                }
            
                 x <- rbind(x,c(current,accept))

                 iteration <- c(iteration,i)          

                par(mfrow=c(2,1))
                plot(iteration,x[,1],type="l",xlim=c(0,niter),ylim=c(-5, 5),main="The random walk")
                hist(x[,1],main="histogram of walk",xlim=c(-5,5),breaks=20)
                }
            return(x)
	}
#########################################################################################
# End function
########################################################################################