Example R programs and commands
	  38.  Metropolis and  Metropolis-Hastings sampling.



##############################################################
# From Seefeld and Linder, p.203  --- Metropolis sampling,
#   target posterior beta(5,3) from proposal unif([0,1])
nchain<-1000;
theta<-vector(length=nchain); theta[1]<-0.04;
target <- function(th) {
   th^(5-1) * (1-th)^(3-1) } # unnormalized beta(5,3)
for(i in 2:nchain){
 thetaold <- theta[i-1];
 thetanew <- runif(1);  # sample the (symmetric) proposal pdf 
 r <- target(thetanew)/target(thetaold);
 u <- runif(1); 
 ifelse( u<r, theta[i]<-thetanew, theta[i]<-thetaold );
}

# Plot the Markov chains
par(mfrow=c(3,1))
plot(1:100,theta[1:100],main="First 100 Runs")
lines(1:100,theta[1:100])
plot(1:1000,theta[1:1000],main="All Runs")
lines(1:1000,theta[1:1000])
plot(901:1000,theta[901:1000],main="Last 100 Runs")
lines(901:1000,theta[901:1000])

# Compare the chain's histogram with target beta(5,3) pdf:
dev.off();
hist(theta[101:1000],nclass=25,prob=T)
xx <- (1:100)/100
lines(xx,dbeta(xx,5,3))


##############################################################
# Modified Seefeld and Linder, p.203 --- Metropolis-Hastings sampling,
#		   target posterior beta(7,12) from proposal beta(2,5)
nchain<-1000;
# Define the proposal pdf and target (unnormalized) pdf
proposal <- function(th) { dbeta(th, 2,5) }
target <- function(th) {
   th^(7-1) * (1-th)^(12-1) } # unnormalized beta(7,12)
theta<-vector(length=nchain);  # put retained samples in this array
theta[1]<-0.04;                # the initial value is not important.
for(i in 2:nchain){
 thetaold <- theta[i-1];    # previous sample
 thetanew <- rbeta(1,2,5);  # new sample from the proposal pdf.
 # proposal pdf is non-symmetric, so include the ratio
 r <- ( target(thetanew)*proposal(thetaold) )/
        ( target(thetaold)*proposal(thetanew) ) ;
 u<-runif(1);
 ifelse( u<r, theta[i]<-thetanew, theta[i]<-thetaold );
}

# Plot the Markov chains
par(mfrow=c(3,1))
plot(1:100,theta[1:100],main="First 100 Runs")
lines(1:100,theta[1:100])
plot(1:1000,theta[1:1000],main="All Runs")
lines(1:1000,theta[1:1000])
plot(901:1000,theta[901:1000],main="Last 100 Runs")
lines(901:1000,theta[901:1000])

# Compare the stationary density with beta(7,12) pdf:
dev.off();
hist(theta[51:nchain],breaks=25,prob=T);
xx<-seq(0,1,by=0.01); lines(xx,dbeta(xx,7,12));