Adaptive Distances

In this example, we show how and when to use the adaptive distances feature of pyabc. “Adaptive distances” means that the distance function is not pre-defined (e.g. after pre-processing), but evolves over time during the ABC run, depending on the observed summary statistics. This can be useful if different summary statistics vary on different scales, but it is not immediately clear how to weight them. For this case, in adaptive distances weights are adjusted in each iteration so as to balance the impact of all summary statistics on the computed distance.

Currently, adaptively weighted p-norm distances (e.g. Euclidean) are implemented in pyABC, but it is easily possible to define arbitrary adaptive distances.

The notebook can be downloaded here.

For illustration, we consider a simple Gaussian model:

In [1]:
import scipy
import tempfile
import os
import matplotlib.pyplot as pyplot
import pyabc.visualization
import logging


# for debugging
df_logger = logging.getLogger('DistanceFunction')
df_logger.setLevel(logging.DEBUG)

# model definition
def model(p):
    return {'ss1': p['theta'] + 1 + 0.1*scipy.randn(),
            'ss2': 2 + 10*scipy.randn()}

# true model parameter
theta_true = 3

# observed summary statistics
observation = {'ss1': theta_true + 1, 'ss2': 2}

# prior distribution
prior = pyabc.Distribution(theta=pyabc.RV('uniform', 0, 10))

Summary statistic ss2 has a high variance compared to summary statistic ss1. In addition, ss1 is informative about the model parameters \(\theta\), ss2 not. We expect that the proposal distribution for \(\theta\) iteratively centers around the true value \(\theta=3\). Thus, the variability for the sampled ss1 decreases iteratively, while the variability of the sampled ss2 stays approximately constant. If both summary statistics are weighted similarly in the calculation of the distance between sample and observation, there is hence an undesirable high impact of ss2, so that convergence can be slowed down. In contrast, if we weight ss1 higher, we may hope that our estimation of \(\theta\) is improved.

These informal expectations being stated, let us continue with the implementation. First, we consider a non-adaptive Euclidean distance:

In [2]:
distance = pyabc.PNormDistance(p=2)

abc = pyabc.ABCSMC(model, prior, distance)
db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "as1.db")
abc.new(db_path, observation)
history1 = abc.run(minimum_epsilon=.1, max_nr_populations=8)
INFO:History:Start <ABCSMC(id=1, start_time=2018-05-12 15:51:28.717308, end_time=None)>
INFO:Epsilon:initial epsilon is 8.761244749689185
INFO:ABC:t:0 eps:8.761244749689185
INFO:ABC:t:1 eps:5.411468424277703
INFO:ABC:t:2 eps:3.8525472217618604
INFO:ABC:t:3 eps:2.5727705794649913
INFO:ABC:t:4 eps:1.7124796517274543
INFO:ABC:t:5 eps:1.258749244114082
INFO:ABC:t:6 eps:0.9035495416544831
INFO:ABC:t:7 eps:0.5837176846030551
INFO:History:Done <ABCSMC(id=1, start_time=2018-05-12 15:51:28.717308, end_time=2018-05-12 15:52:57.608529)>

Let us visualize the results for the non-adaptive distance:

In [3]:
fig, ax = pyplot.subplots()
for t in range(history1.max_t + 1):
    df, w = history1.get_distribution(m=0, t=t)
    pyabc.visualization.plot_kde_1d(df, w, xmin=0, xmax=10,
                                   x='theta', ax=ax,
                                   label="PDF t={}".format(t))
ax.axvline(theta_true, color='k', linestyle='dashed')
ax.legend()
pyplot.show()
../_images/examples_adaptive_distances_9_0.png

Second, we consider an adaptive Euclidean distance:

In [4]:
distance = pyabc.AdaptivePNormDistance(p=2, adaptive=True)

abc = pyabc.ABCSMC(model, prior, distance)
db_path = "sqlite:///" + os.path.join(tempfile.gettempdir(), "as2.db")
abc.new(db_path, observation)
history2 = abc.run(minimum_epsilon=.1, max_nr_populations=8)
INFO:History:Start <ABCSMC(id=1, start_time=2018-05-12 15:53:00.365568, end_time=None)>
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.5133804231283783, 'ss2': 0.4866195768716217}
INFO:Epsilon:initial epsilon is 6.247205118760612
INFO:ABC:t:0 eps:6.247205118760612
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.5493899219352016, 'ss2': 0.4506100780647983}
INFO:ABC:t:1 eps:4.060230861877399
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.6756806547880043, 'ss2': 0.3243193452119958}
INFO:ABC:t:2 eps:2.271919498802732
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.7696321015615883, 'ss2': 0.23036789843841163}
INFO:ABC:t:3 eps:1.4842497319988184
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.857883245151725, 'ss2': 0.14211675484827477}
INFO:ABC:t:4 eps:0.8710188136983389
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.8980440690962745, 'ss2': 0.10195593090372569}
INFO:ABC:t:5 eps:0.542371204910505
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.941163988630602, 'ss2': 0.05883601136939804}
INFO:ABC:t:6 eps:0.3062394212285337
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.9529909029512236, 'ss2': 0.047009097048776496}
INFO:ABC:t:7 eps:0.19146509889020086
DEBUG:DistanceFunction:update distance weights = {'ss1': 1.9660609934351712, 'ss2': 0.03393900656482895}
INFO:History:Done <ABCSMC(id=1, start_time=2018-05-12 15:53:00.365568, end_time=2018-05-12 15:53:33.861279)>

In the debug output of abc.run above, it can be seen how the weights evolve over time. Let us visualize the results for the adaptive distance:

In [5]:
fig, ax = pyplot.subplots()
for t in range(history2.max_t + 1):
    df, w = history2.get_distribution(m=0, t=t)
    pyabc.visualization.plot_kde_1d(df, w, xmin=0, xmax=10,
                                   x='theta', ax=ax,
                                   label="PDF t={}".format(t))
ax.axvline(theta_true, color='k', linestyle='dashed')
ax.legend()
pyplot.show()
../_images/examples_adaptive_distances_13_0.png

We observe differences compared to the non-adaptive setting. In particular, the densitities tend to be narrower around the true parameter \(\theta=3\).

In detail, the adaptive distance feature works as follows: In each iteration of the ABCSMC run, after having obtained the desired number of accepted particles (and once at the beginning using a sample from the prior), the method DistanceFunction.update() is called. It is given a set of summary statistics which can be used to e.g. compute weights for the distance measure in the next iteration. In order to avoid bias, via DistanceFunction.configure_sampler(), the distance function can tell the sampler to not only record accepted particles, but all that were generated during the sampling process. So, when you want to define your own adaptive distance function, you will typically only need to overwrite these two methods. For implementation details and an example of how this can look in practice, please inspect the code of AdaptivePNormDistance.