Just to clarify the issue of integrating external code into a pymc3 model (which is different from the question about likelihood-free approaches): I think this is an entirely valid and useful thing to do. But for larger models (say more than ~20 parameters) you pretty much need the gradients. Metropolis just won’t cut it. How hard they are to come by can be very different depending on the external code, and you also need some familiarity with theano and often quite a bit of (multivariate) calculus to get that to work. We are definitely missing good examples for this.