Reference

AdaptiveSMC(; alpha = 0.9, epsilon = .1, ess_min = 0.5, M = 1,
              prior, K = 10^3, proposal = defaultproposal(prior))

Adaptive Sequential Monte Carlo structure with K particles, M calls of the model per parameter value, final epsilon, decrease parameter alpha and resampling threshold ess_min. See also PMC.

Example

using Distributions, LinearAlgebra
asmc = AdaptiveSMC(prior = MultivariateNormal(zeros(2), 4I))
model(theta) = vcat([.3*randn(2) .+ theta for _ in 1:2]...) # 2 i.i.d samples
data = [.4, -.3, .5, -.2]
run!(asmc, model, data)

mean(asmc)
using StatsPlots
p = hcat(particles(asmc)...); StatsPlots.scatter(p[1, :], p[2, :])

Del Moral, P., Doucet, A., and Jasra, A. (2012) An adaptive sequential Monte Carlo method for approximate Bayesian computation. Statistics and Computing, 22, 1009-1020

PMC(; epsilon, K, prior, ess_min = 0.5, proposal = defaultproposal(prior))

Population Monte Carlo ABC structure, with epsilon schedule, K particles, prior, resampling threshold ess_min and proposal distribution. The epsilon schedule should be an iterable, e.g. an array of numbers, EpsilonExponentialDecay or EpsilonLinearDecay. One dimensional problems should be defined in vectorized form ( see run! for an example).

Example

pmc = PMC(epsilon = EpsilonExponentialDecay(1, .01, .9), K = 20,
                 prior = MultivariateUniform([-1, -1], [1, 1]))
model(theta) = .1 * randn(2) .+ theta
data = [.3, -.2]
run!(pmc, model, data)

mean(pmc)
particles(pmc)
weights(pmc)

Beaumont, M. A., Cornuet, J., Marin, J., and Robert, C. P. (2009) Adaptive approximate Bayesian computation. Biometrika,96, 983-990.

run!([rng], method, model, data;
     callback = () -> nothing, maxfevals = Inf, verbose = true)

Run method on model and data. The function callback gets evaluated after every iteration. The method stops after the first iteration that reaches more than maxfevals calls to the model.

Example

using Distributions, LinearAlgebra
pmc = PMC(epsilon = [1, .5, .2, .1, .01, .001], K = 10^3,
                 prior = TruncatedMultivariateNormal([1.], 2I,
                                                     lower = [0], upper = [3]))
model(theta) = theta * randn() .+ 1.2
data = [.5]
callback() = @show mean(pmc)
run!(pmc, model, data, callback = callback, maxfevals = 10^5)

using StatsBase, StatsPlots
h = fit(Histogram, vcat(particles(pmc)...), nbins = 50)
StatsPlots.plot(h)

K2ABC(; kernel, epsilon, prior, K)

Creates a K2ABC structure.

Example

using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = K2ABC(kernel = Kernel(),
          epsilon = .1,
          prior = MultivariateNormal([0.], [5.]),
          K = 10^3)
result = run!(k, model, data)
mean(k)

Park, M., Jitkrittum, W., and Sejdinovic, D. (2016), K2-ABC: Approximate Bayesian Computation with Kernel Embeddings, Proceedings of Machine Learning Research, 51:398-407

KernelABC(; kernel, prior, K, delta)

Creates a KernelABC structure.

Example

using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = KernelABC(kernel = Kernel(),
              delta = 1e-8,
              prior = MultivariateNormal([0.], [5.]),
              K = 10^3)
result = run!(k, model, data)
mean(k)

Fukumizu, K., Song, L. and Gretton, A. (2013), Kernel Bayes' Rule: Bayesian Inference with Positive Definite Kernels, Journal of Machine Learning Research, 14:3753-3783, http://jmlr.org/papers/v14/fukumizu13a.html

See also Muandet, K., Fukumizu, K., Sriperumbudur, B. and Schölkopf, B. (2017), Kernel Mean Embedding of Distributions: A Review and Beyond, Foundations and Trends® in Machine Learning, 10:1–141

KernelLoss(; loss = L, K = 50,
             kerneltype = ModifiedGaussian,
             gamma = .1,
             heuristic = Nothing,
             prior = nothing,
             distance = euclidean,
             kernel = Kernel(type = kerneltype,
                             heuristic = heuristic,
                             gamma = gamma,
                             distance = distance))

