Code
function log_gaussian(x, μ, σ)
-0.5 * ((x - μ) / σ)^2 + log(1 / (sqrt(2 * π) * σ))
end
log_gaussian (generic function with 1 method)
Implementation of MCMC- with the Metropolis algorithm.
We will follow the steps described in the excellent article in introduction to the MCMC of David W. Hogg and Daniel Foreman-Mackey, available here. Section 3 will be particularly useful for this example.
For this first exercise, we will implement the Metropolis algorithm and apply it to a unidimensional normal distribution.
Use the following information:
The equation of normal distribution is
However, to avoid digital errors, you use your logarithm. Code directly a function for
log_gaussian (generic function with 1 method)
We can now implement the Metropolis algorithm.
We want our algorithm to be applicable to any density of (log-) probability which accepts an argument $ Theta $ scalar. We can therefore give log_Density
(our above probability function) in the function argument.
"""
Metropolis MCMC algorithm for sampling from a probability distribution.
Arguments:
- log_density: log-density function, accepts a theta argument
- θ₀: initial value of theta for the MCMC
- n: number of steps to take in the MCMC
- q_scale: standard deviation of the proposal distribution.
"""
function mcmc_metropolis(log_density, θ₀, n, q_scale=1.0, verbose=false)
chain = zeros(n + 1)
chain[1] = θ₀
p_c = log_density(θ₀)
n_accepted = 0
for i in 1:n
proposal = chain[i] + q_scale * randn()
p_i = log_density(proposal)
log_alpha = p_i - p_c
if log(rand()) < log_alpha # Accept the proposal
chain[i+1] = proposal
p_c = p_i
n_accepted += 1
else # Reject the proposal, stay at current position
chain[i+1] = chain[i]
end
end
# Print acceptance rate for diagnostics
verbose && println("Acceptance rate: $(n_accepted / n)")
return chain
end
Main.Notebook.mcmc_metropolis
Apply the algorithm to obtain 10,000 samples.
Then display a histogram and compare it with the analytical PDF.
Then display the temporal evolution (
using CairoMakie
μ = 2.0
σ = sqrt(2.0)
log_density(θ) = log_gaussian(θ, μ, σ)
# Initial conditions
θ₀ = 0.0
nsteps = 10000
q_scale = 1.0
# Run the MCMC algorithm
samples = mcmc_metropolis(log_density, θ₀, nsteps, q_scale)
# Create a figure for the histogram and analytical PDF comparison
f = Figure()
ax1 = Axis(f[1, 1], xlabel="θ", ylabel="Density", title="MCMC Sampling vs Analytical PDF")
hist!(ax1, samples, bins=50, normalization=:pdf, color=(:blue, 0.5), label="MCMC Samples")
# Overlay the analytical PDF
θ_range = range(extrema(samples)..., length=1000)
analytical_pdf = [exp(log_density(θ)) for θ in θ_range]
lines!(ax1, θ_range, analytical_pdf, color=:red, linewidth=2, label="Analytical PDF")
axislegend(ax1)
f
---
title: Lecture lab
subtitle: Implementation of MCMC- with the Metropolis algorithm.
engine: julia
---
We will follow the steps described in the excellent article in introduction to the MCMC of David W. Hogg and Daniel Foreman-Mackey, available [here](https://ui.adsabs.harvard.edu/abs/2018ApJS..236...11H/abstract).
Section 3 will be particularly useful for this example.
## One-dimensional density function (Problems 2 and 3 of the article)
### Gaussian density
For this first exercise, we will implement the Metropolis algorithm and apply it to a unidimensional normal distribution.
Use the following information:
- The density function $P(\Theta)$ is a Gaussian with a dimension with average of $μ = 2 $ and a variance $σ^2 = 2 $.
- The distribution of proposal $q (θ' | \Theta) $ is a Gaussian for $θ'$ with an average $μ = \Theta $ and a standard deviation $σ = 1$.
- The initial point of the MCMC is $ \Theta = 0 $.
- The MCMC must perform $10^4 $ iterations.
The equation of normal distribution is
$$
p(θ) = \frac{1}{\sqrt{2 \pi σ^2}} \exp\left[ -\frac{(θ - μ)^2}{2 σ^2}\right].
$$
However, to avoid digital errors, you use your logarithm. **Code directly a function for $\ln p(θ)$**.
```{julia}
function log_gaussian(x, μ, σ)
-0.5 * ((x - μ) / σ)^2 + log(1 / (sqrt(2 * π) * σ))
end
```
We can now implement the Metropolis algorithm.
We want our algorithm to be applicable to any density of (log-) probability which accepts an argument $ \ Theta $ scalar. We can therefore give `log_Density` (our above probability function) in the function argument.
```{julia}
"""
Metropolis MCMC algorithm for sampling from a probability distribution.
Arguments:
- log_density: log-density function, accepts a theta argument
- θ₀: initial value of theta for the MCMC
- n: number of steps to take in the MCMC
- q_scale: standard deviation of the proposal distribution.
"""
function mcmc_metropolis(log_density, θ₀, n, q_scale=1.0, verbose=false)
chain = zeros(n + 1)
chain[1] = θ₀
p_c = log_density(θ₀)
n_accepted = 0
for i in 1:n
proposal = chain[i] + q_scale * randn()
p_i = log_density(proposal)
log_alpha = p_i - p_c
if log(rand()) < log_alpha # Accept the proposal
chain[i+1] = proposal
p_c = p_i
n_accepted += 1
else # Reject the proposal, stay at current position
chain[i+1] = chain[i]
end
end
# Print acceptance rate for diagnostics
verbose && println("Acceptance rate: $(n_accepted / n)")
return chain
end
```
Apply the algorithm to obtain 10,000 samples.
Then display a histogram and compare it with the analytical PDF.
Then display the temporal evolution ($Θ$ vs $k$) of the MCMC.
```{julia}
using CairoMakie
μ = 2.0
σ = sqrt(2.0)
log_density(θ) = log_gaussian(θ, μ, σ)
# Initial conditions
θ₀ = 0.0
nsteps = 10000
q_scale = 1.0
# Run the MCMC algorithm
samples = mcmc_metropolis(log_density, θ₀, nsteps, q_scale)
# Create a figure for the histogram and analytical PDF comparison
f = Figure()
ax1 = Axis(f[1, 1], xlabel="θ", ylabel="Density", title="MCMC Sampling vs Analytical PDF")
hist!(ax1, samples, bins=50, normalization=:pdf, color=(:blue, 0.5), label="MCMC Samples")
# Overlay the analytical PDF
θ_range = range(extrema(samples)..., length=1000)
analytical_pdf = [exp(log_density(θ)) for θ in θ_range]
lines!(ax1, θ_range, analytical_pdf, color=:red, linewidth=2, label="Analytical PDF")
axislegend(ax1)
f
```
```{julia}
ax2 = Axis(f[2, 1], xlabel="Iteration (k)", ylabel="θ", title="MCMC Trace Plot")
lines!(ax2, 0:nsteps, samples)
f
```
<!-- ```{julia}
burnin = 1000 # Discard the first 1000 samples as burn-in
samples_post_burnin = samples[burnin+1:end]
function autocorrelation(x, lags)
n = length(x)
x_centered = x .- mean(x)
var_x = sum(x_centered .^ 2) / n
ac = zeros(length(lags))
for (i, lag) in enumerate(lags)
ac[i] = sum(x_centered[1:n-lag] .* x_centered[lag+1:n]) / ((n - lag) * var_x)
end
return ac
end
lags = 0:100
ac = autocorrelation(samples_post_burnin, lags)
ax3 = Axis(f[3, 1],
xlabel="Lag",
ylabel="Autocorrelation",
title="Autocorrelation Plot")
lines!(ax3, lags, ac)
f
``` -->