Bespoke MCMC priors & likelihoods, & feeding a posterior joint pdf back in as the prior next time

Question

We're looking at PyStan, PyMC3 and emcee. (switching to R could also be an option, if need be).

We have a lot of bespoke priors and bespoke likelihood functions: they are bespoke in the sense that they are not one of the standard distributions: they are unique to our work. These are used to estimate the joint distribution of several real parameters; and ideally we'd like the estimator to spit out a posterior in a format that it could later read in as a prior for further estimation: e.g. as an empirical distribution function (that would be then interpolated when used as a prior), perhaps with something to ensure the tails go out appropriately far.

So, I need an expert compare-and-contrast: is it easier to do that with one of those libraries (PyStan, PyMC3, emcee) rather than the others, and what additional programming languages are needed / would be helpful to do so? (e.g. AIUI it's possible to extend STAN with additional C++ programming)

Answer 1

The BayesianTools R package (disclaimer: I am one of the authors) allows specifying arbitrary priors, and there is a function createPriorDensity for summarising the posterior as a new prior.

The limitation of this is that we currently only implement a multivariate normal density estimator, so you will have a loss of information if your posterior is not approximately multivariate normal, which is likely the case if you have weak data. Extending the function to include more flexible density estimators (e.g. gaussian processes) is on our todo list, but it's quite tricky to get this stable, so I'm not sure when we will have this working. In general, creating reliable empirical density estimates in high-dimensional parameter spaces is extremely tricky, regardless of the software you are using.

If possible, I would therefore recommend

1. For a few updates to just accept the computational costs and compute the model again with the full / updated data

2. For frequent updates, consider using SMC instead of MCMC sampling, which doesn't require a density estimate for updating the posterior with new data (the BT package has a rudimentary SMC sampler on CRAN, and a more evolved sampler in a development branch)

Answer 2

See this previous answer of mine to another question: In Stan is there a way to use parameter posterior from old analysis as prior in new analysis? There I mentioned the idea to approximate the posterior samples by some distribution - if different parameters are correlated this might have to be some multivariate distribution.

If one goes down that route, then a general purpose MCMC sampler that let's you specify your own prior distributions without any real restrictions (such as the PyStan python package or the R rstan package, which only limit you by not allowing discrete prior parameters) is a good options. In that case no matter what you approximate your posterior by, you can specify it as a prior.

In terms of fitting some multivariate mixture distribution to posterior samples, there is a considerable difference in the quality of different implementations in terms of numerical stability and speed. E.g. in R the mclust package (page 53 here has an example of using mclust for a bivariate normal approximation to a posterior) or the RBeST package have good algorithms at least for the bivariate normal mixture case that you could use or take a look at. I would assume other programming languages also have good packages for this sort of thing, but I do not know them.

Answer 3

A few things to note (I could be miss-interpreting your problem so please correct me if I am taking this in the wrong direction):

Define your prior (before seeing any data) as $\pi_0(\theta)$. Define a first data set as $Y_1$ and a second data set as $Y_2$. Then define the first posterior $$ \pi_1(\theta|Y_1) = \frac{f(Y_1|\theta)\pi_0(\theta)}{m(Y_1)} $$

You can simulate from $\pi_1(\theta|Y_1)$ with STAN, pyMC, etc. (the tools you mention)

To restate what you said in your question, you want to use $\pi_1(\theta|Y_1)$ as the prior for $\pi_2(\theta|Y_2)$ i.e.

$$ \pi_2(\theta|Y_2) = \frac{f(Y_2|\theta)\pi_1(\theta|Y_1)}{m(Y_2)} $$ but notice $$ \frac{f(Y_2|\theta)\pi_1(\theta|Y_1)}{m(Y_2)}= \frac{f(Y_2|\theta)f(Y_1|\theta)\pi_0(\theta)}{m(Y_1)m(Y_2)} $$

So under the assumption $Y_1$ and $Y_2$ are independent (this generally goes along with iid assumptions, and if it wasn't true you would still use the final result below)

$$ \pi_2(\theta|Y_2)=\pi_2(\theta|Y_2,Y_1) = \frac{f(Y_2,Y_1|\theta)\pi_0(\theta)}{m(Y_2,Y_1)} $$

Which can be simulated using STAN, pyMC, etc.

My point: You don't have to empirically estimate the prior with the previous posterior, You can just run the model over again with all the data

I understand that the above recommendation may not be possible in some situations, it could be computationally costly or the data may not be in a convenient form it. In these situations you can use a self-normalized importance sampler (which can be implemented in Python, R, or something similar).

So let's back-up and suppose you drew $\theta^{(1)},...,\theta^{(G)} \sim \pi_1(\theta|Y_1)$. Now you want to sample from $\pi_2(\theta|Y_2,Y_1)$. To do this you can re-weight $\theta^{(1)},...,\theta^{(G)}$ with the importance weights $$ w^{(g)} = \frac{f(Y_2|\theta)}{\sum_{g=1}^G f(Y_2|\theta) } $$ so they effectively become a sample from $\pi_2(\theta|Y_2,Y_1)$.

For importance sampling to work, $\pi_1(\theta|Y_1)$, and $\pi_2(\theta|Y_2,Y_1)$ must be similar. If not, the variance of the sampler will be very high. You can approximate the "effective sample size" of the importance sampler via; $$ G.\mathrm{eff} \approx \frac{1}{\sum_{g=1}^G [w^{(g)}]^2} $$ There are other related techniques that reduce the variance of the procedure (i.e. sequential importance sampling/particle filtering and rejection sampling). It would be a good idea for you to research these algorithms on your own and decide for yourself how you would like to use them if you are interested in this type pof approach.

Also see @Xian answer to my own question Does this Monte Carlo Technique Have a Name? ...he knows more about this stuff than I do.