DEMetropolis Documentation

Run the Differential Evolution Markov Chain (DE-MC) sampler proposed by ter Braak (2006)

This sampler uses the de_update step. It can optionally switch between two γ values (γ₁ and γ₂) with probability p_γ₂. This is so that the sampler can occasionally move between modes, by having γ₂ = 1 while γ₁ remains the optimal value based on the dimension of the problem.

This algorithm varies slightly from the original. Updates within a population occur on the previous position of that population. i.e. if chain 1 has been updated (a₁ → a₂) and chain 2 picks chain 1 to update from, then the value of chain 1 used by chain 2 is the pre-update version of chain 1 (a₁). This change allows the algorithm to be easily parallelised.

See doi.org/10.1007/s11222-006-8769-1 for more information on sampler.

Arguments

• ld: The log-density function to sample from, intended to be a LogDensityProblem.
• n_its: The number of sampling iterations per chain.

Keyword Arguments

• n_burnin: Number of burn-in iterations. Defaults to n_its * 5.
• n_chains: Number of chains. Defaults to dimension(ld) * 2.
• initial_state: Initial states for the chains. Defaults to randn(rng, n_chains, dimension(ld)).
• N₀: Size of the initial population (must be >= nchains + 3), only used if using memory-based sampling. Defaults to `nchains * 2`.
• memory: Use memory-based sampling (true) or memoryless (false). Defaults to false.
• save_burnt: Save burn-in samples. Defaults to false.
• parallel: Run chains in parallel. Defaults to false.
• rng: Random number generator. Defaults to default_rng().
• diagnostic_checks: Diagnostic checks to run during burn-in. Defaults to nothing.
• check_epochs: Splits n_burnin into check_epochs + 1 epochs and applies the diagnostic checks at the end of each epoch, other than the final epoch. Defaults to 1.
• thin: Thinning interval. Defaults to 1.
• γ₁: Primary scaling factor for DE update. Defaults to 2.38 / sqrt(2 * dim).
• γ₂: Secondary scaling factor for DE update. Defaults to 1.0.
• p_γ₂: Probability of using γ₂. Defaults to 0.1.
• β: Noise distribution for DE update. Defaults to Uniform(-1e-4, 1e-4).

Returns

• A named tuple containing the samples, sampler scheme, and potentially burn-in samples.

Example

julia> deMC(ld, 1000; n_chains = 10)

See also composite_sampler, deMCzs, DREAMz.

Run the Differential Evolution Markov Chain with snooker update and historic sampling (DE-MCzs) sampler proposed by ter Braak and Vrugt (2008).

This sampler runs until a stopping_criteria (default: Rhat of the last 50% of the chains is <1.2) is met. The sampler can occasionally propose snooker updates which can sample areas away from the current chains. Combined with the adaptive memory-based sampling this sampler can efficiently sample from a problem where n_chains < the dimension of the problem.

See: doi.org/10.1007/s11222-008-9104-9 for more information

Arguments

• ld: The log-density function to sample from, intended to be a LogDensityProblem.
• epoch_size: The number of saved iterations per chain per epoch.

Keyword Arguments

• warmup_epochs: Number of warm-up epochs. Defaults to 5.
• epoch_limit: Maximum number of sampling epochs. Defaults to 20.
• n_chains: Number of chains. Defaults to dimension(ld) * 2.
• N₀: Size of the initial population (must be >= nchains + 3). Defaults to `nchains * 2`.
• initial_state: Initial population state. Defaults to randn(rng, N₀, dimension(ld)).
• memory: Use memory based sampling? Defaults to true.
• save_burnt: Save warm-up samples. Defaults to true.
• parallel: Run chains in parallel with multithreading. Defaults to false.
• rng: Random number generator. Defaults to default_rng().
• diagnostic_checks: Diagnostic checks during warm-up. Defaults to nothing.
• stopping_criteria: Criterion to stop sampling. Defaults to R̂_stopping_criteria().
• γ: Scaling factor for the DE update. Defaults to 2.38 / sqrt(2 * dim).
• γₛ: Scaling factor for the Snooker update. Defaults to 2.38 / sqrt(2).
• p_snooker: Probability of using the Snooker update. Defaults to 0.1.
• β: Noise distribution for DE update. Defaults to Uniform(-1e-4, 1e-4).
• thin: Thinning interval. Defaults to 10.

