# Random Walk Metropolis Algorithm
# Scott M. Lynch, Introduction to Applied Bayesian Statistics
# and Estimation for Social Scientists, 2007, p.16, p.21, p.115f
#
print(sprintf("-----------------------------------------------------------------------"))
## [1] "-----------------------------------------------------------------------"
print(sprintf("| R-code (modified and expanded by) PCM, 2016/02/18 |"))
## [1] "| R-code (modified and expanded by) PCM, 2016/02/18 |"
print(sprintf("-----------------------------------------------------------------------"))
## [1] "-----------------------------------------------------------------------"
linear_density_graph <- function(r=0, s=5, m=2, b=3) {
# "Linear Density", Lynch, 2007, ch2.2, p.16f
c <- 2/((s-r)*(m*(s+r)+2*b))
curve(c*(m*x+b), from=r, to=s, type="l", xlab="x", ylab="f(x)", lwd=3,
main="Linear Density f(x)=(1/40)*(2*x+3)", col="steelblue")
lines(c(r,r),c(0,c*(m*r+b)),type="l",col="steelblue",lwd=1)
lines(c(s,s),c(0,c*(m*s+b)),type="l",col="steelblue",lwd=1)
legend(x=2,y=0.15, legend="area under density f(x) is 1 ")
}
linear_density_graph()
library(lattice)
## Warning: package 'lattice' was built under R version 3.2.3
m_to_b <- function (m, s=5, r=0) {(2-m*(s-r)^2) / (2*(s-r))} # m -> b
likelihood_ratio <- function() {
log_likelihood_ratio <- sum(log(z[]*m[i]+b[i])) - sum(log(z[]*m[i-1]+b[i-1])) # log of ratio
return(exp(log_likelihood_ratio))}
pr_of_accept <- function() {min(1,likelihood_ratio())} # probability of acceptance
# like flip in Church = samples from Bernoulli distribution
flip <- function(theta)
{if (rbinom(n=1,size=1,prob=theta)) TRUE else FALSE}
T <- 100000 # number of samples
T_BIn <- T*0.10 # burn-in period is 10% of T
N <- 1337 # number of subjects
r <- 0 # minimal x limit
s <- 5 # maximal x limit
thin <- 25 # thin factor (only each thin_th sample will be retained)
# initials
m=matrix(0,T); b=matrix(.2,T); z=matrix(0,1377); acc=matrix(1,T); pr_of_acceptance=matrix(0.5,T)
# data
z[1:17]=0; z[18:32]=1; z[33:103]=2 # data table 2.1, p.21
z[104:425]=3; z[426:826]=4; z[827:N]=5 # data table 2.1, p.21
hist(z,breaks=seq(r-0.5, s+0.5, 1), main="Importance of Expressing Unpopular Views", col="steelblue",
xlab="less important ... most important")
# MH-Algorithm
for(i in 2:T) { # Lynch's value is sd=0.002
m[i] <- rnorm(1,mean = m[i-1],sd = 0.001) # proposal q(x'|x) for m
b[i] <- m_to_b(m[i]) # dependent proposal for b
pr_of_acceptance[i] <- pr_of_accept()
if(! flip(pr_of_acceptance[i]) || b[i] < 0)
{ m[i] <- m[i - 1]; b[i] <- b[i - 1];
acc[i] <- 0; pr_of_acceptance[i] <- pr_of_acceptance[i-1] } # refuse proposal
# else
# { m[i] <- m[i]; b[i] <- b[i]; acc[i] <- 1; pr_of_acceptance[i] <- pr_of_acc() } # accept proposal
}
#
print(sprintf("Percentage(AcceptanceDecisions) = %f (Burn-In included)", mean(acc)*100.0))
## [1] "Percentage(AcceptanceDecisions) = 76.549000 (Burn-In included)"
print(sprintf("AcceptanceRate = %f (Burn-In included)", mean(pr_of_acceptance)))
## [1] "AcceptanceRate = 0.871779 (Burn-In included)"
plot(m, type="l", main="MH-Trace for m (Burn-In included)", col="steelblue")
lines(x=c(T_BIn,T_BIn),y=c(0,max(m)),lty=4,col="red",lwd=2) # vertical red line terminating burn-in period
legend(x=200,y=0.02, legend="Proposal Function q(x->x') = q(x'|x):
m' <- rnorm(1,mean = m[i-1],sd = 0.001)
")
acf(m, main="ACF(m) (Burn-In included)")
densityplot(m, main="Posterior Distribution for m (Burn-In included)", xlim=c(0, 0.1), col="steelblue")
# excluding samples from burn-in period
m <- m[-1: -T_BIn]
b <- b[-1: -T_BIn]
acc <- acc[-1: -T_BIn]
pr_of_acceptance <- pr_of_acceptance[-1: -T_BIn]
T <- length(m)
print(sprintf("Percentage(AcceptanceDecisions) = %f (Burn-In included)", mean(acc)*100.0))
## [1] "Percentage(AcceptanceDecisions) = 76.515556 (Burn-In included)"
print(sprintf("AcceptanceRate = %f (Burn-In included)", mean(pr_of_acceptance)))
## [1] "AcceptanceRate = 0.871724 (Burn-In included)"
plot(m, type="l", main="MH-Trace for m (Burn-In excluded)", col="steelblue")
acf(m, main="ACF(m) (Burn-In excluded)")
densityplot(m, main="Posterior Distribution for m (Burn-In excluded)", xlim=c(0, 0.1), col="steelblue")
m_est <- mean(m)
sdm_est <- sd(m)
print(sprintf("mean(m) = %f (Burn-In excluded)", m_est))
## [1] "mean(m) = 0.069693 (Burn-In excluded)"
print(sprintf("sd(m) = %f (Burn-In excluded)", sdm_est))
## [1] "sd(m) = 0.001307 (Burn-In excluded)"
plot(b, type="l", main="MH-Trace for b (Burn-In excluded)", col="steelblue")
acf(b, main="ACF(b) (Burn-In excluded)")
densityplot(b, main="Posterior Distribution for b (Burn-In excluded)", xlim=c(0, 0.1), col="steelblue")
b_est <- mean(b)
sdb_est <- sd(b)
print(sprintf("mean(b) = %f (Burn-In excluded)", b_est))
## [1] "mean(b) = 0.025767 (Burn-In excluded)"
print(sprintf("sd(b) = %f2 (Burn-In excluded)", sdb_est))
## [1] "sd(b) = 0.0032682 (Burn-In excluded)"
c_est <- 2/((s-r)*(m_est*(s+r)+2*b_est))
print(sprintf("const c = %f (Burn-In excluded)", c_est))
## [1] "const c = 1.000000 (Burn-In excluded)"
x <- c(0:5)
freq_est <- (m_est*x + b_est)*N
hist(z,breaks=seq(-0.5, 5.5, 1), main="MH-based Frequency Estimates (Burn-In excluded)", col="steelblue",
xlab="less important ... most important")
legend(x=-0.5,y=500,legend=sprintf("Percentage(AcceptanceDecisions) = %f", mean(acc)*100.0))
legend(x=-0.5,y=430,legend=sprintf("AcceptanceRate = %f",mean(pr_of_acceptance)))
mtext("Proposal Function q(x->x'): m' <- rnorm(1,mean = m[i-1],sd = 0.001)", side=3)
lines(x, freq_est,lwd=3, col="red")
lines(c(r,r),c(0,N*c_est*(m_est*r + b_est)),type="l",col="red",lwd=1)
lines(c(s,s),c(0,N*c_est*(m_est*s + b_est)),type="l",col="red",lwd=1)
# excluding samples by thinning
T <- length(m)
m <- m[seq(from=1,to=T,by=thin)]
b <- b[seq(from=1,to=T,by=thin)]
acc <- acc[seq(from=1,to=T,by=thin)]
pr_of_acceptance <- pr_of_acceptance[seq(from=1,to=T,by=thin)]
print(sprintf("Percentage(AcceptanceDecisions) = %f (Burn-In included)", mean(acc)*100.0))
## [1] "Percentage(AcceptanceDecisions) = 76.944444 (Burn-In included)"
print(sprintf("AcceptanceRate = %f (Burn-In included)", mean(pr_of_acceptance)))
## [1] "AcceptanceRate = 0.867402 (Burn-In included)"
plot(m, type="l", main="MH-Trace for m (Thinned)", col="steelblue")
acf(m, main="ACF(m) (Thinned)")
densityplot(m, main="Posterior Distribution for m (Thinned)", xlim=c(0, 0.1), col="steelblue")
m_est <- mean(m)
sdm_est <- sd(m)
print(sprintf("mean(m) = %f (Thinned)", m_est))
## [1] "mean(m) = 0.069686 (Thinned)"
print(sprintf("sd(m) = %f (Thinned)", sdm_est))
## [1] "sd(m) = 0.001310 (Thinned)"
plot(b, type="l", main="MH-Trace for b (Thinned)", col="steelblue")
acf(b, main="ACF(b) (Thinned)")
densityplot(b, main="Posterior Distribution for b (Thinned)", xlim=c(0, 0.1), col="steelblue")
b_est <- mean(b)
sdb_est <- sd(b)
print(sprintf("mean(b) = %f (Thinned)", b_est))
## [1] "mean(b) = 0.025785 (Thinned)"
print(sprintf("sd(b) = %f2 (Thinned)", sdb_est))
## [1] "sd(b) = 0.0032762 (Thinned)"
c_est <- 2/((s-r)*(m_est*(s+r)+2*b_est))
print(sprintf("const c = %f (Thinned)", c_est))
## [1] "const c = 1.000000 (Thinned)"
x <- c(0:5)
freq_est <- (m_est*x + b_est)*N
hist(z,breaks=seq(-0.5, 5.5, 1), main="MH-based Frequency Estimates (Thinned)", col="steelblue",
xlab="less important ... most important")
legend(x=-0.5,y=500,legend=sprintf("Percentage(AcceptanceDecisions) = %f",mean(acc)*100.0))
legend(x=-0.5,y=430,legend=sprintf("AcceptanceRate = %f",mean(pr_of_acceptance)))
mtext("Proposal Function q(x->x'): m' <- rnorm(1,mean = m[i-1],sd = 0.001)", side=3)
lines(x, freq_est,lwd=3, col="red")
lines(c(r,r),c(0,N*c_est*(m_est*r + b_est)),type="l",col="red",lwd=1)
lines(c(s,s),c(0,N*c_est*(m_est*s + b_est)),type="l",col="red",lwd=1)
#