[ANN-RFC] KissABC.jl- Approximate Bayesian Computation

Package Announcements
,
42

yes i’m planning to improve docs, “AIS” stands for affine invariant sampler, since it implements mcmc moves which are affine invariant, like the DE move (which was the only move of ABCDE).

as a start to imitate old behavior you could try to use it like:

sample(model, AIS(N), N, burnin=iterations)
# N = nparticles parameter of ABCDE
# iterations = generations parameter of ABCDE

from there you can start exploring the new options, let me know what you find, 2.0 is a complete redesign and very breaking release, if it really performs worse on your problems you can momentarily use 1.3 via:

]add KissABC@1.3
]pin KissABC

and i will work on a fix and a new release :slight_smile:
btw is there a chance you can send me a toy version of the kind of problems you feed to kissabc?

1 Like
43

Thanks, I’ll give that a try and let you know how it goes!

It’s mostly parameter estimation for differential equation models from biology.

1 Like
44

I would love a toy example on that, it would make for a nice tutorial and benchmark, now i’m almost completely satisfied with the API of KissABC, so the next release will focus on perf improvements and new algorithms (i removed ABCSMCPR, but surely a new SMC algorithm will come into play, possibly exploiting the efficient proposal mechanism of AIS)

45

Ah ok, sure when I have time I will.

Could you explain the output of sample quickly? What is the third index of that resulting array? The previous output was just the nparticles particles that were below the threshold.

Secondly, when fitting ODE models, some parameter combinations tend to cause the solver to stop because they’re too stiff. Usually I just return a cost of Inf, is there a better way to do it?

edit: it seems like returning threshold+1 works, but is there a less hacky way?

Thanks!

1 Like
46

the output of sample is a type AISChain which behaves as an array with 3 axis, samples, parameters, and chains thus the code

results = sample(...)
print(results[1,2,3])