Returns

• A named tuple containing the samples, sampler scheme, and potentially burn-in samples.

Example

julia> deMCzs(ld, 1000; n_chains = 3)

See also composite_sampler, deMC, DREAMz.

Run the Differential Evolution Adaptive Metropolis (DREAMz) sampler

This sampler runs until a stopping_criteria (default: Rhat of the last 50% of the chains is <1.2) is met. The sampler uses subspace_sampling, where the cross-over probability is adapted during the burn-in period. It can optionally switch between two γ values, so that γ₁ can be the optimal value (based on sampled parameters) and γ₂ can be some fixed value (i.e. 1) so that the sampler can switch modes. This sampler also checks for outlier chains (where the mean log-density falls outside the IQR) and replaces then with the position of the chain with the highest log-density. This step breaks detailed balance its not performed in the last epoch of the warm-up period.

By default this is a memory-based sampler (DREAMz). Setting memory = false makes this the DREAM sampler.

See doi.org/10.1515/IJNSNS.2009.10.3.273 for more info.

Arguments

• ld: The log-density function to sample from, intended to be a LogDensityProblem.
• epoch_size: The number of saved iterations per chain per epoch.

Keyword Arguments

• warmup_epochs: Number of warm-up epochs. Defaults to 5. Crossover probabilities are adapted in this period.
• epoch_limit: Maximum number of sampling epochs. Defaults to 20.
• n_chains: Number of chains. Defaults to dimension(ld) * 2.
• N₀: Size of the initial population (must be >= nchains). Defaults to `nchains. Only the firstn_chainswill be used ifmemory = false`.
• initial_state: Initial population state. Defaults to randn(rng, N₀, dimension(ld)).
• memory: Use memory-based sampling (true) or memoryless (false). Defaults to true.
• save_burnt: Save warm-up samples. Defaults to true.
• parallel: Run chains in parallel. Defaults to false.
• rng: Random number generator. Defaults to default_rng().
• diagnostic_checks: Diagnostic checks during warm-up. Defaults to [ld_check()].
• stopping_criteria: Criterion to stop sampling. Defaults to R̂_stopping_criteria().
• γ₁: Primary scaling factor for subspace update. Defaults to nothing (uses 2.38 / sqrt(2 * δ * d)). Can also be a Real value.
• γ₂: Secondary scaling factor for subspace update. Defaults to 1.0. Can also be a Real value
• p_γ₂: Probability of using γ₂. Defaults to 0.2.
• n_cr: Number of crossover probabilities to adapt if cr₁/cr₂ are nothing. Defaults to 3.
• cr₁: Crossover probability distribution/value for γ₁. Defaults to nothing (adaptive). Can also be a Real value (<1) or a Distributions.UnivariateDistribution, in either case it is not adapted.
• cr₂: Crossover probability distribution/value for γ₂. See above.
• ϵ: Additive noise distribution. Defaults to Uniform(-1e-4, 1e-4).
• e: Multiplicative noise distribution. Defaults to Normal(0.0, 1e-2).
• δ: Number of difference vectors distribution. Defaults to DiscreteUniform(1, 3).
• thin: Thinning interval. Defaults to 1.

Returns

• A named tuple containing the samples, sampler scheme, and potentially burn-in samples.

Example

julia> DREAMz(ld, 1000; n_chains = 10)

See also composite_sampler, deMC, deMCzs.

Create a sampler scheme defining the update steps to be used in composite_sampler.

The update used in each iteration for each chain is randomly selected from the updates given here.

