# Example R code for the Metropolis and Metropolis-Hastings algorithms


##############################################################
# From Seefeld and Linder, p.203  --- Metropolis sampling, beta(a,b)
k<-3; n<-5;       # experiment results
a1<-k; b1<-(n-k); # posterior shape alpha-1,beta-1, after experiment
nsteps<-1000;   # number of steps in the Markov chain
X<-vector(length=nsteps); # allocate storage for the steps
x<-0.04; X[1]<-x;               # initial state, somewhat arbitrary
gx<-x**a1 * (1-x)**b1;          # initial posterior value g(x)
for(i in 2:nsteps){         # run the Markov process
 y<-runif(1);       # proposal: sample the (symmetric) prior
 gy<-y**a1 * (1-y)**b1; # posterior value g(y), at proposal
 r<-gy/gx;        # ratio of posteriors
 u<-runif(1);     # acceptance function: compare r with U(1)
 if(u<r) {
    X[i]<-y;  # accept the new posterior (g) value
    gx<-gy; x<-y;
 } else {
    X[i]<-x; # retain the old posterior (g) value
    # ...and retain "gx" and "x"
 }
}
# Plot the start, middle, and end of the Markov chain
par(mfrow=c(3,1))
plot(1:100,X[1:100],main="First 100 Runs")
lines(1:100,X[1:100])
plot(1:1000,X[1:1000],main="All Runs")
lines(1:1000,X[1:1000])
plot(901:1000,X[901:1000],main="Last 100 Runs")
lines(901:1000,X[901:1000])
# Compare the histogram with the beta(a,b) pdf:
dev.off();
hist(X[101:1000],nclass=25,prob=T)
xx <- (1:100)/100
lines(xx,dbeta(xx,a1+1,b1+1))


##############################################################
# Modified Seefeld and Linder, p.203 --- Metropolis-Hastings sampling,
#                  posterior beta(7,12) from prior beta(2,5)
nsteps<-1000;
k<-5; n<-12;                        # binomial experiment variables
alpha0<-2; beta0<-5;                # prior shape parameters
alpha<-alpha0+k; beta<-beta0+(n-k); # posterior shape parameters
a1<-alpha-1; b1<-beta-1;

X<-vector(length=nsteps);       # allocate
x<-0.04;  X[1]<-x;              # ...and initialize
gx<-x**a1 * (1-x)**b1;          # initial posterior value g(x)
qyx<-dbeta(x,alpha0,beta0);     # (for all y)q(y,x)=proposal(x)
for(i in 2:nsteps){
 y<-rbeta(1,alpha0,beta0); # proposal: sample the non-symmetric prior 
 gy<-y**a1 * (1-y)**b1;    # posterior value at proposal                
 qxy<-dbeta(y,alpha0,beta0);    # (for all x)q(x,y)=proposal(y)
 r<-(gy*qyx)/(gx*qxy);          # corrected ratio of posteriors
 ##### NOTE:
 # If g is a posterior pdf from Bayes' rule g(y)=experiment(y)*prior(y),
 # and the proposal density is the prior, then the corrected ratio of
 # posteriors simplifies to r=experiment(y)/experiment(x).
 #####
 u<-runif(1);                   # used for acceptance test
 if(u<r) {
    X[i]<-y;  # accept the new posterior (g) value
    gx<-gy; x<-y; qyx<-qxy;
 } else {
    X[i]<-x; # retain the old posterior (g) value
    # ...and retain "gx" and "x" and "qyx"
 }
}
# Plot the start, middle, and end of the Markov chain
par(mfrow=c(3,1))
plot(1:100,X[1:100],main="First 100 Runs")
lines(1:100,X[1:100])
plot(1:1000,X[1:1000],main="All Runs")
lines(1:1000,X[1:1000])
plot(901:1000,X[901:1000],main="Last 100 Runs")
lines(901:1000,X[901:1000])
# Compare the histogram with the beta(7,12) pdf:
dev.off();
hist(X[51:nsteps],breaks=25,prob=T);
xx<-(1:100)/100; lines(xx,dbeta(xx,alpha,beta));