This tutorial is generated from a Jupyter notebook that can be found here.

# ELFI tutorial¶

This tutorial covers the basics of using ELFI, i.e. how to make models, save results for later use and run different inference algorithms.

Let’s begin by importing libraries that we will use and specify some settings.

```
import time
import numpy as np
import scipy.stats
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
%precision 2
# Set an arbitrary seed and a global random state to keep the randomly generated quantities the same between runs
seed = 20170530
np.random.seed(seed)
```

## Inference with ELFI: case MA(2) model¶

Throughout this tutorial we will use the 2nd order moving average model MA(2) as an example. MA(2) is a common model used in univariate time series analysis. Assuming zero mean it can be written as

where \(\theta_1, \theta_2 \in \mathbb{R}\) and \((w_k)_{k\in \mathbb{Z}} \sim N(0,1)\) represents an independent and identically distributed sequence of white noise.

### The observed data and the inference problem¶

In this tutorial, our task is to infer the parameters \(\theta_1, \theta_2\) given a sequence of 100 observations \(y\) that originate from an MA(2) process. Let’s define the MA(2) simulator as a Python function:

```
def MA2(t1, t2, n_obs=100, batch_size=1, random_state=None):
# Make inputs 2d arrays for numpy broadcasting with w
t1 = np.asanyarray(t1).reshape((-1, 1))
t2 = np.asanyarray(t2).reshape((-1, 1))
random_state = random_state or np.random
w = random_state.randn(batch_size, n_obs+2) # i.i.d. sequence ~ N(0,1)
x = w[:, 2:] + t1*w[:, 1:-1] + t2*w[:, :-2]
return x
```

Above, `t1`

, `t2`

, and `n_obs`

are the arguments specific to the
MA2 process. The latter two, `batch_size`

and `random_state`

are
ELFI specific keyword arguments. The `batch_size`

argument tells how
many simulations are needed. The `random_state`

argument is for
generating random quantities in your simulator. It is a
`numpy.RandomState`

object that has all the same methods as
`numpy.random`

module has. It is used for ensuring consistent results
and handling random number generation in parallel settings.

### Vectorization¶

What is the purpose of the `batch_size`

argument? In ELFI, operations
are vectorized, meaning that instead of simulating a single MA2 sequence
at a time, we simulate a batch of them. A vectorized function takes
vectors as inputs, and computes the output for each element in the
vector. Vectorization is a way to make operations efficient in Python.
Above we rely on numpy to carry out the vectorized calculations.

In this case the arguments `t1`

and `t2`

are going to be vectors of
length `batch_size`

and the method returns a 2d array with the
simulations on the rows. Notice that for convenience, the funtion also
works with scalars that are first converted to vectors.

Note

There is a built-in tool (elfi.tools.vectorize) in ELFI to vectorize operations that are not vectorized. It is basically a for loop wrapper.

Important

In order to guarantee a consistent state of pseudo-random number generation, the simulator must have random_state as a keyword argument for reading in a numpy.RandomState object.

Let’s now use this simulator to create toy observations. We will use parameter values \(\theta_1=0.6, \theta_2=0.2\) as in *Marin et al. (2012)* and then try to infer these parameter values back based on the toy observed data alone.

```
# true parameters
t1_true = 0.6
t2_true = 0.2
y_obs = MA2(t1_true, t2_true)
# Plot the observed sequence
plt.figure(figsize=(11, 6));
plt.plot(y_obs.ravel());
# To illustrate the stochasticity, let's plot a couple of more observations with the same true parameters:
plt.plot(MA2(t1_true, t2_true).ravel());
plt.plot(MA2(t1_true, t2_true).ravel());
```

### Approximate Bayesian Computation¶

Standard statistical inference methods rely on the use of the
*likelihood* function. Given a configuration of the parameters, the
likelihood function quantifies how likely it is that values of the
parameters produced the observed data. In our simple example case above
however, evaluating the likelihood is difficult due to the unobserved
latent sequence (variable `w`

in the simulator code). In many real
world applications the likelihood function is not available or it is too
expensive to evaluate preventing the use of traditional inference
methods.