Arguments

• updates...: One or more update_struct objects (e.g., created by setup_de_update, setup_snooker_update, setup_subspace_sampling or your own custom sampler).

Keyword Arguments

• w: A vector of weights corresponding to each update step. If nothing, updates are chosen with equal probability. Weights must be non-negative.

Examples

# only snooker updates
julia> setup_sampler_scheme(setup_snooker_update())

# DE and Snooker
julia> setup_sampler_scheme(setup_snooker_update(), setup_de_update())

# With weights, snookers 10% of the time
julia> setup_sampler_scheme(setup_snooker_update(), setup_de_update(), w = [0.9, 0.1])

Run a composite MCMC sampler for a fixed number of iterations.

For sampling with your own sampling scheme from setup_sampling_scheme

Arguments

• ld: The log-density function to sample from, intended to be a LogDensityProblem.
• n_its: The desired number of samples per chain.
• n_chains: The number of chains to run.
• memory: Boolean indicating whether to sample from historic values instead sampling from current chains.
• initial_state: An array containing the initial states for the chains and the initial population if memory=true.
• sampler_scheme: A sampler_scheme_struct defining the update steps to use.

Keyword Arguments

• thin: Thinning interval for storing samples. Defaults to 1. If using memory this also effects which samples are added to the memory.
• save_burnt: Boolean indicating whether to save burn-in samples. Defaults to false. Does not save thinned samples
• rng: Random number generator. Defaults to default_rng().
• n_burnin: Number of burn-in iterations. Defaults to n_its * 5. Any adaptions occur over this period.
• parallel: Boolean indicating whether to run chains in parallel using multithreading. Defaults to false.
• diagnostic_checks: A vector of diagnostic_check_struct to run during burn-in. Defaults to nothing.
• check_epochs: Splits n_burnin into check_epochs + 1 epochs and applies the diagnostic checks at the end of each epoch, other than the final epoch. Defaults to 1.

Returns

• A named tuple containing the samples, sampler scheme, and potentially burn-in samples.

Example

# DE and Snooker sample scheme with memory
julia> composite_sampler(
    ld, 1000, 10, true, initial_state, setup_sampler_scheme(setup_snooker_update(), setup_de_update())
)

See also deMC, deMCzs, DREAMz.

Run a composite MCMC sampler until a stopping criterion is met.

For sampling with your own sampling scheme from setup_sampling_scheme.

Arguments

• ld: The log-density function to sample from, intended to be a LogDensityProblem.
• epoch_size: The number of saved iterations per chain per epoch.
• n_chains: The number of chains to run.
• memory: Boolean indicating whether to sample from historic values instead sampling from current chains.
• initial_state: An array containing the initial states for the chains and the initial population if memory==true.
• sampler_scheme: A sampler_scheme_struct defining the update steps to use.
• stopping_criteria: A stopping_criteria_struct defining when to stop sampling.

Keyword Arguments

• thin: Thinning interval for storing samples. Defaults to 1. If using memory this also effects which samples are added to the memory.
• save_burnt: Boolean indicating whether to save burn-in samples. Defaults to false. Does not save thinned samples.
• rng: Random number generator. Defaults to default_rng().
• warmup_epochs: Number of warm-up epochs before we begin checking the stopping criteria. Defaults to 5. Samples from these epochs won't be included in the final sample. Sampler adaptation only occurs in this period.
• parallel: Boolean indicating whether to run chains in parallel using multithreading. Defaults to false.
• epoch_limit: Maximum number of sampling epochs to run. Defaults to 20.
• diagnostic_checks: A vector of diagnostic_check_struct to run during warm-up. Defaults to nothing.

Returns

• A named tuple containing the samples, sampler scheme, and potentially burn-in samples.

Example

# DE and Snooker sample scheme with memory until Rhat ≤ 1.05
julia> composite_sampler(
    ld, 1000, 10, true, initial_state, setup_sampler_scheme(setup_snooker_update(), setup_de_update()), R̂_stopping_criteria(1.05)
)

