I’m sure this question is anathema to many of you.

Is there a way to obtain maximum likelihood estimates, or even just the likelihood given some parameter values, of a model in PyMC3?

For context, I have an idea about how to model some data which works as a Bayesian hierarchical model in PyMC3. However, I work in a field that is still very much dominated by *p* values, and my concern is that readers, and even some of my coauthors, will not accept this model if it is Bayesian, at least not initially. I would like the initial presentation of the model to be as accessible as possible for my audience, and unfortunately that means using MLE and *p* values.

I know how to code this model in PyMC3 but not outside of it (e.g., using just SciPy).