# Problem 1:

tot.draws <- 40000 # number of sampled values from each full conditional

# Initial values for quantities of interest:
delta.init <- 2
theta.init <- 0

mu.delta <- c(-3,0,3)

# Dummy Vectors that will hold values for samples of interest:
delta.vec <- c(delta.init, rep(NULL, times = tot.draws) )
theta.vec <- c(theta.init, rep(NULL, times = tot.draws) )

for (j in 2:(tot.draws+1) ) {
 theta.vec[j] <- rnorm(n=1, mu.delta[delta.vec[j-1]], sd=sqrt(1/3))
 my.prob.vec.numerators <- c(.45*dnorm(theta.vec[j], -3, sd=sqrt(1/3)), 
                             .10*dnorm(theta.vec[j], 0, sd=sqrt(1/3)),
                             .45*dnorm(theta.vec[j], 3, sd=sqrt(1/3)) )
 my.prob.vec <- my.prob.vec.numerators / sum(my.prob.vec.numerators)
 delta.vec[j] <- sample(1:3, size=1, prob=my.prob.vec)
}

theta.post <- theta.vec[-1]

## Posterior summary for theta:

hist(theta.post,freq=F)  

my.seq <- seq(-6,6,length=1000)

lines(my.seq, 0.45*dnorm(my.seq,-3, 1/3) + 0.10*dnorm(my.seq,0, 1/3) + 0.45*dnorm(my.seq,3, 1/3))

# Redo with num.draws = 40000

Part(e):

#trace plot:
plot(theta.post,type='l')

# ACF plot:
acf(theta.post)

#########################################################################

# Problem 2:

x<-c(0.92,0.42,3.62,0.89,-0.69,0.45,-0.11,-0.14,-0.47,1.09,-0.34,0.62,0.27)
y<-c(0.26,1.65,2.10,0.62,-1.16,1.29,-0.82,-0.36,-0.29,0.86,0.19,1.25,0.33)

sum.xsq <- sum(x^2)
sum.xy <- sum(x*y)
sum.ysq <- sum(y^2)
n <- length(x)

S <- 30000
rho.current <- 0.5 # initial value for M-H algorithm
acs <- 0 # will be to track "acceptance rate"

rho.values <- rep(0,times=S)  # will store sampled values of rho

for (s in 1:S) {
 rho.proposed <- runif(1,min=rho.current-0.2, max=rho.current+0.2)
 if (rho.proposed < 0) rho.proposed <- abs(rho.proposed)
 if (rho.proposed > 1) rho.proposed <- (2 - rho.proposed)

log.accept.ratio <-  (-0.5*n*log(1-rho.proposed^2) - (1/(2*(1-rho.proposed^2))*(sum.xsq - 2*rho.proposed*sum.xy + sum.ysq))) - (-0.5*n*log(1-rho.current^2) - (1/(2*(1-rho.current^2))*(sum.xsq - 2*rho.current*sum.xy + sum.ysq)))

  if (log.accept.ratio > log(runif(1)) ) {
    rho.current <- rho.proposed
    acs <- acs + 1
  }

rho.values[s] <- rho.current
}

acs/S  # gives the acceptance rate

acf(rho.values)

plot(rho.values,type='l')


# Thinning out the sampled values by taking every 10th rho:

rho.values.thin <- rho.values[10*(1:(S/10) )]

acf(rho.values.thin)

plot(rho.values.thin,type='l')

### Posterior summary:

plot(density(rho.values.thin))

quantile(rho.values.thin, probs=c(0.025,0.5,0.975))

library(TeachingDemos)
emp.hpd(rho.values.thin)

#########################################################################



## Problem 7.6 from book:

#(a)
set.seed(84735)
rnorm(n=1,mean=4.6,sd=2)

#(b)
set.seed(84735)
rnorm(n=1,mean=2.1,sd=7)

#(c)
set.seed(84735)
runif(n=1,min=8.9-2,max=8.9+2)

#(d)
set.seed(84735)
runif(n=1,min=1.2-0.5,max=1.2+0.5)

#(e)
set.seed(84735)
runif(n=1,min=7.7-3,max=7.7+3)


## Problem 7.7 from book:

proposed<-2.1
current<-2

#(a)

top<-((proposed)^(-2))*dnorm(current,mean=proposed, sd=1)
bottom<-((current)^(-2))*dnorm(proposed,mean=current, sd=1)
accept.ratio <- top/bottom
print(accept.ratio)

#(b)

top<-(exp(proposed))*dnorm(current,mean=proposed, sd=0.2)
bottom<-(exp(current))*dnorm(proposed,mean=current, sd=0.2)
accept.ratio <- top/bottom
print(accept.ratio)

#(c)

top<-(exp(-10*proposed))*dunif(current,min=proposed-0.5, max=proposed+0.5)
bottom<-(exp(-10*current))*dunif(proposed,min=current-0.5, max=current+0.5)
accept.ratio <- top/bottom
print(accept.ratio)

#(d)

top<-(exp(-(proposed^4)))*dexp(current,rate=proposed)
bottom<-(exp(-(current^4)))*dexp(proposed,rate=current)
accept.ratio <- top/bottom
print(accept.ratio)