N <- 100
RWsample <- rep(0.5, N)
for (i in 2:N) {
  prop <- rnorm(1, RWsample[i-1], 0.1) 
  prop <- abs(prop)     # This "reflects" the proposal function at 0, but not at 1. 
  accept <- (prop >= 0) & (prop <= 1)   #Discuss why the acceptance prob becomes this! 
  if (runif(1) < accept) RWsample[i] <- prop
  else RWsample[i] <- RWsample[i-1]
}
#plot(RWsample, type="l")
#hist(RWsample)


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

result <- matrix(c(0, 2, 0, 5), N, 4, byrow=T)
predict <- rep(0, N)
count <- 0
for (i in 2:N) {
  prop <- result[i-1,] + rnorm(4, 0, c(0.5, 0.2, 0.01, 0.2))
  accept <- min(1, f(prop)/f(result[i-1,]))
  if (runif(1) < accept) {
    count <- count + 1
    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])
}
print("Acceptance rate:")
print(count/N)
print("Answer to question:")
print(sum(predict>70)/length(predict))
par(mfrow=c(2,2))
#for (j in 1:4) plot(result[,j], type="l")
for (j in 1:4) acf(result[,j])


MAT <- cov(result)
print(MAT)
library(mgcv)

result <- matrix(c(0, 2, 0, 5), N, 4, byrow=T)
predict <- rep(0, N)
# ADD LATER: 
 count   <- 0
for (i in 2:N) {
#  prop <- result[i-1,] + rnorm(4, 0, c(0.2, 0.2, 0.01, 0.2))
   prop <- result[i-1,] + rmvn(1, rep(0, 4), 0.1*MAT)
    accept <- min(1, f(prop)/f(result[i-1,]))
  if (runif(1) < accept) {
    count <- count + 1
    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])
}
print("Acceptance rate:")
print(count/N)
print("Answer to question:")
print(sum(predict>70)/length(predict))
par(mfrow=c(2,2))


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

##for (j in 1:4) acf(result[,j])
##for (j in 1:4) acf(result[,j])

N <- 1000
data <- c(4, 6, 7, 5, 5)
ndata <- length(data)

#data <- data - mean(data)

result <- matrix(c(mean(data), sd(data)), N, 2, byrow=TRUE)
for (i in 2:N) {
  result[i,1] <- rnorm(1, (ndata*result[i-1,2]*mean(data))/
                         (1+ndata*result[i-1,2]), 
                       1/sqrt(1+ndata*result[i-1,2]))
  result[i,2] <- rgamma(1, ndata/2+3, 0.5*sum((data-result[i,1])^2)+4)
}

par(mfrow=c(1,1))
plot(result[,1], sqrt(1/result[,2]))


###################
# Heart transplants! 
###################

# install.packages("LearnBayes")
library(LearnBayes)
data("hearttransplants")
attach(hearttransplants)
e = hearttransplants$e
y = hearttransplants$y
ysum <- sum(y)
esum <- sum(e)

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

plot(e,y)

print(ysum/esum)

x <- seq(0, 0.002, length.out=1001)
plot(x, dgamma(x, ysum, esum), type="l")

plot(sort(ppois(y, e*ysum/esum)))