By “non-analytical” and “numerical”, do you mean it requires a non-standard embedded iterative solver of some kind (i.e., not just an ODE solver, root finder, 1D integrator, etc.)? Even if you need an iterative algorithm, If the energy is continuously differentiable and can be defined in terms of elementary distributions and transforms and matrix operations, you should be able to code it in PyMC and let that take care of the derivatives.