Usage

IPython notebook examples

Detailed examples are available in the project and online: A 2D gauss example and a toy model including a comparison to theoretical predictions.

General usage of abcpmc

In generall the use of abcpmc in simple:

import abcpmc
import numpy as np

#create "observed" data set
size = 5000
sigma = np.eye(4) * 0.25
means = np.array([1.1, 1.5, 1.1, 1.5])
data = np.random.multivariate_normal(means, sigma, size)
#-------

#distance function: sum of abs mean differences
def dist(x, y):
    return np.sum(np.abs(np.mean(x, axis=0) - np.mean(y, axis=0)))

#our "model", a gaussian with varying means
def postfn(theta):
    return np.random.multivariate_normal(theta, sigma, size)

eps = abcpmc.LinearEps(20, 5, 0.075)
prior = abcpmc.GaussianPrior(means*1.1, sigma*2) #our best guess

sampler = abcpmc.Sampler(N=10, Y=data, postfn=postfn, dist=dist)

for pool in sampler.sample(prior, eps):
    print("T: {0}, eps: {1:>.4f}, ratio: {2:>.4f}".format(pool.t, pool.eps, pool.ratio))
    for i, (mean, std) in enumerate(zip(np.mean(pool.thetas, axis=0), np.std(pool.thetas, axis=0))):
        print(u"    theta[{0}]: {1:>.4f} \u00B1 {2:>.4f}".format(i, mean,std))

the resulting output would look something like this:

T: 0, eps: 2.0000, ratio: 0.2439
    theta[0]: 0.9903 ± 0.5435
    theta[1]: 1.6050 ± 0.4912
    theta[2]: 1.0567 ± 0.4548
    theta[3]: 1.2859 ± 0.5213
T: 1, eps: 1.8987, ratio: 0.3226
    theta[0]: 1.1666 ± 0.4129
    theta[1]: 1.6597 ± 0.5227
    theta[2]: 1.1263 ± 0.3366
    theta[3]: 1.4711 ± 0.2150
T: 2, eps: 1.7974, ratio: 0.3030
    theta[0]: 1.1263 ± 0.2505
    theta[1]: 1.4832 ± 0.5057
    theta[2]: 1.0585 ± 0.3387
    theta[3]: 1.4782 ± 0.2808
T: 3, eps: 1.6961, ratio: 0.4167
    theta[0]: 1.1265 ± 0.1845
    theta[1]: 1.2032 ± 0.4470
    theta[2]: 1.0248 ± 0.2074
    theta[3]: 1.4689 ± 0.4250

...

T: 19, eps: 0.0750, ratio: 0.0441
    theta[0]: 1.1108 ± 0.0172
    theta[1]: 1.4832 ± 0.0166
    theta[2]: 1.0895 ± 0.0202
    theta[3]: 1.5016 ± 0.0097

Parallelisation on cluster with MPI

abcpmc has an built-in support for massively parallelized sampling on a cluster using MPI.

To make use of this parallelization the abcpmc Sampler need to be initialized with an instance of the MpiPool:

import abcpmc
from abcpmc import mpi_util

...

mpi_pool = mpi_util.MpiPool()
sampler = abcpmc.Sampler(N, Y, postfn, dist, pool=mpi_pool) #pass the mpi_pool

if mpi_pool.isMaster(): print("Start sampling")

for pool in sampler.sample(prior, eps):

...

If the threshold is dynamically adapted the user has to make sure that the state is synchonized among all MPI task with a broadcast:

eps.eps = mpi_util.mpiBCast(new_threshold)

Finally, the job has to be launched as follows to run on N tasks in parallel (might depend on your system):

$ mpirun -np N python <your-abc-script.py>