# \[ANN-RFC\] KissABC.jl- Approximate Bayesian Computation
**URL:** https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668
**Category:** Package Announcements
**Tags:** announcement, bayesian-inference
**Created:** [June 3, 2020, 1:25pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668 "2020-06-03T13:25:12Z")
**Posts on this page:** 20
**Page:** 3
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [July 6, 2020, 3:57am UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/42 "2020-07-06T03:57:02Z")
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:
```julia
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:
```julia
]add KissABC@1.3
]pin KissABC
```
and i will work on a fix and a new release 🙂
btw is there a chance you can send me a toy version of the kind of problems you feed to kissabc?
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [July 8, 2020, 6:25pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/43 "2020-07-08T18:25:33Z")
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.
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [July 8, 2020, 7:06pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/44 "2020-07-08T19:06:41Z")
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`)
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [July 8, 2020, 9:08pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/45 "2020-07-08T21:08:47Z")
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!
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [July 8, 2020, 11:05pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/46 "2020-07-08T23:05:45Z")
the output of sample is a type `AISChain` which behaves as an array with 3 axis, `samples`, `parameters`, and `chains` thus the code
```julia
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](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
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [July 9, 2020, 6:40pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/47 "2020-07-09T18:40:21Z")
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!
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [July 9, 2020, 11:01pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/48 "2020-07-09T23:01:48Z")
@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,
> **[A Markov Chain Monte Carlo version of the genetic algorithm Differential...](https://link.springer.com/article/10.1007/s11222-006-8769-1)**
>
> Differential Evolution (DE) is a simple genetic algorithm for numerical optimization in real parameter spaces. In a statistical context one would not just want the optimum but also its uncertainty. The uncertainty distribution can be obtained by a...
[https://msp.org/camcos/2010/5-1/p04.xhtml](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.
---
### Author:  [@mcreel](https://discourse.julialang.org/u/mcreel)
#### Post date: [July 20, 2020, 4:42pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/49 "2020-07-20T16:42:10Z")
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
```julia
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!
---
### Author:  [@Juan](https://discourse.julialang.org/u/Juan)
#### Post date: [July 20, 2020, 5:18pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/50 "2020-07-20T17:18:19Z")
Is KissABC similar to INLA?
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [September 7, 2020, 8:24pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/51 "2020-09-07T20:24:55Z")
**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:
- Adaptive Sequential Monte Carlo for ABC targets based on [An adaptive sequential Monte Carlo method for approximate Bayesian computation | SpringerLink](https://link.springer.com/article/10.1007/s11222-011-9271-y)
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! 🙂
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [September 7, 2020, 9:16pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/52 "2020-09-07T21:16:11Z")
> [@francesco.alemanno](#):
>
> The package exports all symbols from `MonteCarloMeasurements` and `Distributions` , for ease of use
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
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [September 8, 2020, 4:05am UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/53 "2020-09-08T04:05:24Z")
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
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [September 11, 2020, 4:24pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/54 "2020-09-11T16:24:19Z")
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`
```julia
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)
```
with `KissABCv3.0`
```julia
res = smc(prior_dist,simulate; nparticles = 500,verbose = true,parallel = true,M = 10).P
```
using
```julia
approx_density = ApproxKernelizedPosterior(prior_dist,simulate,0.01)
res = sample(approx_density,AIS(30),MCMCThreads(),300,12;progress = true,ntransitions = 1600)
```
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
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [September 12, 2020, 8:25am UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/55 "2020-09-12T08:25:16Z")
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 😃 , 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
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [September 12, 2020, 5:17pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/56 "2020-09-12T17:17:10Z")
> 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:
```julia
(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`.
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [September 12, 2020, 7:32pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/57 "2020-09-12T19:32:04Z")
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)
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [September 12, 2020, 8:31pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/58 "2020-09-12T20:31:14Z")
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:
```julia
(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:
```julia
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.
---
### Author:  [@francesco.alemanno](https://discourse.julialang.org/u/francesco.alemanno)
#### Post date: [September 12, 2020, 8:45pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/59 "2020-09-12T20:45:17Z")
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.
---
### Author:  [@pjentsch0](https://discourse.julialang.org/u/pjentsch0)
#### Post date: [September 13, 2020, 6:13pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/60 "2020-09-13T18:13:45Z")
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.
---
### Author:  [@mlanghinrichs](https://discourse.julialang.org/u/mlanghinrichs)
#### Post date: [June 2, 2021, 2:59pm UTC](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668/61 "2021-06-02T14:59:49Z")
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](https://github.com/francescoalemanno/KissABC.jl/issues/30)). Glad for any help!
[Previous page](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668.md?page=2)
[Next page](https://discourse.julialang.org/t/ann-rfc-kissabc-jl-approximate-bayesian-computation/40668.md?page=4)