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

  • Misalkan sebaran posterior bagi \(\theta\) adalah: \(g(\theta \mid y) = e^{-\theta}\) dengan \(\theta>0\)

  • 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 = 1000000

set.seed(1111)
N = 1000000
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\), Ragam, dan Galat Bakunya

# Penduga titik bagi theta
round(mean(theta),7)
## [1] 1.005546
# Ragam dari penduga theta
round(var(theta),7)
## [1] 1.020186
# Galat baku (standard error) bagi theta
round(sqrt(var(theta)),7)
## [1] 1.010042

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.02484 3.71061]

Perhatian:

  • Buktikan secara analitik (kalkulus): penduga bagi \(\theta\), ragam, dan galat bakunya.
  • Kemudian bandingkan dengan hasil dari metode MCMC.