### L14 exercise

p <- 0.4
mu <- c(-1, 2)
sd <- c(0.5, 2)
f <- function(x) {
  p * dnorm(x, mu[1], sd[1]) + (1 - p) * dnorm(x, mu[2], sd[2])
}
curve(f(x), -4, 8, n = 301, las = 1)

q1 <- function(x) rnorm(1, x, 4)
q2 <- function(x) rnorm(1, x, 50)
q3 <- function(x) rnorm(1, x, 0.1)

step <- function(x, f, q) {
  ## Pick new point
  xp <- q(x)
  ## Acceptance probability:
  alpha <- min(1, f(xp) / f(x))
  ## Accept new point with probability alpha:
  if (runif(1) < alpha)
    x <- xp
  ## Returning the point:
  x
}

run <- function(x, f, q, nsteps) {
  res <- matrix(NA, nsteps, length(x))
  for (i in seq_len(nsteps))
    res[i,] <- x <- step(x, f, q)
  return(res)
}


#1.Try 10000 samples
res <- run(-10,f,q1,10000)

plot(res, type = "l",main = "Trace plot (10000 samples, with the first 1000 burn-in)")
hist(res[1001:10000], 50, freq=FALSE, main="", ylim=c(0, .4), las=1,
     xlab="x", ylab="Probability density")
z1 <- integrate(f, -Inf, Inf)$value
curve(f(x) / z1, add=TRUE, col="red", n=200)


#2.Try different initial values 
res1 <- run(-2,f,q1,1000)
res2 <- run(10,f,q1,1000)
res3 <- run(100,f,q1,1000) # error. Log transform is needed.
plot(res1, type = "l",main = "Trace plot (1000 samples, with the first 100 burn-in)",ylim = c(-10,10))
lines(res2, col = 2 , lty = 1)


#3.Try different values of standard deviation for q
res4 <- run(-5, f, q1, 1000)
res5 <- run(-5, f, q2, 1000)
res6 <- run(-5, f, q3, 1000)

plot(res4, type = "l",ylim = c(-10,10),main = "Trace plot (1000 samples, with the first 100 burn-in)")
lines(res5, col = 2 )
lines(res6, col = 4)
legend("bottomright", c( "sd = 4 ","sd = 50","sd = 0.1"), 
       col = c(1,2,4),bty = "n",lty = c(1,1,1))

plot(density(res4),col = 2,ylim = c(0,0.7))
lines(density(res5), col = 3 , lty = 1)
lines(density(res6), col = 4,lty = 1)
curve(f(x) , add=TRUE,  n=200)
legend("topright", c("f", "sd = 4 ","sd = 50","sd = 0.1"), 
       col = c(1,2,3,4),lty = c(1,1,1,1))



new_run <- function(x, f, q, nsteps) {
  recept = 0
  res <- matrix(NA, nsteps, length(x))
  for (i in seq_len(nsteps)){
    x_new <- step(x, f, q)
    if(x_new!=x){
      recept = recept+1
    }
    res[i,] = x_new
    x = x_new
  }
  recept_rate = recept/nsteps
  return(list(res=res, recept_rate=recept_rate))
}
res7 <- new_run(-5, f, q1, 1000)
res8 <- new_run(-5, f, q2, 1000)
res9 <- new_run(-5, f, q3, 1000)
recepted_p <- c(res7$recept_rate,res8$recept_rate,res9$recept_rate)
type <- c("sd = 4","sd = 50","sd = 0.1")
data <- data.frame(type,recepted_p)
data



#4.Based on 3, plot the autocorrelation plots of the three traces.
par(mfrow = c(1,3))
acf(res4,main = "sd = 4")
acf(res5,main = "sd = 50")
acf(res6,main = "sd = 0.1")
par(mfrow = c(1,1))