One way to approach this problem is to use Approximate Bayesian Computation (ABC) which is a statistically based method replacing the use of the likelihood function with a simulator of the data. Loosely speaking, it is based on the intuition that similar data is likely to have been produced by similar parameters. Looking at the picture above, in essence we would keep simulating until we have found enough sequences that are similar to the observed sequence. Although the idea may appear inapplicable for the task at hand, you will soon see that it does work. For more information about ABC, please see e.g.

- Lintusaari, J., Gutmann, M. U., Dutta, R., Kaski, S., and Corander, J. (2016). Fundamentals and recent developments in approximate Bayesian computation. *Systematic Biology*, doi: 10.1093/sysbio/syw077.
- Marin, J.-M., Pudlo, P., Robert, C. P., and Ryder, R. J. (2012). Approximate Bayesian computational methods. *Statistics and Computing*, 22(6):1167–1180.
- https://en.wikipedia.org/wiki/Approximate_Bayesian_computation

## Defining the model¶

ELFI includes an easy to use generative modeling syntax, where the generative model is specified as a directed acyclic graph (DAG). This provides an intuitive means to describe rather complex dependencies conveniently. Often the target of the generative model is a distance between the simulated and observed data. To start creating our model, we will first import ELFI:

```
import elfi
```

As is usual in Bayesian statistical inference, we need to define *prior*
distributions for the unknown parameters \(\theta_1, \theta_2\). In
ELFI the priors can be any of the continuous and discrete distributions
available in `scipy.stats`

(for custom priors, see
below). For simplicity, let’s start by assuming
that both parameters follow `Uniform(0, 2)`

.

```
# a node is defined by giving a distribution from scipy.stats together with any arguments (here 0 and 2)
t1 = elfi.Prior(scipy.stats.uniform, 0, 2)
# ELFI also supports giving the scipy.stats distributions as strings
t2 = elfi.Prior('uniform', 0, 2)
```

Next, we define the *simulator* node with the `MA2`

function above,
and give the priors to it as arguments. This means that the parameters
for the simulations will be drawn from the priors. Because we have the
observed data available for this node, we provide it here as well:

```
Y = elfi.Simulator(MA2, t1, t2, observed=y_obs)
```

But how does one compare the simulated sequences with the observed sequence? Looking at the plot of just a few observed sequences above, a direct pointwise comparison would probably not work very well: the three sequences look quite different although they were generated with the same parameter values. Indeed, the comparison of simulated sequences is often the most difficult (and ad hoc) part of ABC. Typically one chooses one or more summary statistics and then calculates the discrepancy between those.

Here, we will apply the intuition arising from the definition of the MA(2) process, and use the autocovariances with lags 1 and 2 as the summary statistics:

```
def autocov(x, lag=1):
C = np.mean(x[:,lag:] * x[:,:-lag], axis=1)
return C
```

As is familiar by now, a `Summary`

node is defined by giving the
autocovariance function and the simulated data (which includes the
observed as well):

```
S1 = elfi.Summary(autocov, Y)
S2 = elfi.Summary(autocov, Y, 2) # the optional keyword lag is given the value 2
```

Here, we choose the discrepancy as the common Euclidean L2-distance.
ELFI can use many common distances directly from
`scipy.spatial.distance`

like this:

```
# Finish the model with the final node that calculates the squared distance (S1_sim-S1_obs)**2 + (S2_sim-S2_obs)**2
d = elfi.Distance('euclidean', S1, S2)
```

One may wish to use a distance function that is unavailable in
`scipy.spatial.distance`

. ELFI supports defining a custom
distance/discrepancy functions as well (see the documentation for
`elfi.Distance`

and `elfi.Discrepancy`

).

Now that the inference model is defined, ELFI can visualize the model as a DAG.

```
elfi.draw(d) # just give it a node in the model, or the model itself (d.model)
```

Note