will print the second parameter of the first sample for the third chain, if you call sample without MCMCThreads enabled, you will only have one chain.
you can also slice results (https://juliaapproxinference.github.io/KissABC.jl/stable/reference/#KissABC.AISChain) like you would do for any other array of size nsamples × nparameters × nchains.

the nice thing is that the type AISChain also provides a summary of the inferred parameters, and it can also be fed directly to the more feature complete Chains type defined in MCMCChains.jl by simply calling Chains(results)

if instead you want access to raw samples like you did in past versions you can simply write results.samples, samples is the field of AISChain containing the raw data produced by the method sample

edit. for the stiffness problem, if you can find a way to define a cost which directly targets stiffness you can include a penalization for stiff solutions, in any case if you have the time open an issue on KissABC for supporting Inf as a cost value

47

Thanks for the reply,

I’ll get you a test example for your fix on the issue I posted when I have the time!

I still have some questions about the AISChain sample output. I think I might be misunderstanding some things. With the past version, calling ABCDE on 8 threads still gave me nparticles samples. However, here it appears that I now get 8*nparticles samples, and some of them are not within threshold. Am I misunderstanding the output of this? I am pretty new to bayesian approximation in general so I think that’s likely.

Also, do you have a citation for the algorithm used by AIS()? I understand that it’s a generalization of the ABCDE algorithm?

Thanks again!

48

@pjentsch0 are you using ApproxPosterior or ApproxKernelizedPosterior?

In the second case, It Is normal (no pun intended) since thresholds are soft and errors are assumed to be gaussian distributed

In the first case it can be fixed by increasing burnin

You i do have the citations,

https://msp.org/camcos/2010/5-1/p04.xhtml

Regarding the difference in parallelism, It Is a byproduct of the use of AbstractMCMC, which on one hand gives both threaded and distributed parallelism, on the other hand It multiplies the number of particles. The good thing Is that it’s approach Is guaranteed to give reproducible results if a seed on the rng is set.

1 Like
49

I am working with a problem where drawing a sample takes about 15 seconds. The model has 8 parameters. I am running on a multiprocessor machine with 16 cores. Currently, I’m experimenting with KissABC using something like

approx_density = ApproxPosterior(prior, D, 0.01)
res = sample(
    approx_density,
    AIS(100),
    MCMCThreads(),
    10000,
    10,
    burnin = 100,
    ntransitions = 10,
    progress = true,
)

So, I’m using 10 threads. This does give reasonable results, but it takes quite a while to run, and tuning experimentally is consequently somewhat difficult. I would appreciate any recommendations for algorithm choice, and tuning parameters for a problem like this. Thanks!

1 Like
50

Is KissABC similar to INLA?

51

Announcement: KissABC 3.0 Released

Breaking Changes:

  • the return type of inference algorithms is now a vector of MonteCarloMeasurements particles (this makes the results readily accessible for further processing)

New algorithms:

the new algorithm can reach much lower tolerances in much less function evaluations (keep in mind that the default parameters are a bit conservative).

Changes:

  • The package exports all symbols from MonteCarloMeasurements and Distributions, for ease of use

Hope you guys enjoy it! :slight_smile:

4 Likes
52

Doesn’t this mean that every time MonteCarloMeasurements or Distributions have a new breaking release you need to bump the major version of your own package? Distributions alone has had 22 breaking releases so far

3 Likes
53

Priors for algorithms are written using Distributions
and results are Particles object
so to use KissABC without the rexported symbols, users should always be using Distributions, MonteCarloMeasurements.
Furthermore the user is free to import KissABC and not pollute the namespace.

I can see the downsides too, maybe It was a bad decision, what would you suggest?

If the sentiment Is shared, i will remove the reexport

2 Likes
54

Thanks for the update, I’ve been using your package for my work since 1.3.0, and unfortunately I am still on 1.3 as the ABCDE method seems to do much better on my problem. I’ve been trying to come up with a MWE example for you but am not able to replicate the issue outside of my code.

I have looked into some other ABC packages, but it doesn’t seem there are any other julia packages that implement algorithms to handle ABC parameter estimation with ~30 variables, and I’d rather not move to python.

I can show some results from fitting here:
with KissABCv1.3

fakedist(x,y)=x+y #workaround for combined cost/simulate function you mentioned
plan = ABCplan(prior_dist, simulate, 0.0, fakedist)
res,_ = ABCDE(plan,16.0,nparticles = nparticles,generations = 0.5e3,earlystop = true,verbose=true,parallel = true)

kissabc1_3
kissabc1_3955×554 33.8 KB

with KissABCv3.0

res = smc(prior_dist,simulate; nparticles = 500,verbose = true,parallel = true,M = 10).P

smc
smc990×554 33.8 KB

using

approx_density = ApproxKernelizedPosterior(prior_dist,simulate,0.01)
res = sample(approx_density,AIS(30),MCMCThreads(),300,12;progress = true,ntransitions = 1600)

approx_kernalized_posterior
approx_kernalized_posterior977×554 23.2 KB

I have tried many parameter combinations to get the resulting distributions as tight as v1.3 with ABCDE, do you have any suggestions?

edit2: added examples

1 Like
55

The quality of results from smc Is Indeed troubling, can you share (privately) the code that produces this results?

atm i’m hoping it can be solved via parameter tuning :smiley: , i suggest trying r_epstol=0, mcmc_tol=0, epstol=16.0, verbose=true, min_r_ess=1.0, M=1, and paste the output here

if the algorithm looks slow, alpha can be reduced, alpha=0.8 perhaps

1 Like
56

The quality of results from smc Is Indeed troubling, can you share (privately) the code that produces this results?

I will discuss with my supervisor for sure. Hopefully we can try solving it with parameter tuning first?

With these parameter the algorithm seems stuck at the same eps for about 700 iterations:

(iteration, ϵ, dest, target) = (1061, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1062, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1063, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1064, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1065, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1066, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1067, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1068, 176.00470826426275, 499.9999999999995, 400.0)
(iteration, ϵ, dest, target) = (1069, 176.00470826426275, 499.9999999999995, 400.0)

This occurs also with alpha = 0.8.

57

Hmm this seems to be a problem of degeneration of particles, are there infinite costs involved?

In such case i would suggest using more particles (try doubling maybe)

58

There are no infinite costs.

I tried 2x the particles which didn’t have much effect. 4x the particles does indeed improve things, but we run into the same issue:

(iteration, ϵ, dest, target) = (222, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (223, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (224, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (225, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (226, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (227, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (228, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (229, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (230, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (231, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (232, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (233, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (234, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (235, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (236, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (237, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (238, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (239, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (240, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (241, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (242, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (243, 83.2524411348279, 2000.0000000000043, 1600.0)
(iteration, ϵ, dest, target) = (244, 83.2524411348279, 2000.0000000000043, 1600.0)

the function call, for reference:

res = smc(prior_dist,simulate; nparticles = 2000,verbose = true,parallel = true, r_epstol=0, mcmc_tol=0, epstol=16.0, min_r_ess=1.0, M=1, alpha = 0.8).P

If the results from 1.3.0 are correct, I can keep using that version for the time being.

59

i would really love a piece of code to help improve things… ABCDE (with ‘earlystop=true’) results are biased… this is the reason i cannot simply use the past algorithm and be done with it.

1 Like
60

Oh, the earlystop=true was only for quickly testing the code, usually I have earlystop=false. I Pmd you with a link to the repo.

1 Like
61

hey! First, huge thanks for making this package!

I want to try it out for my problems as well and saw that the ABCDE method is on GitHub master, but not yet in the latest release (KissABC v3.0.1). Is it still considered “work-in-progress” or somewhat safe to use? Or is the smc method the only recommended one?

Happy to hear about any info on this.

EDIT: Found a potential bug when using the (unreleased) ABCDE method from master branch with discrete priors (also here). Glad for any help!