Examples

One-dimensional Metropolis-walk

Note

Even in the one-dimensional case, the initial value must be defined as a vector.

Minimal Working Example for Single-Walker

This is a minimal working example for 1-walker x 1000-steps, 1-dimensional Metropolis-walk on the Gaussian distribution function:

\[p(x) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{1}{2}x^2\right).\]

# sampling Gaussian distribution
using MetropolisAlgorithm
p(x) = exp(-x[1]^2/2) / sqrt(2*π)
r₀ = [1.0]
R = metropolis(p, r₀, n_steps=1000)

1000-element Vector{Vector{Float64}}:
 [1.0]
 [1.0157365327985188]
 [1.0157365327985188]
 [0.6514845979353355]
 [0.6146693755218664]
 [0.9867345779027923]
 [0.6403736493893666]
 [0.5903570812963379]
 [0.26502028287031665]
 [0.4078178283408803]
 ⋮
 [-0.6655601066911362]
 [-0.35265099558962765]
 [-0.4808902444153368]
 [-0.10984476259431197]
 [0.06203716536661985]
 [0.017678447545500298]
 [-0.051390131923333304]
 [0.017075758215570613]
 [0.24496040309502198]

Convert from Vector{Vector{Float64}} to Vector{Float64} for plotting R. A histogram of the output trajectory data R should be consistent with the input distribution function p. Consistency is confirmed in another example.

# plotting histogram
using CairoMakie
X = [r[1] for r in R]
hist(X)

This is the trajectory of a walker at each step. The histogram above shows the number of these points in each bin.

# plotting trajectory
using CairoMakie
X = [r[1] for r in R]
Y = keys(X)
lines(X, Y)

Using Distributions.jl

Here is an example of sampling the distribution functions in Distributions.jl.

# distribution
using Distributions
d = Normal(0, 1)

# sampling
using MetropolisAlgorithm
R = metropolis(x -> Distributions.pdf(d,x[1]), [1.0], n_steps=10000, d=1.0)

# reshape for plotting
X = [r[1] for r in R]
Y = keys(X)  .- 1

# figure
using CairoMakie
fig = Figure(
  size = (420,600),
  fontsize = 11,
  backgroundcolor = :transparent,
)

# histogram
axis = Axis(
  fig[1,1],
  limits = (-5, 5, 0, 1.1*Distributions.pdf(d,d.μ)),
  titlesize = 16.5,
  xlabelsize = 16.5,
  ylabelsize = 16.5,
  title = "Histogram",
  xlabel = "x",
  ylabel = "PDF(x)",
  backgroundcolor = :transparent,
)
hist!(axis, [first(r) for r in R], label = "Metropolis", bins = 50, normalization = :pdf)
lines!(axis, -50..50, x -> Distributions.pdf(d,x), label = "Exact", color=:black)
axislegend(axis, position = :rt, framevisible = false)

# trajectory
axis = Axis(
  fig[2,1],
  limits = (-5, 5, 0, length(R)),
  titlesize = 16.5,
  xlabelsize = 16.5,
  ylabelsize = 16.5,
  title = "Trajectory",
  xlabel = "x",
  ylabel = "steps",
  backgroundcolor = :transparent,
)
lines!(axis, X, Y, linewidth = 0.3, label = "Metropolis")
axislegend(axis, position = :rt, framevisible = false)

# display
fig

The output histograms are consistent with the input distribution functions.

# packages
using CairoMakie
using Distributions
using MetropolisAlgorithm

# initialize
fig = Figure(
  size = (1260,600),
  fontsize = 11,
  backgroundcolor = :transparent,
)

for n in 1:6

    # distribution
    d = [
        Normal(0, 1)
        SymTriangularDist(0, 1)
        Uniform(0, 1)
        Gamma(7.5, 1)
        TriangularDist(0, 1, 0.2)
        Semicircle(1)
    ][n]
    μ = Distributions.mean(d)
    σ = Distributions.std(d)

    # sampling
    R = metropolis(x -> Distributions.pdf(d,x[1]), [1.0], n_steps=100000, d=σ)

    # plot
    axis = Axis(
      fig[div(n-1,3)+1,rem(n-1,3)+1],
      limits = (-5*σ+μ, 5*σ+μ, 0, 1.2*maximum(Distributions.pdf(d,x) for x in -5*σ+μ:0.1:5*σ+μ)),
      titlesize = 16.5,
      xlabelsize = 16.5,
      ylabelsize = 16.5,
      title = "$d",
      xlabel = "x",
      ylabel = "PDF(x)",
      backgroundcolor = :transparent,
    )
    hist!(axis, [first(r) for r in R], label = "Metropolis", bins = 50, normalization = :pdf)
    lines!(axis, -50..50, x -> Distributions.pdf(d,x), label = "PDF", color=:black)
    axislegend(axis, position = :rt, framevisible = false)

end

fig

Minimal Working Example for Multiple-Walkers

Allocate memory by yourself for multiple walkers.

using MetropolisAlgorithm
p(x) = exp(-x[1]^2/2) / sqrt(2*π)
R = fill([0.0], 10000)

10000-element Vector{Vector{Float64}}:
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 ⋮
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]
 [0.0]

Each step is run without memory allocation for walkers, overwriting the second argument.

metropolis!(p, R)
R

