#################### Assignment 6.1 ###################
####### Assignment 6.1.(a) Obtain Mixture Data #######
# simulate the mixture data
set.seed(50)
options(digits=15)
N = 500
# sample N random uniform U
U =runif(N)
# variable to store the samples from the mixture distribution
rand.samples = rep(NA,N)
# sampling from the mixture
for(i in 1:N){
if(U[i]<0.7){
rand.samples[i] = rnorm(1,7,0.5)
}else{
rand.samples[i] = rnorm(1,10,0.5)
}
}
mixture.data <- as.data.frame(rand.samples)
y = mixture.data$rand.samples
# draw histogram and plot of mixture distribution
x = seq(5,15,by=0.01);delta = 0.7
d = delta * dnorm(x,7,0.5) + (1-delta) * dnorm(x,10,0.5)
hist(y,breaks=40,freq=FALSE,main="Histogram of Mixture Data",ylab="Density")
points(x,d,type="l")
plot(x,d,type="l")
###### Assignment 6.1.(b) Independence chain MCMC procedure ######
# functions
f = function(delta,data){prod(delta*dnorm(data,7,0.5)+(1-delta)*dnorm(data,10,0.5))}
MH.ratio = function(delta.t,delta,data){(f(delta,data)*g(delta.t))/(f(delta.t,data)*g(delta))}
# initial Values
set.seed(50)
niter = 10500
delta.set = NULL
# main chain
g = function(delta){dbeta(delta,1,1)}
delta.set[1] = rbeta(1,1,1)
for(i in 1:niter){
delta.t = delta.set[i]
delta.p = rbeta(1,1,1)
p = min(MH.ratio(delta.t,delta.p,y),1)
u = rbinom(1,1,p)
delta.set[i+1] = delta.p * u + delta.t * (1-u)
}
mean(delta.set[501:(niter+1)])
## [1] 0.728231851612344
sd(delta.set[501:niter])
## [1] 0.0194620114865262
# par(mfrow=c(1,1))
plot(delta.set[501:(niter+1)],ylim=c(0,1),type="l",xlab=expression(delta),ylab="t")
title("Sample Path for Beta(1,1) Proposal Distribution")
hist(delta.set[501:(niter+1)],breaks=40,xlab=expression(delta),
main="Histogram for Beta(1,1) Proposal Distribution")
acf(delta.set)
###### Assignment 6.1.(c) Random walk chain ######
# establish Initial Values
set.seed(50)
niter = 10500
delta = rep(0, niter); delta[1] = runif(1,0,1);# starting location for random walk
# functions
f = function(delta,data){prod(delta*dnorm(data,7,0.5)+(1-delta)*dnorm(data,10,0.5))}
# main Function
for(i in 1:(niter-1)){
temp = delta[i] + runif(1,-1,1)
while (temp<0 | temp>1){
temp = delta[i] + runif(1,-1,1)
}
prop = temp
MH.ratio = f(prop,y)/f(delta[i],y)
u <- runif(1,0,1)
if(u < MH.ratio){
delta[i+1] = prop
}else{
delta[i+1] = delta[i]
}
}
mean(delta[501:niter])
## [1] 0.72665715157548
sd(delta[501:niter])
## [1] 0.0206981970312185
plot(delta[501:niter],ylim=c(0,1),type="l",ylab=expression(delta),xlab="t")
hist(delta[500:niter],breaks=40,xlab=expression(delta),main="Hist. for Random Walk")
acf(delta)
###### Assignment 6.1.(d) Reparameterize ######
# establish Initial Values
set.seed(50)
niter = 10500
u1 = rep(0, niter); u1[1] = runif(1,-1,1);
p1 = rep(0, niter); p1[1] = exp(u1[1])/(1+exp(u1[1]));
# functions
loglike = function(delta,data){sum(log(delta*dnorm(data,7,0.5)+(1-delta)*dnorm(data,10,0.5)))}
# main Function
for(i in 1:(niter-1)){
u1[i+1] = u1[i] + runif(1,-1,1)
p1[i+1] = exp(u1[i+1])/(1+exp(u1[i+1]))
MHratio = exp(loglike(p1[i+1],y)-loglike(p1[i],y))*exp(u1[i])/exp(u1[i+1])
if(MHratio < 1)
if(rbinom(1,1,MHratio)==0){p1[i+1] = p1[i]; u1[i+1] = u1[i]}
}
mean(p1[501:niter])
## [1] 0.72534316735217
sd(p1[501:niter])
## [1] 0.0202784100320009
plot(p1[501:niter],ylim=c(0,1),type="l",ylab=expression(delta),xlab="t")
hist(p1[501:niter],breaks=40,xlab=expression(delta),main="Hist. for Unif(-1,1) Walk")
acf(p1)
############### Assignment 6.2 #################
############# Assignment 6.2.(a) ############
Metropolis <- function(F_sample # distribution to sample
, F_prop # proposal distribution
, y0 # starting value
, I = 10000 # iterations
){
y = rep(NA,T)
y[1] = y0 # starting location for random walk
accepted = c(1)
for(t in 2:I) {
y.prop <- F_prop(y[t-1]) # random walk proposal
R = F_sample(y.prop)/F_sample(y[t-1])
u <- runif(1,0,1) # uniform variable to determine acceptance
if(u < R){ # accept the new value
y[t] = y.prop
accepted = c(accepted,1)
}
else{
y[t] = y[t-1] # reject the new value
accepted = c(accepted,0)
}
}
return(list(y, accepted))
}
set.seed(50)
y0 <- c(0,7,15)
test = function(x){0.7*dnorm(x,7,0.5)+0.3*dnorm(x,10,0.5)}
# random walk, sigma = 0.01
# y0 = y0[1]
response1 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, 0.01)}
,y0 = y0[1]
,I = 10000
)
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
## [1] 0.9478
mean(y_trace1)
## [1] 3.24538446465904
plot(density(y_trace1))
plot(y_trace1,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace1,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
# y0 = y0[2]
response2 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, 0.01)}
,y0 = y0[2]
,I = 10000
)
y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # acceptance rate without considering burn-in
## [1] 0.9944
mean(y_trace2)
## [1] 7.18066846768784
plot(density(y_trace2))
plot(y_trace2,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace2,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
# y0 = y0[3]
response3 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, 0.01)}
,y0 = y0[3]
,I = 10000
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # acceptance rate without considering burn-in
## [1] 0.9655
mean(y_trace3)
## [1] 12.3650049156097
plot(density(y_trace3))
plot(y_trace3,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace3,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
########### Assignment 6.2.(b) ###########
# sigma = 1
response1 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x,1)}
,y0 = y0[1]
,I = 10000
)
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
## [1] 0.5322
mean(y_trace1)
## [1] 7.9544925737704
plot(density(y_trace1))
plot(y_trace1,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace1,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
# y0 = y0[2]
response2 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, 1)}
,y0 = y0[2]
,I = 10000
)
y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # acceptance rate without considering burn-in
## [1] 0.5304
mean(y_trace2)
## [1] 7.86929380558813
plot(density(y_trace2))
plot(y_trace2,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace2,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
# y0 = y0[3]
response3 <- Metropolis(F_sample = test
,F_prop = function(x){rnorm(1, x, 1)}
,y0 = y0[3]
,I = 10000
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # acceptance rate without considering burn-in
## [1] 0.5276
mean(y_trace3)
## [1] 7.86246774396089
plot(density(y_trace3))
plot(y_trace3,ylim=c(0,15),type="l",xlab="x",ylab="t")
title("Sample Path for Normal Proposal Distribution")
par(mfrow=c(1,1))
s = seq(5,15,by=0.01)
hist(y_trace3,breaks=50,freq=FALSE,xlab="x",main="Histogram of Mixture Data",ylab="Density")
points(s,test(s),type="l")
abline(v=7);abline(v=10) # Highlight the approximate modes of the true distribution
############### Assignment 6.3 #################
#initialize constants and parameters
N <- 5000 # length of chain
burn <- 1000 # burn-in length
rho <- c(0.3,0.5,0.7) #correlation
# initialization
x <- rep(0, N); y <- rep(0, N)
x[1]=1; y[1]=1
# generate the chain
gibbs <- function(rho){
for (i in 2:N) {
# sample the x direction conditional on y
x[i] <- rnorm(1, rho*y[i-1], sqrt(1-rho^2))
# sample the y direction conditional on x
y[i] <- rnorm(1, rho*x[i], sqrt(1-rho^2))
}
return(cbind(x,y))
}
b <- burn + 1 # burn-in
# rho = 0.3
X1 <- gibbs(rho[1]) # the chain, a bivariate sample
x1 <- X1[b:N, ]
# compare sample statistics to parameters
colMeans(x1)
## x y
## -0.00887840782043451 -0.01808631928779589
cov(x1)
## x y
## x 1.025797394049380 0.298339687433479
## y 0.298339687433479 1.003749177053653
cor(x1)
## x y
## x 1.000000000000000 0.294013744738008
## y 0.294013744738008 1.000000000000000
plot(x1, main="Bivariate Normal", xlab="X", ylab="Y", ylim=range(x1[,2]))
# rho = 0.5
X2 <- gibbs(rho[2])
x2 <- X2[b:N, ]
# compare sample statistics to parameters
colMeans(x2)
## x y
## 0.0251289893083738 0.0158887088165823
cov(x2)
## x y
## x 0.990304091280521 0.475366083033126
## y 0.475366083033126 0.961768477885859
cor(x2)
## x y
## x 1.000000000000000 0.487089352081737
## y 0.487089352081737 1.000000000000000
plot(x2, main="Bivariate Normal", xlab="X", ylab="Y", ylim=range(x2[,2]))
# rho = 0.7
X3 <- gibbs(rho[3])
x3 <- X3[b:N, ]
# compare sample statistics to parameters
colMeans(x3)
## x y
## 0.0318876897951430 0.0268764744128826
cov(x3)
## x y
## x 0.975297251476826 0.693843302139826
## y 0.693843302139826 0.996856672543793
cor(x3)
## x y
## x 1.00000000000000 0.70368215877767
## y 0.70368215877767 1.00000000000000
plot(x3, main="Bivariate Normal", xlab="X", ylab="Y", ylim=range(x3[,2]))
############### Assignment 6.4 #################
#initialize constants and parameters
N <- 5000 # length of chain
burn <- 1000 # burn-in length
# function
f <- function(x){dnorm(x,-3,1)}
# initialization
x1 <- rep(0, N); x2 <- rep(0, N)
x1[1]=runif(1,0,1); x2[1]=runif(1,0,f(x1[1]))
for (i in 2:N) {
# sample the x1 direction conditional on x2
x1[i] <- runif(1, 0, -3+sqrt(-2*log(x2[i-1])))
# sample the x2 direction conditional on x1
x2[i] <- runif(1, 0, f(x1[i]))
}
b <- burn + 1 # burn-in
X <- cbind(x1,x2) # the chain, a bivariate sample
x <- X[b:N, ]
# compare sample statistics to parameters
colMeans(x)
## x1 x2
## 0.496668989000769501 0.000731115676966559
cov(x)
## x1 x2
## x1 0.134358386307548539 -1.88984226856934e-04
## x2 -0.000188984226856934 6.53518812667498e-07
cor(x)
## x1 x2
## x1 1.000000000000000 -0.637770051444662
## x2 -0.637770051444662 1.000000000000000
plot(x, main="Bivariate Normal", xlab="X", ylab="Y", ylim=range(x[,2]))
s1 = seq(-8,2,by=0.01)
s2 = seq(0,2,by=0.01)
plot(s1,f(s1),type="l")
plot(s2,f(s2),type="l")
