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.345286839829905]
 [1.5973517280440164]
 [1.5973517280440164]
 [1.7217537757384416]
 [1.3552008338346235]
 [1.2296768333719443]
 [1.5031779371736174]
 [1.5031779371736174]
 [1.2283229383986507]
 ⋮
 [1.0327768683058869]
 [1.0601341918195426]
 [1.2432469549008571]
 [0.8970601186890979]
 [0.6152672219162981]
 [1.0152197715157667]
 [0.6165087236508786]
 [0.5371294397706108]
 [0.8403069460351947]

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.4636944236329489]
 [0.41322398208421174]
 [0.10814446572044956]
 [0.421704728046172]
 [-0.012346235879550305]
 [-0.3361322637537677]
 [-0.20772433744733587]
 [0.015541414256029396]
 [-0.22804772354139358]
 [0.23162982517249053]
 ⋮
 [-0.001926000288464902]
 [-0.20169828821390234]
 [0.3703927237197855]
 [0.3416334385768399]
 [0.38152590241118445]
 [0.06494657977401075]
 [0.0838206494189927]
 [-0.3948670832002573]
 [0.13180613759407]

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}}:
 [-0.6852511537656505]
 [-1.0753418115664584]
 [-1.5667480827702644]
 [-0.049398570485480464]
 [-3.1529552024192347]
 [-1.8260384493665347]
 [0.8480229226523952]
 [-0.6755543232575186]
 [0.824607998141009]
 [-1.2623667414635253]
 ⋮
 [-0.8555137528145933]
 [-1.6730399538813625]
 [2.4404010476503504]
 [-0.4910469963297216]
 [-1.9081146786413221]
 [-0.17632956562059676]
 [0.8269867914037929]
 [-0.33064694180217713]
 [0.39205195237546564]

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