# Adaptive distance

ABC provides means to sample an approximate posterior distribution over unknown parameters based on comparison between observed and simulated data. This comparison is often based on distance between features that summarise the data and are informative about the parameter values.

Here we assume that the summaries calculated based on observed and simulated data are compared based on weighted distance with weight \(w_i=1/\sigma_i\) calculated based on their standard deviation \(\sigma_i\). This ensures that the selected summaries to have an equal contribution in the distance between observed and simulated data.

This notebook studies adaptive distance SMC-ABC where \(\sigma_i\) and \(w_i\) are recalculated between SMC iterations as proposed in [1].

```
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
%matplotlib inline
```

```
import elfi
```

## Example 1:

Assume we have an unknown parameter with prior distribution \(\theta\sim U(0,50)\) and two simulator outputs \(S_1\sim N(\theta, 1)\) and \(S_2\sim N(\theta, 100)\) whose observed values are 20.

```
def simulator(mu, batch_size=1, random_state=None):
batches_mu = np.asarray(mu).reshape((-1,1))
obs_1 = ss.norm.rvs(loc=batches_mu, scale=1, random_state=random_state).reshape((-1,1))
obs_2 = ss.norm.rvs(loc=batches_mu, scale=100, random_state=random_state).reshape((-1,1))
return np.hstack((obs_1, obs_2))
```

```
observed_data = np.array([20,20])[None,:]
```

Here the simulator outputs are both informative about the unknown model parameter, but \(S_2\) has more observation noise than \(S_1\). We do not calculate separate summaries in this example, but compare observed and simulated data based on these two variables.

Euclidean distance between observed and simulated outputs or summaries
can be used to find parameter values that could produce the observed
data. Here we describe dependencies between the unknown parameter value
and observed distances as an ELFI model `m`

and sample the approximate
posterior distribution with the rejection
sampler.

```
m = elfi.new_model()
theta = elfi.Prior(ss.uniform, 0, 50, model=m)
sim = elfi.Simulator(simulator, theta, observed=observed_data)
d = elfi.Distance('euclidean', sim)
```

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

Let us sample 100 parameters with `quantile=0.01`

. This means that we
sample 10000 candidate parameters from the prior distribution and take
the 100 parameters that produce simulated data closest to the observed
data.

```
sample = rej.sample(100, quantile=0.01)
```

```
Progress [==================================================] 100.0% Complete
```

```
sample
```

```
Method: Rejection
Number of samples: 100
Number of simulations: 10000
Threshold: 6.66
Sample means: theta: 19.6
```

```
plt.hist(sample.samples_array,range=(0,50),bins=20)
plt.xlabel('theta');
```

The approximate posterior sample is concentrated around \(\theta=20\) as expected in this example. However the sample distribution is much wider than we would observe in case the sample was selected based on \(S_1\) alone.

Now let us test adaptive distance in the same example.

First we switch the distance node `d`

to an adaptive distance node and
initialise adaptive distance SMC-ABC. Initialisation is identical to the
rejection sampler, and here we use the same batch size and seed as
earlier, so that the methods are presented with the exact same candidate
parameters.

```
d.become(elfi.AdaptiveDistance(sim))
```

```
ada_smc = elfi.AdaptiveDistanceSMC(d, batch_size=10000, seed=123)
```

Since this is an iterative method, we must decide both sample size
(`n_samples`

) and how many populations are sampled (`rounds`

). In
addition we can decide the \(\alpha\) quantile (`quantile`

) used
in estimation.

Each population with `n_samples`

parameter values is sampled as
follows: 1. `n_samples/quantile`

parameters are sampled from the
current proposal distribution with acceptance threshold determined based
on the previous population and 2. the distance measure is updated based
on the observed sample and `n_samples`

with the smallest distance are
selected as the new population. The first population is sampled from the
prior distribution and all samples are accepted in step 1.

Here we sample one population with `quantile=0.01`

. This means that
the total simulation count will be the same as with the rejection
sampler, but now the distance function is updated based on the 10000
simulated observations, and the 100 parameters included in the posterior
sample are selected based on the new distance measure.

```
sample_ada = ada_smc.sample(100, 1, quantile=0.01)
```

```
ABC-SMC Round 1 / 1
Progress [==================================================] 100.0% Complete
```

```
sample_ada
```

```
Method: AdaptiveDistanceSMC
Number of samples: 100
Number of simulations: 10000
Threshold: 0.462
Sample means: theta: 19.8
```

```
plt.hist(sample_ada.samples_array,range=(0,50),bins=20)
plt.xlabel('theta');
```

We see that the posterior distribution over unknown parameter values is narrower than in the previous example. This is because the simulator outputs are now normalised based on their estimated standard deviation.

