######################################################
#M-H Demo
#2/2/06
######################################################

setwd("U:\\Teaching\\Stat694\\Computing\\")

######################################################
#M-H Demo -- unimodal distribution
######################################################

#Function for calculating the acceptance probability
calculate.alpha.unimodal <- function(value.new,value.old,true.mean,true.sd){
                       alpha <- exp(-.5*((value.new-true.mean)/true.sd)^2)/exp(-.5*((value.old-true.mean)/true.sd)^2)
                       return(min(alpha,1))
                       }
                      
#Set parameters of the true posterior distribution of theta
true.mean <- 5
true.sd <- 2

#Set the standard deviation of the proposal distribution
proposal.sd <- .8

#Set up the demonstration
K <- 20   #K is the number of samples (MH-steps)
k.hist <- 50  #start drawing the histograms
n.modes <- 1

x <- seq(-5,15,length=500)   #Grid of thetas to evaluate the density function                 
theta.old <- 1   #ititialize the parameter     
theta.samples <- rep(0,K) #theta storage vector
accept <- rep(0,K)        #acceptance storage

#Run the demo
source(file="mh-demo-plots.r")

######################################################
#M-H Demo -- biimodal distribution
######################################################

calculate.alpha.bimodal <- function(value.new,value.old,p.1,p.2,true.mean.1,true.sd.1,true.mean.2,true.sd.2){
                       alpha <- (p.1*exp(-.5*((value.new-true.mean.1)/true.sd.1)^2)+p.2*exp(-.5*((value.new-true.mean.2)/true.sd.2)^2))/(p.1*exp(-.5*((value.old-true.mean.1)/true.sd.1)^2)+p.2*exp(-.5*((value.old-true.mean.2)/true.sd.2)^2))
                       return(min(alpha,1))
                       }

#Set parameters of the true posterior distribution of theta
true.mean.1 <- 7
true.sd.1 <- 2
true.mean.2 <- 1
true.sd.2 <- 1

#Set mixture probabilies
p.1 <- .5
p.2 <- 1-p.1

#Set the standard deviation of the proposal distribution
proposal.sd <- 1.5

#Set up the demonstration
K <- 500   #K is the number of samples (MH-steps)
k.hist <- 50  #start drawing the histograms
n.modes <- 2

x <- seq(-5,15,length=500)   #Grid of thetas to evaluate the density function      
theta.old <- 1   #ititialize the parameter     
theta.samples <- rep(0,K) #theta storage vector
accept <- rep(0,K)        #acceptance storage     

#Run the demo
source(file="mh-demo-plots.r")