The main work horse of pymc3 and many more modern Bayesian library is the no-u-turn sampler. Internally, pymc3 build the probabilistic model (everything in the `with pm.Model() ...`

context) in theano by compiling a callable that you can pass a 1d array (all the unknown/free parameters) and get back a scalar log_prob and the gradient (a 1d array). Then NUTS (or any other sampler) takes callable and the input and do the magic in numpy + multiprocess.

It would be great to isolate this part of the code into a standalone library, it will help us testing and benchmarking.