Constructs an structure, used to compare K samples generated from a model to some data with a kernel and loss function loss. The loss function has the signature L(list_of_samples, data, kernel); by default it is the negative log of the average distance between samples and data (see LikelihoodfreeInference.L). The KernelLoss can be used as a loss in a PointEstimator.

The returned structure is callable.

Example

model(x) = randn() .+ x
data = [2.]
k = KernelLoss()
k([1.], model, data)

K can be an integer or a function that returns integers, to implement a schedule.

Example

using Statistics
model(x) = randn() .+ x
data = [2.]
k = KernelLoss(K = (n) -> n > 20 ? 10^3 : 10)
std(k([1.], model, data, 1) for _ in 1:100)   # large std, because K = 10
std(k([1.], model, data, 21) for _ in 1:100)  # small std, because K = 10^3

This KernelLoss is inspired by Bertl, J., Ewing, G., Kosiol, C., and Futschik, A., (2017) Approximate maximum likelihood estimation for population genetic inference, Statistical Applications in Genetics and Molecular Biology. 16:5-6

PointEstimator(optimizer, loss)

PointEstimator(; prior = nothing, lower = nothing, upper = nothing,
                 optimizer = nlopt(prior),
                 losstype = QDLoss, kwargs...)

Creates a PointEstimator with given prior, optimizer and loss = losstype(; prior = prior, kwargs...). If the prior is nothing a uniform prior is assumed. In this case the optimizer has to be initialized with bounds, , e.g. nlopt(lower, upper) where lower and upper are arrays.

Example

model(x) = randn() .+ x
data = [2.]
p = PointEstimator(lower = [-5], upper = [5], losstype = QDLoss, K = 10^3)
res = run!(p, model, data, maxfevals = 10^5)

QDLoss(; K = 50, epsilon = 2/K, prior = nothing, distance = euclidean)

Constructs an structure, used to compare K samples generated from a model to some data by computing the first epsilon-quantile of the distance between samples and data. The QDLoss can be used as a loss in a PointEstimator.

The returned structure is callable.

Example

model(x) = randn() .+ x
data = [2.]
q = QDLoss()
q([1.], model, data)

K can be an integer or a function that returns integers, to implement a schedule. Similarly, epsilon can be a number or a function that returns numbers.

Example

using Statistics
model(x) = randn() .+ x
data = [2.]
q = QDLoss(K = (n) -> n > 20 ? 10^3 : 10)
std(q([1.], model, data, 1) for _ in 1:100)   # large std, because K = 10
std(q([1.], model, data, 21) for _ in 1:100)  # small std, because K = 10^3

KernelRecursiveABC(; kernel, prior, K, delta = 1e-2,
                     herding_options = (maxtime = 10,
                                        method = :LD_LBFGS,
                                        restarts = 10),
                     kernelx = Kernel())

Creates a KernelRecursiveABC structure. KernelRecursiveABC iterates a KernelABC step and a kernel herding step. kernel, prior, K and delta determine the KernelABC step. herding_options options and kernelx determine the herding step.

Example

using Distributions
model(x) = [randn() .+ x for _ in 1:3]
data = [[3.], [2.9], [3.3]]
k = KernelRecursiveABC(kernel = Kernel(),
                       delta = 1e-8,
                       prior = MultivariateNormal([0.], [5.]),
                       K = 10^2)
result = run!(k, model, data, maxfevals = 10^3)
result.x

Kajihara, T., Kanagawa, M., Yamazaki, K. and Fukumizu, K. (2018), Kernel Recursive ABC: Point Estimation with Intractable Likelihood, Proceedings of the 35th International Conference on Machine Learning, 80:2400-2409

L(x, y, k) = -log(mean(k(xi, y) for xi in x))

MultiParticleNormal(particles, weights;
                    lower = nothing, upper = nothing,
                    diagonal = false,
                    scaling = 2, regularization = 0.)

Mixture of k multivariate Gaussians with means given by an array of particles (length(particles) == k), weighted by weights and with covariance matrix given by scaling * cov(particles) + regularization * I for each Gaussian, if diagonal == false; otherwise scaling * std(particles) + regularization. Sampling truncated Gaussians with diagonal covariance matrix can be more efficient. The mean of the distribution is the weighted average of the particles.

If lower and upper are not nothing the Gaussians are truncated to the box confined by the lower and upper bounds. Note that the logpdf of the truncated version is not properly normalized (see also TruncatedMultivariateNormal).

Example

julia> d = MultiParticleNormal([[0, 0], [2, 3], [4, 0]], [1/3, 1/3, 1/3]);