We can see \(w_1\) and \(w_2\):

```
sample_ada.adaptive_distance_w
```

```
[array([0.06940134, 0.0097677 ])]
```

## Example 2:

This is the normal distribution example presented in [1].

Here we have an unknown parameter with prior distribution \(\theta\sim N(0,100)\) and two simulator outputs \(S_1\sim N(\theta, 0.1)\) and \(S_2\sim N(1, 1)\) whose observed values are 0.

```
def simulator(mu, batch_size=1, random_state=None):
batches_mu = np.asarray(mu).reshape((-1,1))
obs_1 = ss.norm.rvs(loc=batches_mu, scale=0.1, random_state=random_state).reshape((-1,1))
obs_2 = ss.norm.rvs(loc=1, scale=1, size=batch_size, random_state=random_state).reshape((-1,1))
return np.hstack((obs_1, obs_2))
```

```
observed_data = np.array([0,0])[None,:]
```

\(S_1\) is now informative and \(S_2\) uninformative about the unknown parameter value, and we note that between the two output variables, \(S_1\) has larger variance under the prior predictive distribution. This means that normalisation estimated based on output data observed in the initial round or based on a separate sample would not work well in this example.

Let us define a new model and initialise adaptive distance SMC-ABC.

```
m = elfi.new_model()
theta = elfi.Prior(ss.norm, 0, 100, model=m)
sim = elfi.Simulator(simulator, theta, observed=observed_data)
d = elfi.AdaptiveDistance(sim)
```

```
ada_smc = elfi.AdaptiveDistanceSMC(d, batch_size=2000, seed=123)
```

Next we sample 1000 parameter values in 5 rounds with the default
`quantile=0.5`

which is recommended in sequential estimation
[1]:

```
sample_ada = ada_smc.sample(1000, 5)
```

```
ABC-SMC Round 1 / 5
Progress [==================================================] 100.0% Complete
ABC-SMC Round 2 / 5
Progress [==================================================] 100.0% Complete
ABC-SMC Round 3 / 5
Progress [==================================================] 100.0% Complete
ABC-SMC Round 4 / 5
Progress [==================================================] 100.0% Complete
ABC-SMC Round 5 / 5
Progress [==================================================] 100.0% Complete
```

```
sample_ada
```

```
Method: AdaptiveDistanceSMC
Number of samples: 1000
Number of simulations: 32000
Threshold: 0.925
Sample means: theta: -0.195
```

```
plt.hist(sample_ada.samples_array, range=(-25,25), bins=20)
plt.xlabel(theta);
```

The sample distribution is concentrated around \(\theta=0\) but wider than could be expected. However we can continue the iterative estimation process. Here we sample two more populations:

```
sample_ada = ada_smc.sample(1000, 2)
```

```
ABC-SMC Round 6 / 7
Progress [==================================================] 100.0% Complete
ABC-SMC Round 7 / 7
Progress [==================================================] 100.0% Complete
```

```
sample_ada
```

```
Method: AdaptiveDistanceSMC
Number of samples: 1000
Number of simulations: 48000
Threshold: 0.868
Sample means: theta: 0.0183
```

```
plt.hist(sample_ada.samples_array, range=(-25,25), bins=20)
plt.xlabel('theta');
```

We observe that the sample mean is now closer to zero and the sample distribution is narrower.

Let us examine \(w_1\) and \(w_2\):

```
sample_ada.adaptive_distance_w
```

```
[array([0.01023228, 1.00584519]),
array([0.00921258, 0.99287166]),
array([0.01201937, 0.99365522]),
array([0.02217631, 0.98925365]),
array([0.04355987, 1.00076738]),
array([0.07863284, 0.9971017 ]),
array([0.13892778, 1.00929049])]
```

We can see that \(w_2\) (second column) is constant across iterations whereas \(w_1\) increases as the method learns more about possible parameter values and the proposal distribution becomes more concentrated around \(\theta=0\).

## Notes

The adaptive distance SMC-ABC method demonstrated in this notebook normalises simulator outputs or summaries calculated based on simulator output based on their estimated standard deviation under the proposal distribution in each iteration. This ensures that all outputs or summaries have an equal contribution to the distance between simulated and observed data in all iterations.

It is important to note that the method does not evaluate whether outputs or summaries are needed or informative. In both examples studied in this notebook, results would improve if inference was carried out based on \(S_1\) alone. Hence one should choose the summaries used in adaptive distance SMC-ABC with the usual care. ELFI tools that aid in the selection process are discussed in the diagnostics notebook available here.

## Reference

[1] Prangle D (2017). Adapting the ABC Distance Function. Bayesian Analysis 12(1): 289-309, 2017. https://projecteuclid.org/euclid.ba/1460641065