bayesian-sampling-methods/bayesian_sampling_methods.R at main · eeshalriz/bayesian-sampling-methods · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
# =============================================================================
# Bayesian Computational Methods: Simulation & Sampling
# Author:  Eeshal Rizwan
# Course:  Bayesian Inference — Heriot-Watt University
# =============================================================================
# Three sampling methods implemented from scratch in R:
#   1. Inverse CDF (Inverse Transform) Sampling
#   2. Rejection Sampling (Accept-Reject method)
#   3. Metropolis-Hastings MCMC Algorithm
# =============================================================================

set.seed(123)
n <- 10000  # number of samples

# =============================================================================
# Method 1: Inverse CDF Sampling
# =============================================================================
# Target distribution: f(x) = 1 / (exp(1) - 1) * exp(x),  x in [0, 1]
# CDF inverse derived analytically: X = log(U * (e - 1) + 1)
# where U ~ Uniform(0, 1)

U <- runif(n)
X_inv <- log(U * (exp(1) - 1) + 1)

# Visualise
hist(X_inv,
     probability = TRUE,
     breaks      = 40,
     main        = "Method 1: Inverse CDF Sampling",
     xlab        = "x",
     col         = "steelblue",
     border      = "white")

# Overlay true density for verification
curve(exp(x) / (exp(1) - 1), from = 0, to = 1,
      col = "red", lwd = 2, add = TRUE)
legend("topleft", legend = c("Sampled", "True density"),
       col = c("steelblue", "red"), lty = c(NA, 1),
       pch = c(15, NA), bty = "n")

cat("Method 1 — Inverse CDF Sampling\n")
cat(sprintf("  Sample mean:     %.4f\n", mean(X_inv)))
cat(sprintf("  Sample variance: %.4f\n\n", var(X_inv)))

# =============================================================================
# Method 2: Rejection Sampling (Accept-Reject)
# =============================================================================
# Target distribution: f(x) = 4x^3,  x in [0, 1]
# Proposal distribution: g(x) = Uniform(0, 1)
# Envelope constant M = 4  (since max of f(x) / g(x) = 4)
# Accept candidate X* if U < f(X*) / (M * g(X*)) = X*^3 / 1 = X*^3
# Acceptance condition: U < 4 * X*^3

X_rej <- numeric(n)
count  <- 0
total_proposals <- 0  # track efficiency

while (count < n) {
  X_star <- runif(1)          # candidate from g(x) = Uniform(0,1)
  U      <- runif(1)          # uniform draw for accept/reject decision
  total_proposals <- total_proposals + 1

  if (U < 4 * X_star^3) {    # acceptance condition
    count        <- count + 1
    X_rej[count] <- X_star
  }
}

# Visualise
hist(X_rej,
     probability = TRUE,
     breaks      = 40,
     main        = "Method 2: Rejection Sampling",
     xlab        = "x",
     col         = "darkorange",
     border      = "white")

curve(4 * x^3, from = 0, to = 1,
      col = "red", lwd = 2, add = TRUE)
legend("topleft", legend = c("Sampled", "True density f(x) = 4x³"),
       col = c("darkorange", "red"), lty = c(NA, 1),
       pch = c(15, NA), bty = "n")

cat("Method 2 — Rejection Sampling\n")
cat(sprintf("  Acceptance rate: %.1f%%\n", 100 * n / total_proposals))
cat(sprintf("  Sample mean:     %.4f\n", mean(X_rej)))
cat(sprintf("  Sample variance: %.4f\n\n", var(X_rej)))

# =============================================================================
# Method 3: Metropolis-Hastings MCMC
# =============================================================================
# Target distribution: f(x) ∝ exp(-x),  x >= 0  (standard Exponential)
# Proposal distribution: Normal(X_current, sigma^2) — truncated to x >= 0
# Acceptance ratio: alpha = min(1, f(X*) / f(X_current))
#                         = min(1, exp(-X* + X_current))

X_mh    <- numeric(n)
X_mh[1] <- 0       # initialise chain
sigma   <- 0.5     # proposal standard deviation (tuning parameter)

accepted <- 0      # track acceptance rate

for (i in 2:n) {
  X_star <- rnorm(1, X_mh[i - 1], sigma)  # propose new value

  if (X_star >= 0) {
    # Log acceptance ratio (numerically stable)
    log_alpha <- -X_star + X_mh[i - 1]
    alpha     <- exp(log_alpha)

    if (runif(1) < alpha) {
      X_mh[i] <- X_star   # accept proposed value
      accepted <- accepted + 1
    } else {
      X_mh[i] <- X_mh[i - 1]  # reject, stay at current value
    }
  } else {
    X_mh[i] <- X_mh[i - 1]    # reject out-of-support proposals
  }
}

# Visualise
hist(X_mh,
     probability = TRUE,
     breaks      = 60,
     main        = "Method 3: Metropolis-Hastings MCMC",
     xlab        = "x",
     col         = "mediumseagreen",
     border      = "white")

curve(exp(-x), from = 0, to = max(X_mh),
      col = "red", lwd = 2, add = TRUE)
legend("topright", legend = c("Sampled", "True density f(x) = exp(-x)"),
       col = c("mediumseagreen", "red"), lty = c(NA, 1),
       pch = c(15, NA), bty = "n")

cat("Method 3 — Metropolis-Hastings MCMC\n")
cat(sprintf("  Acceptance rate: %.1f%%\n", 100 * accepted / (n - 1)))
cat(sprintf("  Sample mean:     %.4f  (theoretical: 1.0000)\n", mean(X_mh)))
cat(sprintf("  Sample variance: %.4f  (theoretical: 1.0000)\n\n", var(X_mh)))

# =============================================================================
# Summary comparison
# =============================================================================
cat("=== Summary ===\n")
cat(sprintf("  Inverse CDF  | mean = %.4f | var = %.4f\n", mean(X_inv), var(X_inv)))
cat(sprintf("  Rejection    | mean = %.4f | var = %.4f\n", mean(X_rej), var(X_rej)))
cat(sprintf("  MH MCMC      | mean = %.4f | var = %.4f\n", mean(X_mh),  var(X_mh)))