# Edpsy 590BAY
# fall 2021
# c.j.anderson
#
#    This tunes metroLogNorm for estimating posterior
#      p(mu|y, sigma) where y is normal .  This produces 
#      a number of acceptance ratios for different values 
#      jumping distribution
#

# input
#   metroLogNorm as source
#   n.iter for metroLogNorm
#   y data
#   min.jump sd
#   max.jump sd
#   n.jumps    number of jumps between min and max
#
# What are acceptance rates
#

setwd("D:/Dropbox/edps 590BAY/Lectures/9 Metropolis Hastings")

anorex <- read.table("anorexia_data.txt",header=TRUE)

# create change 
anorex$change <- anorex$weight2 - anorex$weight1

#############################################
# Sample statistics
#############################################
n=nrow(anorex)
ybar = mean(anorex$change)
s2 = var(anorex$change)
stderr <- sqrt(s2/n)
sample.stat <- c(n,ybar,s2,stderr)
names(sample.stat) <- c("n","ybar","s2","stderr")
sample.stat

# Function with graphing removed
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)          

                }
             return(x)
	}
#########################################################################################
# End function
########################################################################################

par(mfrow=c(2,2))

# For tunings and/or metroLogNorm
nchains <- 20
y <- anorex$change
min.jump <- .01
max.jump <- 10
njumps <- 50

# for metroLogNorm
ybar <- mean(y)
stderr <- sd(y)/sqrt(length(y))
n.iter <-1000

# for tunning
acceptance <- matrix(NA, nrow=njumps, ncol=2)
jumps <- seq(from=min.jump, to=max.jump, length.out=njumps)

# Start of loop through jump standard deviations
i <- 1
for (jump.sd in jumps) {
   chain <- metroLogNorm(niter=n.iter, y, ybar, stderr, jump.sd)

   jsd <- round(jump.sd, digits=4)
   plot(seq(1:n.iter), chain[,1], type="l", main=paste("Jump sd =", jsd))

   accept <- table(chain[,2])
   percent <- accept[2]/n.iter * 100
   acceptance[i, 1] <- percent
   acceptance[i, 2] <- jump.sd
   i <- i + 1
}   

# Graph them
par(mfrow=c(1,1))
plot(acceptance[, 2], acceptance[, 1], type='l', lwd=2, 
     xlab="Jump sd",
	 ylab="Percent Accept",
	 main="Acceptance Rate for Anorexia")
lines(c(min.jump,max.jump),c(50,50), lty=2, col="cyan")
lines(c(min.jump,max.jump),c(20,20), lty=2, col="cyan")