10000-element Vector{Vector{Float64}}:
 [-0.2514748259856592]
 [0.32825036500794813]
 [-0.03583459422794444]
 [0.3532503225132929]
 [-0.27327407899879785]
 [-0.37984673252517]
 [-0.24799443078273364]
 [-0.06817789577916122]
 [0.04901243286515822]
 [-0.4408397518595978]
 ⋮
 [-0.29837274824980076]
 [-0.1404610758358339]
 [0.3211355141120935]
 [0.48062292998437783]
 [-0.019466274768758174]
 [0.09471402774707927]
 [0.45849977929853414]
 [-0.32742668523375174]
 [0.1550412809173819]

Use the For statement to repeat as many times as you like.

for i in 1:100
  metropolis!(p, R)
end
R

10000-element Vector{Vector{Float64}}:
 [1.5286739598275951]
 [-0.600630162229513]
 [-1.3559511305331502]
 [0.2339331246680385]
 [-0.7082985403089351]
 [2.5443828189580713]
 [0.639206146556351]
 [1.4335768674575626]
 [0.24294439012767388]
 [-0.23219187329494484]
 ⋮
 [0.5082769464064576]
 [-0.20212887554452874]
 [0.6794534897162369]
 [2.6511294432070787]
 [0.9404820025817854]
 [-0.34688263823256216]
 [0.24394584490990912]
 [-1.014657930525833]
 [1.0818928172918554]

Time evolution to reach equilibrium. The first several steps are not consistent with the correct distribution.

using MetropolisAlgorithm

# distribution function
p(x) = exp(-x[1]^2/2) / sqrt(2*π)

# figure
using CairoMakie
fig = Figure(size=(1680, 420))

# axis
axis1 = Axis(fig[1,1], limits=(-5, 5, 0, 1.1*p([0])), title="n=0")
axis2 = Axis(fig[1,2], limits=(-5, 5, 0, 1.1*p([0])), title="n=1")
axis3 = Axis(fig[1,3], limits=(-5, 5, 0, 1.1*p([0])), title="n=2")
axis4 = Axis(fig[1,4], limits=(-5, 5, 0, 1.1*p([0])), title="n=100")

# n = 0
R = fill(zeros(1), 10000)
hist!(axis1, [r[1] for r in R], bins=-5:0.1:5, normalization=:pdf)
lines!(axis1, -5..5, p, color=:black)

# n = 1
metropolis!(p, R)
hist!(axis2, [r[1] for r in R], bins=-5:0.1:5, normalization=:pdf)
lines!(axis2, -5..5, p, color=:black)

# n = 2
metropolis!(p, R)
hist!(axis3, [r[1] for r in R], bins=-5:0.1:5, normalization=:pdf)
lines!(axis3, -5..5, p, color=:black)

# n = 100
for i in 3:100
    metropolis!(p, R)
end
hist!(axis4, [r[1] for r in R], bins=-5:0.1:5, normalization=:pdf)
lines!(axis4, -5..5, p, color=:black)

# display
fig

Three-dimensional Metropolis-walk

Here is an example of sampling a function like an atomic orbital (d-orbital).

Single-Walker

# sampling
using MetropolisAlgorithm
ψ(r) = r[1] * r[2] * exp(- r[1]^2 - r[2]^2 - r[3]^2)
p(r) = abs2(ψ(r))
R = metropolis(r -> abs2(ψ(r)), [1.0, 0.0, 0.0], n_steps=50000)

# plot
using CairoMakie
CairoMakie.activate!(type = "png")
fig = Figure(size=(420,420), figure_padding=0)
axis = Axis(fig[1,1], aspect=1, backgroundcolor=:black, limits=(-2,2,-2,2))
hidespines!(axis)
hidedecorations!(axis)
lines!(axis, [r[1] for r in R], [r[2] for r in R], linewidth=0.1, color="#00FFFF")
fig

Multiple-Walkers

using MetropolisAlgorithm

# distribution function
ψ(r) = r[1] * r[2] * exp(- r[1]^2 - r[2]^2 - r[3]^2)
p(r) = abs2(ψ(r))

# figure
using CairoMakie
CairoMakie.activate!(type = "png")
fig = Figure(size=(1680, 420))

# axis
axis1 = Axis(fig[1,1], aspect=1, limits=(-2,2,-2,2), backgroundcolor=:black, title="n=0")
axis2 = Axis(fig[1,2], aspect=1, limits=(-2,2,-2,2), backgroundcolor=:black, title="n=1")
axis3 = Axis(fig[1,3], aspect=1, limits=(-2,2,-2,2), backgroundcolor=:black, title="n=2")
axis4 = Axis(fig[1,4], aspect=1, limits=(-2,2,-2,2), backgroundcolor=:black, title="n=10")
hidespines!(axis1)
hidespines!(axis2)
hidespines!(axis3)
hidespines!(axis4)
hidedecorations!(axis1)
hidedecorations!(axis2)
hidedecorations!(axis3)
hidedecorations!(axis4)

# n = 0
R = fill(zeros(3), 10000)
scatter!(axis1, [r[1] for r in R], [r[2] for r in R], markersize=2, color="#00FFFF")

# n = 1
metropolis!(p, R)
scatter!(axis2, [r[1] for r in R], [r[2] for r in R], markersize=2, color="#00FFFF")

# n = 2
metropolis!(p, R)
scatter!(axis3, [r[1] for r in R], [r[2] for r in R], markersize=2, color="#00FFFF")

# n = 10
for i in 1:10
    metropolis!(p, R)
end
scatter!(axis4, [r[1] for r in R], [r[2] for r in R], markersize=2, color="#00FFFF")

# display
fig