#Lecture 6
alpha <- 10
beta  <- 30

x <- seq(0, 1, length.out=1001)
plot(x, dbeta(x, alpha, beta))

N <- 500
u <- runif(N)
sample <- runif(N)

M <- 3
M <- dbeta((alpha-1)/(alpha+beta-2), alpha, beta)

subsample <- sample[u*M < dbeta(sample, alpha, beta)]

plot(subsample, u[u*M < dbeta(sample, alpha, beta)])

MHsample <- rep(0.5, N)
for (i in 2:N) {
  prop <- runif(1)
  accept <- min(1, dbeta(prop, alpha, beta)/dbeta(MHsample[i-1], alpha, beta))
  if (runif(1)<accept)
    MHsample[i] <- prop
  else 
    MHsample[i]<- MHsample[i-1]
}

hist(MHsample)
plot(MHsample, type="l")
print(sum(MHsample[-1]==MHsample[-N]))

RWsample <- rep(0.5, N)
for (i in 2:N) {
  prop <- rnorm(1, RWsample[i-1], 100)
  # NOTE: The algorithm is correct both with and without the line below. 
  # Try both! 
  prop <- prop - floor(prop)
  accept <- min(1, dbeta(prop, alpha, beta)/dbeta(RWsample[i-1], alpha, beta))
  if (runif(1)<accept)
    RWsample[i] <- prop
  else 
    RWsample[i]<- RWsample[i-1]
}

plot(RWsample, type="l")

#Bimodal distribution: 
f <- function(x) {0.5*dbeta(x, 5, 20)+0.5*dbeta(x, 20, 5)}
plot(x, f(x), type="l")

MHsample <- rep(0.5, N)
for (i in 2:N) {
  prop <- runif(1)
  accept <- min(1, f(prop)/f(MHsample[i-1]))
  if (runif(1)<accept)
    MHsample[i] <- prop
  else 
    MHsample[i]<- MHsample[i-1]
}

plot(MHsample, type="l")
hist(MHsample)

N <- 500
RWsample <- rep(0.5, N)
for (i in 2:N) {
  prop <- rnorm(1, RWsample[i-1], 0.01)
  # NOTE: The algorithm is correct both with and without the line below. 
  # Try both! 
  prop <- prop - floor(prop)
  accept <- min(1, f(prop)/f(RWsample[i-1]))
  if (runif(1)<accept)
    RWsample[i] <- prop
  else 
    RWsample[i]<- RWsample[i-1]
}

plot(RWsample, type="l")

f <- function(x) {prod(dnorm(cars[,2], x[1]+x[2]*cars[,1]+x[3]*cars[,1]^2, 
                             x[4]  ))}
speed <- 21
N <- 10000

result <- matrix(c(0, 2, 0, 5), N, 4, byrow=TRUE)
predict <- rep(0, N)
for (i in 2:N) {
  prop <- result[i-1,] + rnorm(4, 0, c(0.2, 0.2, 0.01, 0.2))
  prop[4] <- abs(prop[4])
  accept <- f(prop)/f(result[i-1,])
  if (runif(1) < accept) result[i,] <- prop
  else result[i,] <- result[i-1,]
  predict[i] <- rnorm(1, result[i,1] + result[i,2]*speed + 
                        result[i,3]*speed^2, result[i,4])
}

par(mfrow=c(2,2))


#for (j in 1:4) plot(result[,j], type="l")

hist(predict)

print(sum(predict>70)/length(predict))