Re-estimating a probability distribution with additional priors

Question

I have a 3D dataset with at least millions of data points (scatter events from atoms, approximately Gaussian). I am modeling this data with a Gaussian Mixture Model. The usual approach would be to optimize the GMM with the EM algorithm. However, I have substantial priors (e.g., the data is produced by a polymer made of light atoms) so the only way to find the correct GMM is via sampling.

Here's the problem: it's very expensive to compare a GMM to the data $D$, and I have to do that evaluation frequently. I would like to introduce an intermediate step: first model $G_d$ without priors, and then introduce priors and re-model $G_m$ where both are GMMs. Then while sampling $G_m$ I could just compare it to $G_d$ instead of $D$, which is a lot cheaper. So $G_d$ is kind of a data reduction. Intuitively the likelihood function for $G_m$ should be:

$$ p(D|G_m)\approx p(D|G_d)S(G_d|G_m) $$

where $S(G_d|G_m)$ is some kind of scoring function (maybe K-L divergence? correlation?). My questions are:

• Can you re-model in this way and still get a good approximation of the likelihood?
• What types of functions can be used for $S(G_d|G_m)$ and still get a proper PDF?
• Does the choice of number of components in $G_d$ affect the final likelihood?

Below I attempt to illustrate this problem. The magenta $G_d$ is computed at the beginning. Then the red $G_m$ with its additional priors (including size and linearity) are fit to the magenta one, much faster than comparing to the data itself. Quick clarification: the number of Gaussians in $G_m$ is much higher than $G_d$. That's because the only prior on $G_d$ is the Dirichlet one (so it's just regular Bayesian GMM fitting) whereas $G_m$ has much more detailed priors which enables us to make a more detailed model than the data itself would indicate.

Answer 1

It sounds like you want to reduce computation by summarizing your data and then running inference on the small summary. This is known as data squashing and there are many variants. The main principle is that you want the summary to capture as much information in the data as possible, in other words, you want it to overfit. Thus a Bayesian GMM is the wrong thing to use as a summary. A greedy algorithm like k-means or density tree would be a better choice.

Similarly, there is no reason to constrain $G_d$ to have a smaller number of clusters than $G_m$. If you do, then you will have to solve an expensive disaggregation problem when fitting $G_m$. Since the whole point of the method is to save computation, it is better to use a cheap over-clustering of the data as the summary. Then while fitting $G_m$ you can constrain the data in each region to have the same cluster assignment, giving a large speedup. For details, see the paper Scalable model-based clustering for large databases based on data summarization.

If instead $G_d$ is part of your sensor model, then the appropriate formula is: $$ p(D|G_m) = \int_{G_d} p(D|G_d) p(G_d|G_m) dG_d $$ where $p(G_d|G_m)$ is a physical model of the "effective" scatterers given the atomic positions. The method in the question corresponds to approximating this integral with a single point. A more accurate result, with similar computational complexity, could be obtained by Laplace's method or variational approximation.

Answer 2

One way of looking at all this is as introducing an intermediate variable $G_d$ from which your data $D$ are drawn. This gives rise to the same integral as @TomMinka gave

$$ p(D|G_m) = \int_{G_d} p(D|G_d) p(G_d|G_m) dG_d. $$

An alternative is to transform/aggregate/summarize your data to some $f(D)$ that is easier to compute with, giving a much simpler question of modeling/maximizing $p(f(D)|G_m)$. The particular example I have in mind is simply binning the data, supposing you know a decent feature size. You could get fancy and bin adaptively too if you wanted. It is an extremely simple way you could make the computation cheaper.

I'm not sure off the top of my head how hard it would be to get a decent approximation of $p(N_\textrm{counts in bin}|G_m)$, but I don't know if there's an analytic form in 3D. This book seems to talk about it, but I don't have access. If it's still too expensive, you might get a decent approximation from pretending all the points in the bin were at its center. That will probably give you more computational savings per loss of information with adaptive binning.

And of course you can pair aggregated/summarized data with an intermediate variable to get both kinds of savings and extra flexibility via

$$ p(f(D)|G_m) = \int_{G_d} p(f(D)|G_d) p(G_d|G_m) dG_d. $$