# Independence Metropolis-Hasting Examples from class on
# February 27 and 28, 2019.

# Target density is the beta(alpha, beta) density
# Candidate density is the Uniform(0,1)


rbeta.mh <- function(N, X0, alpha, beta) {
    Y <- runif(N) # candidates
    U <- runif(N) # Used for acceptance testing
    X <- c(X0, rep(NA,N))
    for (t in 1:N) {
       Xt <- X[t]
       Yt <- Y[t]
       pt <- (dbeta(Yt, alpha, beta)/dunif(Yt))*
             (dunif(Xt)/dbeta(Xt, alpha, beta))
       Ut <- U[t]
       X[t+1] <- ifelse(Ut<pt, Yt, Xt)
    }
    X[2:(N+1)]
}

# Compare this with accept-reject sampling with N the number of candidates

rbeta.ar1 <- function(N, alpha, beta) {
    if ((alpha<1)||(beta<1))
       stop("Both parameters must be at least 1")
    x0 <- (alpha-1)/(alpha+beta-2) 
    M <- dbeta(x0, alpha, beta)
    Y <- runif(N) # candidates
    U <- runif(N) # Used for acceptance testing
    p.accept <- dbeta(Y, alpha, beta)/(M*dunif(Y))
    X <- Y[U<p.accept]
    X
}

# Example 2: Generate from a Standard Cauchy using an 
# independence M-H algorithm with a standard normal proposal.
rcauchy.mh <- function(N, X0) {
    Y <- rnorm(N)
    U <- runif(N)
    X <- c(X0,rep(NA,N))
    for (t in 1:N) {
       Xt <- X[t]
       Yt <- Y[t]
       pt <- (dcauchy(Yt)/dnorm(Yt))*(dnorm(Xt)/dcauchy(Xt))
       Ut <- U[t]
       X[t+1] <- ifelse(Ut<pt, Yt, Xt)
    }
    X[2:(N+1)]
}

# Example 3: Generate from a Standard Cauchy using an 
# Random Walk M-H algorithm with a standard normal error distribution.
rcauchy.rwmh <- function(N, X0, sigma=1) {
    eps <- rnorm(N,0,sigma)
    U <- runif(N)
    X <- c(X0,rep(NA,N))
    for (t in 1:N) {
       Xt <- X[t]
       Yt <- Xt+eps[t]
       pt <- dcauchy(Yt)/dcauchy(Xt)
       Ut <- U[t]
       X[t+1] <- ifelse(Ut<pt, Yt, Xt)
    }
    X[2:(N+1)]
}