See also deMC, deMCzs, DREAMz.

Set up a Differential Evolution (DE) update step.

See doi.org/10.1007/s11222-006-8769-1 for more information.

Arguments

• ld: The log-density function (used to determine dimension if γ is not provided).

Keyword Arguments

• γ: The scaling factor for the difference vector. Can be a Real, a Distributions.UnivariateDistribution, or nothing. If nothing, it defaults based on deterministic_γ and the problem dimension.
• β: Distribution for the small noise term added to the proposal. Defaults to Uniform(-1e-4, 1e-4).
• deterministic_γ: If true and γ is nothing, sets γ to the theoretically optimal 2.38 / sqrt(2 * dim). If false, sets γ to Uniform(0.8, 1.2). Defaults to true.

Example

julia> setup_de_update(ld; β = Normal(0.0, 0.01))

See also setup_snooker_update, setup_subspace_sampling.

Set up a Snooker update step.

See doi.org/10.1007/s11222-008-9104-9 for more information.

Keyword Arguments

• γ: The scaling factor for the projection. Can be a Real, a Distributions.UnivariateDistribution, or nothing. If nothing, it defaults based on deterministic_γ.
• deterministic_γ: If true and γ is nothing, sets γ to the theoretically optimal 2.38 / sqrt(2). If false, sets γ to Uniform(0.8, 1.2). Defaults to true.

Example

julia> setup_snooker_update(γ = Uniform(0.1, 2.0))

See also setup_de_update, setup_subspace_sampling.

Set up a Subspace Sampling (DREAM-like) update step.

See doi.org/10.1515/IJNSNS.2009.10.3.273 for more information.

Keyword Arguments

• γ: Fixed scaling factor for the difference vector sum. If nothing (default), uses 2.38 / sqrt(2 * δ * d) where d is the number of updated dimensions. If a Real is provided, uses that fixed value.
• cr: Crossover probability. Can be a Real (fixed probability), nothing (adaptive probability using n_cr values), or a Distributions.UnivariateDistribution. Defaults to nothing.
• n_cr: Number of crossover probabilities to adapt between if cr is nothing. Defaults to 3.
• δ: Number of difference vectors to add. Can be an Integer or a Distributions.DiscreteUnivariateDistribution. Defaults to DiscreteUniform(1, 3).
• ϵ: Distribution for small noise added to the proposal in the selected subspace. Defaults to Uniform(-1e-4, 1e-4).
• e: Distribution for multiplicative noise (e + 1) applied to the difference vector sum. Defaults to Normal(0.0, 1e-2).

Example

julia> setup_subspace_sampling(cr = Beta(1, 2))

See also setup_de_update, setup_snooker_update.

Create a stopping criterion based on the rank Gelman-Rubin diagnostic (R̂). Sampling stops when the R̂ value for all parameters is below maximum_R̂.

Arguments

• maximum_R̂: The maximum acceptable R̂ value. Defaults to 1.2.

See also MCMCDiagnosticTools.rhat.

Create a diagnostic check that identifies outlier chains based on their mean log-density (of the last 50% of the chain) during burn-in/warm-up. Outlier chains (those with mean log-density below Q1 - 2*IQR) are reset to the state of the chain with the highest mean log-density.

See: Vrugt 2009 doi.org/10.1515/IJNSNS.2009.10.3.273

See also acceptance_check.

Create a diagnostic check that identifies poorly mixing chains based on their acceptance rate (of the last 50% of the chain) during burn-in/warm-up. Chains with acceptance rates below min_acceptance and significantly lower than others (based on log-acceptance rate IQR) are reset to the state of the chain closest to the target_acceptance rate.

Arguments

• min_acceptance: The minimum acceptable acceptance rate. Defaults to 0.1.
• target_acceptance: The target acceptance rate used to select the best chain for resetting outliers. Defaults to 0.24.

See also ld_check.