```
# states
= 30
n_states = 1:n_states
i
# target distribution - proportional to
= i ^ 2 * exp(-0.5 * i) # notice: not probabilities
pi_i
= 10000
n_steps = rep(NA, n_steps)
i_sim 1] = 3 # initialize
i_sim[
for (n in 2:n_steps){
= i_sim[n - 1]
current
# propose next state
= sample(c(current + 1, current - 1), size = 1, prob = c(0.5, 0.5))
proposed
# compute acceptance probability
if (!(proposed %in% i)){ # to correct for proposing moves outside of boundaries
= current
proposed
}= min(1, pi_i[proposed] / pi_i[current])
a
# simulate the next state
= sample(c(proposed, current), size = 1, prob = c(a, 1 - a))
i_sim[n] }
```

# Metropolis-Hastings Algorithm

## Island hopping

First few steps

```
# display the first few steps
data.frame(step = 1:n_steps, i_sim) |> head(10) |> kbl() |> kable_styling()
```

step | i_sim |
---|---|

1 | 3 |

2 | 2 |

3 | 2 |

4 | 3 |

5 | 2 |

6 | 2 |

7 | 3 |

8 | 3 |

9 | 4 |

10 | 5 |

First 100 steps

```
# trace plot
plot(1:100, i_sim[1:100], type = "o",
ylim = range(i), xlab = "Day", ylab = "Island")
```

Relative frequencies of first 100 steps

```
# simulated relative frequencies
plot(table(i_sim[1:100]) / 100,
xlab = "Island", ylab = "Proportion of days")
```

First 10000 steps

```
# trace plot
plot(1:n_steps, i_sim, type = "l",
ylim = range(i), xlab = "Day", ylab = "Island")
```

Relative frequencies of first 10000 steps

```
# simulated relative frequencies
plot(table(i_sim) / n_steps,
xlab = "Island", ylab = "Proportion of days")
# target distribution
points(i, pi_i / sum(pi_i), type = "o",
col = "seagreen")
```

## Explicit M-H transition matrix

```
= rbind(
P c(0.5, 0.5, 0),
c(0.25, 0.25, 0.5),
c(0, 1 / 3, 2 / 3)
)
compute_stationary_distribution(P)
```

```
[,1] [,2] [,3]
[1,] 0.1666667 0.3333333 0.5
```

## Zipf

```
# parameters of Zipf
= 1
alpha = 10
M = 1:M
states
# proposal chain is random walk (with reflection at boundaries)
= 0.5
pup = matrix(rep(0, M * M), nrow = M)
Q
1, 1] = 1 - pup
Q[1, 2] = pup
Q[= pup
Q[M, M] - 1] = 1 - pup
Q[M, M for (i in 2:(M - 1)){
- 1] = 1 - pup
Q[i, i + 1] = pup
Q[i, i
}
# Metropolis-Hastings chain
= 10000
N = rep(NA, N)
Xn 1] = 1
Xn[
for (n in 2:N){
# current state
= Xn[n - 1]
i
# proposed state
= sample(states, 1 , prob = Q[i, ])
j
# acceptance probability
= min(1, (1 / j) ^ alpha / (1 / i) ^ alpha)
aij
# next state - accept j or stay at i
= sample(c(i, j), 1, prob = c(1 - aij, aij))
Xn[n] }
```

```
# trace plot
plot(1:N, Xn, type = "l",
ylim = range(states), xlab = "n", ylab = "X_n")
```

```
# simulated relative frequencies
plot(table(Xn) / N,
xlab = "x", ylab = "Relative frequency")
# target Zipf distribution
= (1 / states) ^ alpha
pi_x = pi_x / sum(pi_x)
pi_x
points(states, pi_x, type = "o",
col = "seagreen")
```

## Gamma

Normal proposals

```
= function(u) {
pi_x if (u > 0) {
^ (7.2 - 1) * exp(-3 * u)
u else {
} 0
}
}
= 0.5
sigma
= 10000
N
= rep(NA, N)
x_n 1] = 4
x_n[
for (i in 2:N){
= x_n[i-1]
current = rnorm(1, current, sigma)
proposed = min(1, pi_x(proposed) / pi_x(current))
a = sample(c(current, proposed), 1, prob = c(1 - a, a))
x_n[i] }
```

`plot(1:N, x_n , type="l", xlab = "n", ylab = "X_n")`

```
hist(x_n, freq = FALSE)
curve(dgamma(x, shape = 7.2, rate = 3), add = TRUE, col = "seagreen")
```

Moving Uniform proposals

```
= 0.5
delta
= 10000
N
= rep(NA, N)
x_n 1] = 4
x_n[
for (i in 2:N){
= x_n[i-1]
current = runif(1, current - delta, current + delta)
proposed = min(1, pi_x(proposed) / pi_x(current))
a = sample(c(current, proposed), 1, prob = c(1 - a, a))
x_n[i] }
```

`plot(1:N, x_n , type="l", xlab = "n", ylab = "X_n")`

```
hist(x_n, freq = FALSE)
curve(dgamma(x, shape = 7.2, rate = 3), add = TRUE, col = "seagreen")
```

## Multivariable

Normal proposals

```
= function(x, y) {
pi_xy if ((x > 0) & (y > 0)) {
exp(-(x + y + x * y))
else {
} 0
} }
```

```
= 0.5
sigma
= 10000
N
= rep(NA, N)
x_n = rep(NA, N)
y_n 1] = 1
x_n[1] = 1
y_n[
for (i in 2:N){
= x_n[i-1]
current_x = y_n[i-1]
current_y
= rnorm(1, current_x, sigma)
proposed_x = rnorm(1, current_y, sigma)
proposed_y
= min(1, pi_xy(proposed_x, proposed_y) / pi_xy(current_x, current_y))
a = sample(c("reject", "accept"), 1, prob = c(1 - a, a))
action
if (action == "accept"){
= proposed_x
x_n[i] = proposed_y
y_n[i] else {
} = current_x
x_n[i] = current_y
y_n[i]
}
}
plot(x_n, y_n)
```

```
# Trace plot of first 100 steps
= data.frame(x_n, y_n)
xy
ggplot(xy[1:100, ] |>
mutate(label = 1:100),
aes(x_n, y_n)) +
geom_path() +
geom_point(size = 2) +
geom_text(aes(label = label, x = x_n + 0.1, y = y_n + 0.01)) +
labs(title = "Trace plot of first 100 steps")
```

```
ggplot(xy, aes(x_n, y_n)) +
geom_point(alpha = 0.4) +
stat_density_2d(color = "grey", size = 1)
```

```
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
```