I recently coauthored this paper applying PyMC3 in a population genetic context:
Bayesian Inference of Joint Coalescence Times for Sampled Sequences
A feature is the use of the CategoricalGibbsMetropoplis step method to sample from extremely large spaces of matrices representing the combinatorial properties of genealogical/coalescent trees. In effect, this is equivalent to sampling from the set of all permutations of a large set.
It also uses the StickBreaking transform to generate prior distributions on a simplex, which are a second order approximation to some sampled distribution.
Code is here