Monte Carlo Markov Chains

Monte Carlo Markov Chain (MCMC) algorithms are essential for sampling from probability distributions when direct sampling is challenging. Let’s dive into implementing one of the most popular MCMC algorithms, the Metropolis-Hastings algorithm, in Julia.

Step-by-Step Tutorial on Implementing Metropolis-Hastings in Julia

1. Setup Julia Environment

First, ensure you have Julia installed. You can download it from here.

using Random
using Distributions

2. Define the Target Distribution

Let’s assume we are sampling from a normal distribution with a mean of 0 and a standard deviation of 1.

target_distribution(x) = pdf(Normal(0, 1), x)

3. Metropolis-Hastings Algorithm Implementation

We’ll define the Metropolis-Hastings function. This function will iterate through the MCMC steps to generate a chain of samples.

function metropolis_hastings(target_distribution, initial_value, num_samples, proposal_std)
    samples = Vector{Float64}(undef, num_samples)
    current_sample = initial_value
    
    for i in 1:num_samples
        proposal = current_sample + rand(Normal(0, proposal_std))
        acceptance_ratio = target_distribution(proposal) / target_distribution(current_sample)
        
        if rand() < acceptance_ratio
            current_sample = proposal
        end
        
        samples[i] = current_sample
    end
    
    return samples
end

4. Run the Algorithm

Now, let’s generate samples using our Metropolis-Hastings implementation.

initial_value = 0.0
num_samples = 10000
proposal_std = 1.0

samples = metropolis_hastings(target_distribution, initial_value, num_samples, proposal_std)

5. Visualize the Results

To visualize the results, we can plot the histogram of the generated samples and compare it with the target distribution.

using Plots

histogram(samples, normalize = true, label = "MCMC Samples")
x = -4:0.1:4
plot!(x, pdf.(Normal(0, 1), x), label = "Target Distribution", linewidth = 2)

6. Interpreting the Results

You should observe that the histogram of the samples closely matches the target distribution, confirming that our Metropolis-Hastings algorithm is working correctly.

7. Next Steps

• Parameter Tuning: Adjust the proposal standard deviation (proposal_std) and the number of samples (num_samples) to observe the effect on the results.
• Different Distributions: Try implementing the algorithm for different target distributions.
• Advanced Algorithms: Explore other MCMC algorithms like Gibbs sampling, Hamiltonian Monte Carlo, etc.

Additional Resources

• MCMC in Julia
• JuliaLang Documentation

---

The Metropolis-Hastings Algorithm

The Metropolis-Hastings algorithm is a popular Monte Carlo Markov Chain (MCMC) method used for generating a sequence of samples from a probability distribution when direct sampling is challenging. It aims to approximate the target distribution by creating a Markov chain that has the desired distribution as its equilibrium distribution.

Key Steps in the Metropolis-Hastings Algorithm:

1. Initialization:
   • Start with an initial value \(x_0\) from the target distribution.
2. Proposal Distribution:
   • Define a proposal distribution \(q(x'|x)\), which suggests a new state \(x'\) given the current state \(x\).
3. Acceptance Ratio:
   • Compute the acceptance ratio \(\alpha\): \[ \alpha = \min\left(1, \frac{p(x') \cdot q(x|x')}{p(x) \cdot q(x'|x)}\right) \] where \(p(x)\) is the target distribution.
4. Transition:
   • Generate a uniform random number \(u\) between 0 and 1.
   • If \(u \leq \alpha\), accept the proposed state \(x'\) and set \(x_{i+1} = x'\).
   • Otherwise, reject the proposed state and set \(x_{i+1} = x_i\).
5. Iteration:
   • Repeat the above steps for a large number of iterations to create a sequence of samples.

Intuition:

• The algorithm ensures that the chain spends more time in regions of higher probability density.
• Over time, the distribution of the samples converges to the target distribution, allowing us to approximate it effectively.

Benefits:

• Flexibility: Can be used for various complex distributions.
• Efficiency: Provides a way to sample from high-dimensional distributions where direct sampling is impractical.

Applications:

• Bayesian Inference: To estimate posterior distributions.
• Statistical Physics: For simulating systems at thermal equilibrium.
• Machine Learning: Used in various probabilistic models and algorithms.

---