DifferentialEvolutionMCMC
DifferentialEvolutionMCMC.jl is a Differential Evolution MCMC sampler written in Julia and uses the AbstractMCMC interface. DifferentialEvolutionMCMC.jl works with any model, provided that it returns an exact or approximate log likeilhood. An annotated example is provided below. Other examples can be found in the examples subfolder.
Example
First, load the required libraries.
using DifferentialEvolutionMCMC, Random, Distributions
Random.seed!(50514)
Define a function that returns the prior log likelihood of parameters μ and σ. Note that order matters for parameters throughout. The algorithm expects parameters to have the same order.
function prior_loglike(μ, σ)
LL = 0.0
LL += logpdf(Normal(0, 1), μ)
LL += logpdf(truncated(Cauchy(0, 1), 0, Inf), σ)
return LL
end
Define a function for the initial sample. Sampling from the prior distribution is a reasonable choice for most applications.
function sample_prior()
μ = rand(Normal(0, 1))
σ = rand(truncated(Cauchy(0, 1), 0, Inf))
return [μ,σ]
end
Next, define a function for the log likelihood which accepts the data followed by the parameters (in the order specififed in the priors).
function loglike(data, μ, σ)
return sum(logpdf.(Normal(μ, σ), data))
end
Specify the upper and lower bounds of the parameters.
bounds = ((-Inf,Inf),(0.0,Inf))
Define the names of parameters. Elements of parameter vectors do not need to be named.
names = (:μ,:σ)
Generate simulated data from a normal distribution
data = rand(Normal(0.0, 1.0), 50)
Now we will create a model object containing the sampling and log likelihood functions, the data and parameter names.
model = DEModel(;
sample_prior,
prior_loglike,
loglike,
data,
names
)
Next, define the DifferentialEvolution sampling object. This requires bounds
, burnin
and Np
, which is the number of particles.
de = DE(;sample_prior, bounds, burnin = 1000, Np = 6)
To run the sampler, pass the model and differential evolution object along with settings for the number of iterations and MCMCMThreads() for multithreading.
n_iter = 2000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)
println(chains)