STA1312: Metode MCMC dengan Algoritma Metropolis–Hastings

par(mar = c(3, 3, 3, 3))

• Misalkan sebaran posterior bagi \(\theta\) adalah \(e^{-\theta}\)

• Akan dilakukan sampling berdasarkan sebaran posterior tersebut

• Sebaran posterior dinyatakan sebagai sebaran target

target <- function(s){
  if (s < 0) {
    return(0)
  } else {
    return(exp(-s))
  }
}

• Algoritma Metropolis–Hastings untuk mendapatkan sampel yang sebarannya proporsional pada target

• Simulasi sebanyak N = 100000

set.seed(1111)
N = 100000
theta <- rep(0, N)
theta[1] <- 3 # Sebagai nilai inisial bagi theta

for (i in 2:N) {
  current.theta <- theta[i - 1]
  proposal.dtheta <- current.theta + rnorm(1, mean = 0, sd = 1) # Random-Walk
  A <- target(proposal.dtheta) / target(current.theta)
  if (runif(1) < A) {
    theta[i] <- proposal.dtheta
  } else {
    theta[i] <- current.theta
  }
}

• Nilai \(\theta\) dari proses di atas merupakan realisasi dari Markov Chain

plot(theta)

hist(theta)

Inferensi Bayesian: Penduga Theta

round(mean(theta),5)

## [1] 0.98065

round(var(theta),5)

## [1] 0.91798

Inferensi Bayesian: Menentukan 95% Credible Interval

nilai.kuantil <- quantile(theta, probs = c(0.025, 0.975))
cat(paste("Pendekatan 95% credible interval : ["
          , round(nilai.kuantil[1], 5), " ", round(nilai.kuantil[2], 5), "]\n", sep = ""))

## Pendekatan 95% credible interval : [0.02481 3.55896]