Sampling Trajectories

Examples of the sample_trajectory command with different samplers.

RWM Example

using Plots
using Printf
using Random
using BasicMD

function V(x)
    return (x[1]^2 -1)^2
end

β = 5.0;
x0 = [-1.0];
seed = 100;
Δt = 1e-1;
n_iters = 10^4; # number of samples

sampler = RWM(V, β, Δt);

Random.seed!(100);

xvals, avals = sample_trajectory(x0, sampler, options=MDOptions(n_iters=n_iters));

histogram([x_[1] for x_ in xvals],label="RWM Samples",normalize=true)
xlabel!("x")
ylabel!("Frequency")
Example block output

HMC Example

using Plots
using Printf
using Random
using BasicMD
using ForwardDiff

function V(x)
    return (x[1]^2 -1)^2
end

gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x);


β = 5.0;
x0 = [-1.0];
seed = 100;
M = 1.;
nΔt = 10^1; # number of Verlet steps per HMC iteration
Δt = 1e-1;
n_iters = 10^4; # number of samples

sampler = HMC(V, gradV!, β, M, Δt, nΔt);

Random.seed!(100);

xvals, avals = sample_trajectory(x0, sampler, options=MDOptions(n_iters=n_iters));

histogram([x_[1] for x_ in xvals],label="HMC Samples",normalize=true)
xlabel!("x")
ylabel!("Frequency")
Example block output

Euler-Maruyama Example

using Plots
using Printf
using Random
using BasicMD
using ForwardDiff

function V(x)
    return (x[1]^2 -1)^2
end

gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x);


β = 5.0;
x0 = [-1.0];
seed = 100;
Δt = 1e-1;
n_iters = 10^4; # number of samples

sampler = EM(gradV!, β, Δt);

Random.seed!(100);

xvals = sample_trajectory(x0, sampler, options=MDOptions(n_iters=n_iters));

histogram([x_[1] for x_ in xvals],label="EM Samples",normalize=true)
xlabel!("x")
ylabel!("Frequency")
Example block output

ABOBA Example

Coordinates for Inertial Samplers

For inertial samplers, which approximate the underdamped Langevin equation, since the trajectory $x(t) = (q(t), p(t))$, it is neccessary to provide both an initial position and an initial momentum.

using Plots
using Printf
using Random
using BasicMD
using ForwardDiff

function V(x)
    return (x[1]^2 -1)^2
end

gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x);


β = 5.0;
γ = 1.;
M = 1.;
q0 = [-1.0];
p0 = [0.0]; # as ABOBA is an inertial sampler, we need to specify a momentum
x0 = [copy(q0), copy(p0)];
seed = 100;
Δt = 1e-1;
n_iters = 10^4; # number of samples

sampler = ABOBA(gradV!, β, γ, M, Δt);
opts = MDOptions(n_iters=n_iters);

Random.seed!(100);
xvals = sample_trajectory(x0, sampler, options= opts);

histogram([x_[1][1] for x_ in xvals],label="ABOBA Samples",normalize=true)
xlabel!("x")
ylabel!("Frequency")
Example block output