julia> mean(d)
2-element Array{Float64,1}:
 2.0
 1.0

MultiParticleNormal(; kwargs...)

Returns a anonymous function that accepts particles as input and creates a MultiParticleNormal distributions with uniform weights. It accepts the same keyword arguments kwargs... as the normal MultiParticleNormal constructor.

Example

julia> constructor = MultiParticleNormal(scaling = 3, lower = [-1, 0], upper = [5, 5]);

julia> d = constructor([[0, 0], [2, 3], [4, 0], [5, 1]])
MultiParticleNormal(4 particles in 2 dimensions, scaling = 3.0, regularization = 0.0, with bounds)

julia> d.weights
4-element Array{Float64,1}:
 0.25
 0.25
 0.25
 0.25

MultivariateUniform(lower, upper)

Uniform distribution in d dimensions, where d == length(lower) == length(upper).

Example

julia> d = MultivariateUniform([-2., -3], [4, 1])
MultivariateUniform{Float64} in 2 dimensions

julia> logpdf(d, [-3, 0])
-Inf

julia> logpdf(d, [-1, 0])
-3.1780538303479453

TruncatedMultivariateNormal(mvnormal, lower, upper)

Multivariate normal distribution mvnormal truncated to a box given by lower and upper bounds. Simple rejection sampling is implemented. IMPORTANT: logpdf is not properly normalized.

Example

julia> using Distributions

julia> d = TruncatedMultivariateNormal(MultivariateNormal([2., 3.], .3*I), [0, 0], [4, 5]);

julia> logpdf(d, [2, 3])
-0.6339042620834094

julia> logpdf(d, [2, 3]) == logpdf(d.mvnormal, [2, 3])
true

julia> logpdf(d, [-1, 1])
-Inf

TruncatedMultivariateNormal(m, cov; lower, upper)

Constructs a TruncatedMultivariateNormal with mean m covariance matrix cov and bounds lower and upper.

Example

julia> d = TruncatedMultivariateNormal([0, 0], 3I;
                                       lower = [-1, -4], upper = [4, 5])

bbo(lower, upper; kwargs...)
bbo(prior; kwargs...) = bbo(lower(prior), upper(prior); kwargs...)

BlackBoxOptim optimizer for PointEstimator, available when using BlackBoxOptim. kwargs are passed to BlackBoxOptim.bboptimize.

bo(lower, upper;
        acquisition = BO.UpperConfidenceBound(),
        capacity = 10^3,
        mean = BO.GaussianProcesses.MeanConst(0.),
        kernel = BO.GaussianProcesses.SEArd(zeros(length(lower)), 5.),
        logNoise = 0,
        modeloptimizer_options = NamedTuple(),
        kwargs...)
bo(prior; kwargs...) = bo(lower(prior), upper(prior); kwargs...)

BayesianOptimization optimizer for PointEstimator, available when using BayesianOptimization. kwargs are passed to BayesianOptimization.BOpt.

cma(lower, upper;
    sigma0 = 0.25,
    fminoptions = NamedTuple(),
    filename = "tmp-",
    kwargs...)
 cma(prior; kwargs...) = cma(lower(prior), upper(prior); kwargs...)

CMA optimizer for PointEstimator, available when using PyCMA. kwargs are passed to PyCMA.cma.CMAOptions, fminoptions to PyCMA.cma.fmin.

spsa(lower, upper; kwargs...)
spsa(prior; kwargs...) = spsa(lower(prior), upper(prior); kwargs...)

SimultaneousPerturbationStochasticApproximation optimizer for PointEstimator, available when using SimultaneousPerturbationStochasticApproximation. kwargs are passed to SimultaneousPerturbationStochasticApproximation.SPSA.

EpsilonExponentialDecay(init, last, decay)

Exponentially decreasing sequence starting at init, decaying with factor decay, ending as soon as it equals (or drops below) last.

Example

julia> collect(EpsilonExponentialDecay(1, .2, .5))
4-element Array{Any,1}:
 1.0
 0.5
 0.25
 0.125

EpsilonLinearDecay(init, last, decay)

Linearly decreasing sequence starting at init, decaying with step size decay, ending as soon as it equals (or drops below) last.

Example

julia> collect(EpsilonLinearDecay(10, 5, 2))
4-element Array{Any,1}:
 10.0
  8.0
  6.0
  4.0

defaultproposal(prior)

Returns a function that takes a list of particles and returns a proposal distribution. This function is called during initialization of PMC and AdaptiveSMC.