### Metroplis-Hastings Algorithm

## Load data
data <- read.table("data.txt", as.is = TRUE, header = TRUE);
data <- data$y;

## Plot of the data 
hist(data, xlab = "y", ylab = "Density", main = "Mixture data", freq = FALSE, breaks = 20, col = "gray");
lines(density(data, bw = 0.5), lwd = 2, col = "red");

## Density function for the MH ratio
dmistnorm <- function (x, d, m1, m2, s1, s2) {
  d*dnorm(x, m1, s1) + (1 - d)*dnorm(x, m2, s2);
}

## Independent Metropolis
## Objective: build a chain with stationary distribution the posterior of d
## Set the parameters
m1 <- 7;
s1 <- 0.5;
m2 <- 10;
s2 <- 0.5;

## We implement the IMH with 2 different independent proposals
par(mfrow = c(2, 2));
## Case 1: Uniform[0,1]=Beta(1,1) proposal
set.seed(1027675)
BB <- 10000; ## number of iterations of the chain
dt <- rep(NA, BB); ## storage for the values of the chain
dt[1] <- 0.5; ## initial value
accept <- 0 ; ## acceptance count
for (b in 2:BB) {
  ## cat("iteration ", b, "\n")
  dt[b] <- runif(1); ## set dt[b]=proposed value

  ## MH ratio
  logMH <- sum(log(dmistnorm(data, dt[b], m1, m2, s1, s2))) - sum(log(dmistnorm(data, dt[b - 1], m1, m2, s1, s2)));
  MH <- exp(logMH);

  ## Acceptance
  u <- runif(1); ## uniform variable to determine acceptance
  if (u < MH) {
    accept <- accept + 1; ## update count, dt[b] is already fixed
  } else {
    dt[b] <- dt[b - 1];
  }
}

## Plot of the chain
plot(dt, type = "l", xlab = paste("t, acceptance ratio=", accept/BB), ylab = "delta", main = "Uniform proposal");
hist(dt, breaks = 50, main = "Histogram of the sampled values", xlab = "delta");

## Case 2: beta(2,10)
set.seed(1027675);
BB <- 10000; ## number of iterations of the chain
dt <- rep(NA, BB); ## storage for the values of the chain
dt[1] <- 0.5; ## initial value
accept <- 0;  ## acceptance count
for (b in 2:BB) {
  ## cat("iteration ", b, "\n")
  dt[b] <- rbeta(1, 2, 10); ## set dt[b]=proposed value

  ## MH ratio
  logMH <- sum(log(dmistnorm(data, dt[b], m1, m2, s1, s2))) - sum(log(dmistnorm(data, dt[b - 1], m1, m2, s1, s2))) + dbeta(dt[b - 1], 2, 10, log = TRUE) - dbeta(dt[b], 2, 10, log = TRUE);
  MH <- exp(logMH);

  ## Acceptance
  u <- runif(1); ## uniform variable to determine acceptance
  if (u < MH) {
    accept <- accept + 1; ## update count, dt[b] is already fixed
  } else {
    dt[b] <- dt[b - 1];
  }
}

## Plot of the chain
plot(dt, type = "l", xlab = paste("t, acceptance ratio=", accept/BB), ylab = "delta", main = "Beta proposal");
hist(dt, breaks = 50, main = "Histogram of the sampled values", xlab = "delta");

par(mfrow = c(1, 1));

plot(density(rbeta(100000, 2, 10)))