MonAPMC

Parallel ABC algorithm based on APMC by Lenormand, Jabot and Deffuant (2012). MonAPMC stands for "Monoid APMC".

Given a stochastic function f: Vector[Double] => Vector[Double], and a value y: Vector[Double], the algorithm aims at finding what input vectors xs are likely to have resulted in the vector y through f. The objective is to estimate the probability distribution P(X|Y = y), where X represent f's input values, and Y its output values. The algorithm returns a sample of vectors that are distributed according to this distribution.

The simplest way to run the algorithm is to use the MonAPMC.run and MonAPMC.scan functions. The first returns the final algorithm state, and the second the sequence of states corresponding to the successive steps of the alorithm. The posterior sample is in the final state s: s.thetas contains the particles, and s.weights their weights. For example, the expected value of the particle according to the posterior is the weighted average of the particles.

Controlling the function f evaluation: During its course, the algorithm evaluates the user-provided function f many times. It can be useful to gain control of the function evaluation. For this purpose, use the ExposedEval object returned by MonAPMC.exposedStep. Its members, ExposedEval.pre and ExposedEval.post decompose the computation around the function evaluation. ExposedEval.pre returns an intermediate state and a org.apache.commons.math3.linear.RealMatrix: (pass, xs). The rows of the matrix xs are the input vectors with which to evaluate the function f. It is up to you to use these values to construct a new matrix, ys, such that its rows are the output values of f with the corresponding row in xs. The row order must be kept. Then, give (pass, ys) back to ExposedEval.post to get the new algorithm state, and reiterate. For example, to run the algorithm sequentially until it stops:

def f(x: Array[Double]): Array[Double] = ...

val step = MonAPMC.exposedStep(p)

var state = MonAPMC.Empty()
while (!MonAPMC.stop(p, state) {
  (pass, xs) = step.pre(state)
  ys = f(xs)
  state = step.post(pass, ys)
}

To run the algorithm in parallel, use the functions MonAPMC.split to split a current algorithm state into 2, such that the algorithm can be continued in parallel. You can step over the two in parallel for as many steps as you like, and merge them back together with MonAPMC.append.

val s0 = MonAPMC.Empty()
val (s0a, s0b) = MonAPMC.split(s0)
val s1a = MonAPMC.step(s0a)
val s2b = MonAPMC.step(MonAPMC.step(s0b))
val s3 = MonAPMC.append(s1a, s2b)

You can split states recursively as many times as you like to run more than 2 parallel threads:

var (s1,s2) = MonAPMC.split(MonAPMC.Empty())
var (s3,s4) = MonAPMC.split(s2)
/* step over all 4 states */
s = MonAPMC.append(s1,MonAPMC.append(s2,MonAPMC.append(s3, s4))),