# Load packages
library(ggplot2)
library(patchwork)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Exercise 7.11 (Changing the proposal model)

one_mh_iteration_normal <- function(s, current){
  # STEP 1: Propose the next chain location
  proposal <- rnorm(1,current,s)
  
  # STEP 2: Decide whether or not to go there
  proposal_plaus <- dnorm(proposal, 0, 1) * dnorm(6.25, proposal, 0.75)
  current_plaus  <- dnorm(current, 0, 1) * dnorm(6.25, current, 0.75)
  alpha <- min(1, proposal_plaus / current_plaus)
  next_stop <- sample(c(proposal, current), 
                      size = 1, prob = c(alpha, 1-alpha))
  
  # Return the results
  return(data.frame(proposal, alpha, next_stop))
}

#corremos la función con los valores que se propone en cada parte del ejercicio

#a)
set.seed(1)
one_mh_iteration_normal(s = 0.01, current = 3)
##   proposal     alpha next_stop
## 1 2.993735 0.9826955  2.993735
#b)
set.seed(1)
one_mh_iteration_normal(s = 0.5, current = 3)
##   proposal     alpha next_stop
## 1 2.686773 0.3655544         3
#c)
set.seed(1)
one_mh_iteration_normal(s = 1, current = 3)
##   proposal     alpha next_stop
## 1 2.373546 0.1017526         3
#d)
set.seed(1)
one_mh_iteration_normal(s = 3, current = 3)
##   proposal        alpha next_stop
## 1 1.120639 4.002513e-05         3

Exercise 7.12 (Metropolis-Hastings tour with Normal proposals)

mh_tour_normal <- function(N, s){
  # 1. Start the chain at location 3
  current <- 3
  
  # 2. Initialize the simulation
  mu <- rep(0, N)
  
  # 3. Simulate N Markov chain stops
  for(i in 1:N){    
    # Simulate one iteration
    sim <- one_mh_iteration_normal(s = s, current = current)
    
    # Record next location
    mu[i] <- sim$next_stop
    
    # Reset the current location
    current <- sim$next_stop
  }
  
  # 4. Return the chain locations
  return(data.frame(iteration = c(1:N), mu))
}
  
#a)20 iterations, s=0.01
set.seed(84735)
mh_simulation_1 <- mh_tour_normal(N = 20, s = 0.01)

t<-ggplot(mh_simulation_1, aes(x = iteration, y = mu)) + 
  geom_line() 

h<-ggplot(mh_simulation_1, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 20) + 
  stat_function(fun = dnorm, args = list(4,0.6), color = "blue")
t+h
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#b)20 iterations, s=10
set.seed(84735)
mh_simulation_1 <- mh_tour_normal(N = 20, s = 10)

t<-ggplot(mh_simulation_1, aes(x = iteration, y = mu)) + 
  geom_line()

h<-ggplot(mh_simulation_1, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 20) + 
  stat_function(fun = dnorm, args = list(4,0.6), color = "blue")
t+h

#c)1000 iterations, s=0.01
set.seed(84735)
mh_simulation_1 <- mh_tour_normal(N = 1000, s = 0.01)

t<-ggplot(mh_simulation_1, aes(x = iteration, y = mu)) + 
  geom_line()

h<-ggplot(mh_simulation_1, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 20) + 
  stat_function(fun = dnorm, args = list(4,0.6), color = "blue")
t+h

#d)1000 iterations, s=10
set.seed(84735)
mh_simulation_1 <- mh_tour_normal(N = 1000, s = 10)

t<-ggplot(mh_simulation_1, aes(x = iteration, y = mu)) + 
  geom_line()

h<-ggplot(mh_simulation_1, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 20) + 
  stat_function(fun = dnorm, args = list(4,0.6), color = "blue")
t+h

#f)
set.seed(84735)
mh_simulation_1 <- mh_tour_normal(N = 1000, s = 1)

t<-ggplot(mh_simulation_1, aes(x = iteration, y = mu)) + 
  geom_line()

h<-ggplot(mh_simulation_1, aes(x = mu)) + 
  geom_histogram(aes(y = ..density..), color = "white", bins = 20) + 
  stat_function(fun = dnorm, args = list(4,0.6), color = "blue")
t+h

## Exercise 7.13 (Change the Normal prior)

new_mh_iteration <- function(w, current,m,s){
  # STEP 1: Propose the next chain location
  proposal <- runif(1, min = current - w, max = current + w)
  
  # STEP 2: Decide whether or not to go there
  proposal_plaus <- dnorm(proposal, m, s) * dnorm(6.25, proposal, 0.75)
  current_plaus  <- dnorm(current, m, s) * dnorm(6.25, current, 0.75)
  alpha <- min(1, proposal_plaus / current_plaus)
  next_stop <- sample(c(proposal, current), 
                      size = 1, prob = c(alpha, 1-alpha))
  
  # Return the results
  return(data.frame(proposal, alpha, next_stop))
}
#a)
set.seed(84735)
new_mh_iteration(w = 1, current = 3, m = 0, s = 10)
##   proposal alpha next_stop
## 1 3.495377     1  3.495377
#b)
set.seed(84735)
new_mh_iteration(w = 1, current = 3, m = 20, s = 1)
##   proposal alpha next_stop
## 1 3.495377     1  3.495377
#c)
set.seed(84735)
new_mh_iteration(w = 0.1, current = 3, m = 20, s = 1)
##   proposal alpha next_stop
## 1 3.049538     1  3.049538
#d)
set.seed(84735)
new_mh_iteration(w = 0.1, current = 3, m = -15, s = 10)
##   proposal alpha next_stop
## 1 3.049538     1  3.049538