KissABC 3.0

ApproxKernelizedPosterior(
    prior::Distribution,
    cost::Function,
    target_average_cost::Real
)

this function will return a type which can be used in the sample function as an ABC density, this type works by assuming Gaussianly distributed errors 𝒩(0,ϵ), ϵ is specified in the variable target_average_cost.

ApproxPosterior(
    prior::Distribution,
    cost::Function,
    max_cost::Real
)

this function will return a type which can be used in the sample function as an ABC density, this type works by assuming uniformly distributed errors in [-ϵ,ϵ], ϵ is specified in the variable max_cost.

CommonLogDensity(nparameters, sample_init, lπ)

this function will return a type for performing classical MCMC via the sample function.

nparameters: total number of parameters per sample.

sample_init: function which accepts an RNG::AbstractRNG and returns a sample for lπ.

lπ: function which accepts a sample, and returns a log-density float value.

Factored{N} <: Distribution{Multivariate, MixedSupport}

a Distribution type that can be used to combine multiple UnivariateDistribution's and sample from them. Example: it can be used as prior = Factored(Normal(0,1), Uniform(-1,1))

length(p::Factored) = begin

returns the number of distributions contained in p.

rand(rng::AbstractRNG, factoreddist::Factored)

function to sample one element from a Factored object

logpdf(d::Factored, x) = begin

Function to evaluate the logpdf of a Factored distribution object

pdf(d::Factored, x) = begin

Function to evaluate the pdf of a Factored distribution object

Inverse cdf of g-pdf, see eq. 10 of Foreman-Mackey et al. 2013.

Sample from g using inverse transform sampling. a=2.0 is recommended.

Adaptive SMC from P. Del Moral 2012, with Affine invariant proposal mechanism, faster that AIS for ABC targets.

function smc(
    prior::Distribution,
    cost::Function;
    rng::AbstractRNG = Random.GLOBAL_RNG,
    nparticles::Int = 100,
    M::Int = 1,
    alpha = 0.95,
    mcmc_retrys::Int = 0,
    mcmc_tol = 0.015,
    epstol = 0.0,
    r_epstol = (1 - alpha)^1.5 / 50,
    min_r_ess = alpha^2,
    max_stretch = 2.0,
    verbose::Bool = false,
    parallel::Bool = false,
)

• prior: a Distribution object representing the parameters prior.
• cost: a function that given a prior sample returns the cost for said sample (e.g. a distance between simulated data and target data).
• rng: an AbstractRNG object which is used by SMC for inference (it can be useful to make an inference reproducible).
• nparticles: number of total particles to use for inference.
• M: number of cost evaluations per particle, increasing this can reduce the chance of rejecting good particles.
• alpha - used for adaptive tolerance, by solving ESS(n,ϵ(n)) = α ESS(n-1, ϵ(n-1)) for ϵ(n) at step n.
• mcmc_retrys - if set > 0, whenever the fraction of accepted particles drops below the tolerance mcmc_tol the MCMC step is repeated (no more than mcmc_retrys times).
• mcmc_tol - stopping condition for SMC, if the fraction of accepted particles drops below mcmc_tol the algorithm terminates.
• epstol - stopping condition for SMC, if the adaptive cost threshold drops below epstol the algorithm has converged and thus it terminates.
• min_r_ess - whenever the fractional effective sample size drops below min_r_ess, a systematic resampling step is performed.
• max_stretch - the proposal distribution of smc is the stretch move of Foreman-Mackey et al. 2013, the larger the parameters the wider becomes the proposal distribution.
• verbose - if set to true, enables verbosity.
• parallel - if set to true, threaded parallelism is enabled, keep in mind that the cost function must be Thread-safe in such case.

Example

using KissABC
prior=Factored(Normal(0,5), Normal(0,5))
cost((x,y)) = 50*(x+randn()*0.01-y^2)^2+(y-1+randn()*0.01)^2
results = smc(prior, cost, alpha=0.5, nparticles=5000).P

output:

2-element Array{Particles{Float64,5000},1}:
 1.0 ± 0.029
 0.999 ± 0.012

sample(model, AIS(N), Ns[; optional args])
sample(model, AIS(N), MCMCThreads(), Ns, Nc[; optional keyword args])
sample(model, AIS(N), MCMCDistributed(), Ns, Nc[; optional keyword args])

Generalities

This function will run an Affine Invariant MCMC sampler, and will return a Particles object for each parameter, the mandatory parameters are:

model: a subtype of AbstractDensity, look at ApproxPosterior, ApproxKernelizedPosterior, CommonLogDensity.

N: number of particles in the ensemble, this particles will be evolved to generate new samples.

Ns: total number of samples which must be recorded.

Nc: total number of chains to run in parallel if MCMCThreads or MCMCDistributed is enabled.

the optional arguments available are:

discard_initial: number of mcmc particles to discard before saving any sample.

ntransitions: number of mcmc steps per particle between each sample.

retry_sampling: number of maximum attempts to resample an initial particle whose cost (or log-density) is ±∞ or NaN.

progress: a boolean to disable verbosity

Minimal Example for CommonLogDensity:

using KissABC
D = CommonLogDensity(
    2, #number of parameters
    rng -> randn(rng, 2), # initial sampling strategy
    x -> -100 * (x[1] - x[2]^2)^2 - (x[2] - 1)^2, # rosenbrock banana log-density
)
res = sample(D, AIS(50), 1000, ntransitions = 100, discard_initial = 500, progress = false)
println(res)

output:

Particles{Float64,1000}[1.43 ± 1.4, 0.99 ± 0.67]

Minimal Example for ApproxKernelizedPosterior (ApproxPosterior)

using KissABC
prior = Uniform(-10, 10) # prior distribution for parameter
sim(μ) = μ + rand((randn() * 0.1, randn())) # simulator function
cost(x) = abs(sim(x) - 0.0) # cost function to compare simulations to target data, in this case simply '0'
plan = ApproxPosterior(prior, cost, 0.01) # Approximate model of log-posterior density (ABC)
#                                           ApproxKernelizedPosterior can be used in the same fashion
res = sample(plan, AIS(100), 2000, discard_initial = 10000, progress = false)
println(res)

output:

0.0 ± 0.46