ABC rejection comparing two time series in deterministic Lotka-Volterra model about elfi HOT 5 CLOSED

Hi, the problem is that ELFI expects the outputs of these operations to be numpy arrays of length batch_size (even with batch_size=1). Currently your distance function SSE returns a numpy scalar. So you can either vectorize it with elfi.tools.vectorize, or as a quick solution to write e.g.

def SSE(x, y):
   return np.atleast_1d(np.sum((x-y)**2))

However your simulator should also wrap the results in numpy arrays with length batch_size (in this case it would be a 2D numpy array), so that it would work with other batch sizes besides 1. You can either do that yourself or use the elfi.tools.vectorize. This in turn then requires that your distance handles them correctly (assuming your observed data is also 2D array):

def SSE(x, y):
    return np.sum((x-y)**2, axis=1)

We should think a way to detect these cases and report a better error message since these seem to be quite common.

from elfi.

Thanks for the quick response and suggested adjustment.

I made the adjustment to the distance function and continued on with the example. However, now the rejection sampler hangs and if I kill the process (ctrl-C) and plot the parameter space that's being explored then it looks like it's searching in the correct area but it seems to take a while to reach the result (roughly a,b = 1). I was curious about ways to speed this up?

# (as example above ...)

def SSE(x, y):
   return np.atleast_1d(np.sum((x-y)**2))

rej = elfi.Rejection(d_lv, batch_size = 1, seed = seed)
result = rej.sample(1000, threshold = 30.0)

# kill the process after about a minute (ctrl-C)

# Plot the parameter space that's been explored
fig, ax = plt.subplots()
rej.plot_state(ax = ax)
ax.set_xlim([-11, 11])
ax.set_ylim([-11, 11])
plt.savefig("parameter_space.png")
plt.close()

Additionally, when outputting the result with rej.extract_result() the threshold parameter seems to have been ignored (it was set to 30.0).

rej.extract_result()

Method: Rejection
Number of samples: 1000
Number of simulations: 21063
Threshold: 484
Sample means: a: -0.518, b: 2

from elfi.

There are several things you can try to speed things up, for example:

1. Is requesting threshold=30 reasonable? ELFI will keep sampling until it finds 1000 samples below that threshold, which for a very small threshold may never happen. To be safe, you could initially request e.g. quantile=0.01 instead of the threshold. ELFI will then draw 1000 * 1 / 0.01 samples, return the best 1000 and report the resulting threshold value. That said, it seems ELFI ran only 21063 simulations in your case, so even this would apparently take several minutes.

2. Enable parallelization. See the documentation.

3. Use a more advanced inference method, e.g. BOLFI.

4. This is probably the most effective way, but also the most difficult (and not always even possible): try to vectorize your simulator so that you could use and benefit from a larger batch_size.

Since you killed the process, what you get using extract_result may not be what you want. I suspect there were not enough samples below the given threshold.

from elfi.

Thanks for those ideas.

I used a threshold of 30.0 because that is what is quoted as the first threshold in ABC-SMC in the paper by Toni et al., (2009) (assuming the first threshold of ABC-SMC would be a suitable threshold for ABC-rejection). They report ~26,000 data generation steps to find 1000 accepted samples (with ABC-SMC).

from elfi.

Also ABC-SMC speeds up the inference compared to Rejection ABC. You might be interested in the iterative advancing of the algorithm:

http://elfi.readthedocs.io/en/latest/usage/tutorial.html#iterative-advancing

Basically you set the objective and then you can advance the inference step by step and investigate the result at any point. Internally the sample method runs the iterative advancing loop until the objective is reached.

from elfi.