You will need the Graphviz software as well as the graphviz Python package (https://pypi.python.org/pypi/graphviz) for drawing this. The software is already installed in many unix-like OS.

## Modifying the model¶

Although the above definition is perfectly valid, let’s use the same priors as in *Marin et al. (2012)* that guarantee that the problem will be identifiable (loosely speaking, the likelihood willl have just one mode). Marin et al. used priors for which \(-2<\theta_1<2\) with \(\theta_1+\theta_2>-1\) and \(\theta_1-\theta_2<1\) i.e. the parameters are sampled from a triangle (see below).

### Custom priors¶

In ELFI, custom distributions can be defined similar to distributions in
`scipy.stats`

(i.e. they need to have at least the `rvs`

method
implemented for the simplest algorithms). To be safe they can inherit
`elfi.Distribution`

which defines the methods needed. In this case we
only need these for sampling, so implementing a static `rvs`

method
suffices. As was in the context of simulators, it is important to accept
the keyword argument `random_state`

, which is needed for ELFI’s
internal book-keeping of pseudo-random number generation. Also the
`size`

keyword is needed (which in the simple cases is the same as the
`batch_size`

in the simulator definition).

```
# define prior for t1 as in Marin et al., 2012 with t1 in range [-b, b]
class CustomPrior_t1(elfi.Distribution):
def rvs(b, size=1, random_state=None):
u = scipy.stats.uniform.rvs(loc=0, scale=1, size=size, random_state=random_state)
t1 = np.where(u<0.5, np.sqrt(2.*u)*b-b, -np.sqrt(2.*(1.-u))*b+b)
return t1
# define prior for t2 conditionally on t1 as in Marin et al., 2012, in range [-a, a]
class CustomPrior_t2(elfi.Distribution):
def rvs(t1, a, size=1, random_state=None):
locs = np.maximum(-a-t1, t1-a)
scales = a - locs
t2 = scipy.stats.uniform.rvs(loc=locs, scale=scales, size=size, random_state=random_state)
return t2
```

These indeed sample from a triangle:

```
t1_1000 = CustomPrior_t1.rvs(2, 1000)
t2_1000 = CustomPrior_t2.rvs(t1_1000, 1, 1000)
plt.scatter(t1_1000, t2_1000, s=4, edgecolor='none');
# plt.plot([0, 2, -2, 0], [-1, 1, 1, -1], 'b') # outlines of the triangle
```

Let’s change the earlier priors to the new ones in the inference model:

```
t1.become(elfi.Prior(CustomPrior_t1, 2))
t2.become(elfi.Prior(CustomPrior_t2, t1, 1))
elfi.draw(d)
```

Note that `t2`

now depends on `t1`

. Yes, ELFI supports hierarchy.

## Inference with rejection sampling¶

The simplest ABC algorithm samples parameters from their prior distributions, runs the simulator with these and compares them to the observations. The samples are either accepted or rejected depending on how large the distance is. The accepted samples represent samples from the approximate posterior distribution.

In ELFI, ABC methods are initialized either with a node giving the
distance, or with the `ElfiModel`

object and the name of the distance
node. Depending on the inference method, additional arguments may be
accepted or required.

A common optional keyword argument, accepted by all inference methods,
`batch_size`

defines how many simulations are performed in each
passing through the graph.

Another optional keyword is the seed. This ensures that the outcome will be always the same for the same data and model. If you leave it out, a random seed will be taken.

```
rej = elfi.Rejection(d, batch_size=10000, seed=seed)
```

Note

In Python, doing many calculations with a single function call can potentially save a lot of CPU time, depending on the operation. For example, here we draw 10000 samples from t1, pass them as input to t2, draw 10000 samples from t2, and then use these both to run 10000 simulations and so forth. All this is done in one passing through the graph and hence the overall number of function calls is reduced 10000-fold. However, this does not mean that batches should be as big as possible, since you may run out of memory, the fraction of time spent in function call overhead becomes insignificant, and many algorithms operate in multiples of batch_size. Furthermore, the batch_size is a crucial element for efficient parallelization (see the notebook on parallelization).

After the ABC method has been initialized, samples can be drawn from it.
By default, rejection sampling in ELFI works in `quantile`

mode i.e. a
certain quantile of the samples with smallest discrepancies is accepted.
The `sample`

method requires the number of output samples as a
parameter. Note that the simulator is then run `(N/quantile)`

times.
(Alternatively, the same behavior can be achieved by saying
`n_sim=1000000`

.)

The IPython magic command `%time`

is used here to give you an idea of
runtime on a typical personal computer. We will turn interactive
visualization on so that if you run this on a notebook you will see the
posterior forming from a prior distribution. In this case most of the
time is spent in drawing.

```
N = 1000
vis = dict(xlim=[-2,2], ylim=[-1,1])
# You can give the sample method a `vis` keyword to see an animation how the prior transforms towards the
# posterior with a decreasing threshold.
%time result = rej.sample(N, quantile=0.01, vis=vis)
```

```
CPU times: user 2.59 s, sys: 100 ms, total: 2.69 s
Wall time: 2.68 s
```

The `sample`

method returns a `Sample`

object, which contains
several attributes and methods. Most notably the attribute `samples`

contains an `OrderedDict`

(i.e. an ordered Python dictionary) of the
posterior numpy arrays for all the model parameters (`elfi.Prior`

s
in the model). For rejection sampling, other attributes include e.g. the
`threshold`

, which is the threshold value resulting in the requested
quantile.

```
result.samples['t1'].mean()
```

```
0.56
```

The `Sample`

object includes a convenient `summary`

method:

```
result.summary()
```

```
Method: Rejection
Number of samples: 1000
Number of simulations: 100000
Threshold: 0.117
Sample means: t1: 0.556, t2: 0.219
```

Rejection sampling can also be performed with using a threshold or total number of simulations. Let’s define here threshold. This means that all draws from the prior for which the generated distance is below the threshold will be accepted as samples. Note that the simulator will run as long as it takes to generate the requested number of samples.

```
%time result2 = rej.sample(N, threshold=0.2)
print(result2) # the Sample object's __str__ contains the output from summary()
```

```
CPU times: user 248 ms, sys: 12 ms, total: 260 ms
Wall time: 255 ms
Method: Rejection
Number of samples: 1000
Number of simulations: 40000
Threshold: 0.185
Sample means: t1: 0.555, t2: 0.223
```

## Iterative advancing¶

Often it may not be practical to wait to the end before investigating
the results. There may be time constraints or one may wish to check the
results at certain intervals. For this, ELFI provides an iterative
approach to advance the inference. First one sets the objective of the
inference and then calls the `iterate`

method.

Below is an example how to run the inference until the objective has been reached or a maximum of one second of time has been used.

```
# Request for 1M simulations.
rej.set_objective(1000, n_sim=1000000)
# We only have 1 sec of time and we are unsure if we will be finished by that time.
# So lets simulate as many as we can.
time0 = time.time()
time1 = time0 + 1
while not rej.finished and time.time() < time1:
rej.iterate()
# One could investigate the rej.state or rej.extract_result() here
# to make more complicated stopping criterions
# Extract and print the result as it stands. It will show us how many simulations were generated.
print(rej.extract_result())
```

```
Method: Rejection
Number of samples: 1000
Number of simulations: 150000
Threshold: 0.0968
Sample means: t1: 0.561, t2: 0.217
```

```
# We will see that it was not finished in 1 sec
rej.finished
```

```
False
```

We could continue from this stage just by continuing to call the
`iterate`

method. The `extract_result`

will always give a proper
result even if the objective was not reached.

Next we will look into how to store all the data that was generated so far. This allows us to e.g. save the data to disk and continue the next day, or modify the model and reuse some of the earlier data if applicable.

## Storing simulated data¶

As the samples are already in numpy arrays, you can just say e.g.
`np.save('t1_data.npy', result.samples['t1'])`

to save them. However,
ELFI provides some additional functionality. You may define a *pool* for
storing all outputs of any node in the model (not just the accepted
samples). Let’s save all outputs for `t1`

, `t2`

, `S1`

and `S2`

in our model:

```
pool = elfi.OutputPool(['t1', 't2', 'S1', 'S2'])
rej = elfi.Rejection(d, pool=pool)
%time result3 = rej.sample(N, n_sim=1000000)
result3
```

```
CPU times: user 9.95 s, sys: 0 ns, total: 9.95 s
Wall time: 9.95 s
```

```
Method: Rejection
Number of samples: 1000
Number of simulations: 1000000
Threshold: 0.0362
Sample means: t1: 0.554, t2: 0.233
```

The benefit of the pool is that you may reuse simulations without having to resimulate them. Above we saved the summaries to the pool, so we can change the distance node of the model without having to resimulate anything. Let’s do that.

```
# Replace the current distance with a cityblock (manhattan) distance and recreate the inference
d.become(elfi.Distance('cityblock', S1, S2, p=1))
rej = elfi.Rejection(d, pool=pool)
%time result4 = rej.sample(N, n_sim=1000000)
result4
```

```
CPU times: user 700 ms, sys: 0 ns, total: 700 ms
Wall time: 699 ms
```

```
Method: Rejection
Number of samples: 1000
Number of simulations: 1000000
Threshold: 0.0453
Sample means: t1: 0.555, t2: 0.235
```

Note the significant saving in time, even though the total number of considered simulations stayed the same.

We can also continue the inference by increasing the total number of simulations and only have to simulate the new ones:

```
%time result5 = rej.sample(N, n_sim=1200000)
result5
```

```
CPU times: user 2.07 s, sys: 16 ms, total: 2.09 s
Wall time: 2.09 s
```

```
Method: Rejection
Number of samples: 1000
Number of simulations: 1200000
Threshold: 0.0421
Sample means: t1: 0.554, t2: 0.239
```

Above the results were saved into a python dictionary. If you store a lot of data to dictionaries, you will eventually run out of memory. ELFI provides an alternative pool that, by default, saves the outputs to standard numpy .npy files:

```
arraypool = elfi.ArrayPool(['t1', 't2', 'Y', 'd'])
rej = elfi.Rejection(d, pool=arraypool)
%time result5 = rej.sample(100, threshold=0.3)
```

```
CPU times: user 40 ms, sys: 0 ns, total: 40 ms
Wall time: 37.1 ms
```

This stores the simulated data in binary `npy`

format under
`arraypool.path`

, and can be loaded with `np.load`

.

```
# Let's flush the outputs to disk (alternatively you can save or close the pool) so that we can read the .npy files.
arraypool.flush()
import os
print('Files in', arraypool.path, 'are', os.listdir(arraypool.path))
```

```
Files in pools/arraypool_3615052699 are ['Y.npy', 't2.npy', 'd.npy', 't1.npy']
```

Now lets load all the parameters `t1`

that were generated with numpy:

```
np.load(arraypool.path + '/t1.npy')
```

```
array([-0.82, -0.03, 0.27, ..., 1.03, 0.44, -0.56])
```

We can also close (or save) the whole pool if we wish to continue later:

```
arraypool.close()
name = arraypool.name
print(name)
```

```
arraypool_3615052699
```

And open it up later to continue where we were left. We can open it using its name:

```
arraypool = elfi.ArrayPool.open(name)
print('This pool has', len(arraypool), 'batches')
# This would give the contents of the first batch
# arraypool[0]
```

```
This pool has 3 batches
```

You can delete the files with:

```
arraypool.delete()
# verify the deletion
try:
os.listdir(arraypool.path)
except FileNotFoundError:
print("The directry is removed")
```

```
The directry is removed
```

## Visualizing the results¶

Instances of `Sample`

contain methods for some basic plotting (these
are convenience methods to plotting functions defined under
`elfi.visualization`

).

For example one can plot the marginal distributions:

```
result.plot_marginals();
```

Often “pairwise relationships” are more informative:

```
result.plot_pairs();
```

Note that if working in a non-interactive environment, you can use e.g.
`plt.savefig('pairs.png')`

after an ELFI plotting command to save the
current figure to disk.

## Sequential Monte Carlo ABC¶

Rejection sampling is quite inefficient, as it does not learn from its
history. The sequential Monte Carlo (SMC) ABC algorithm does just that
by applying importance sampling: samples are *weighed* according to the
resulting discrepancies and the next *population* of samples is drawn
near to the previous using the weights as probabilities.

For evaluating the weights, SMC ABC needs to be able to compute the
probability density of the generated parameters. In our MA2 example we
used custom priors, so we have to specify a `pdf`

function by
ourselves. If we used standard priors, this step would not be needed.
Let’s modify the prior distribution classes:

```
# define prior for t1 as in Marin et al., 2012 with t1 in range [-b, b]
class CustomPrior_t1(elfi.Distribution):
def rvs(b, size=1, random_state=None):
u = scipy.stats.uniform.rvs(loc=0, scale=1, size=size, random_state=random_state)
t1 = np.where(u<0.5, np.sqrt(2.*u)*b-b, -np.sqrt(2.*(1.-u))*b+b)
return t1
def pdf(x, b):
p = 1./b - np.abs(x) / (b*b)
p = np.where(p < 0., 0., p) # disallow values outside of [-b, b] (affects weights only)
return p
# define prior for t2 conditionally on t1 as in Marin et al., 2012, in range [-a, a]
class CustomPrior_t2(elfi.Distribution):
def rvs(t1, a, size=1, random_state=None):
locs = np.maximum(-a-t1, t1-a)
scales = a - locs
t2 = scipy.stats.uniform.rvs(loc=locs, scale=scales, size=size, random_state=random_state)
return t2
def pdf(x, t1, a):
locs = np.maximum(-a-t1, t1-a)
scales = a - locs
p = scipy.stats.uniform.pdf(x, loc=locs, scale=scales)
p = np.where(scales>0., p, 0.) # disallow values outside of [-a, a] (affects weights only)
return p
# Redefine the priors
t1.become(elfi.Prior(CustomPrior_t1, 2, model=t1.model))
t2.become(elfi.Prior(CustomPrior_t2, t1, 1))
```

### Run SMC ABC¶

In ELFI, one can setup a SMC ABC sampler just like the Rejection sampler:

```
smc = elfi.SMC(d, batch_size=10000, seed=seed)
```

For sampling, one has to define the number of output samples, the number
of populations and a *schedule* i.e. a list of quantiles to use for each
population. In essence, a population is just refined rejection sampling.

```
N = 1000
schedule = [0.7, 0.2, 0.05]
%time result_smc = smc.sample(N, schedule)
```

```
CPU times: user 6.46 s, sys: 132 ms, total: 6.59 s
Wall time: 1.83 s
```

We can have summaries and plots of the results just like above:

```
result_smc.summary(all=True)
```

```
Method: SMC
Number of samples: 1000
Number of simulations: 190000
Threshold: 0.0491
Sample means: t1: 0.556, t2: 0.22
Population 0:
Method: Rejection within SMC-ABC
Number of samples: 1000
Number of simulations: 10000
Threshold: 0.488
Sample means: t1: 0.547, t2: 0.232
Population 1:
Method: Rejection within SMC-ABC
Number of samples: 1000
Number of simulations: 20000
Threshold: 0.185
Sample means: t1: 0.556, t2: 0.236
Population 2:
Method: Rejection within SMC-ABC
Number of samples: 1000
Number of simulations: 160000
Threshold: 0.0491
Sample means: t1: 0.556, t2: 0.22
```

Or just the means:

```
result_smc.sample_means_summary(all=True)
```

```
Sample means for population 0: t1: 0.547, t2: 0.232
Sample means for population 1: t1: 0.556, t2: 0.236
Sample means for population 2: t1: 0.556, t2: 0.22
```

```
result_smc.plot_marginals(all=True, bins=25, figsize=(8, 2), fontsize=12)
```

Obviously one still has direct access to the samples as well, which allows custom plotting:

```
n_populations = len(schedule)
fig, ax = plt.subplots(ncols=n_populations, sharex=True, sharey=True, figsize=(16,6))
for i, pop in enumerate(result_smc.populations):
s = pop.samples
ax[i].scatter(s['t1'], s['t2'], s=5, edgecolor='none');
ax[i].set_title("Population {}".format(i));
ax[i].plot([0, 2, -2, 0], [-1, 1, 1, -1], 'b')
ax[i].set_xlabel('t1');
ax[0].set_ylabel('t2');
ax[0].set_xlim([-2, 2])
ax[0].set_ylim([-1, 1]);
```

It can be seen that the populations iteratively concentrate more and more around the true parameter values.

Note that for the later populations some of the samples lie outside allowed region. This is due to the SMC algorithm sampling near previous samples. However, the outliers carry zero weight and have no effect on the estimates.

That’s it! See the other documentation for more advanced topics on e.g. BOLFI, external simulators and parallelization.