Want to fit mixture of multivariate normals; maybe Ben and Marcus have some thoughts on this one

[Andrew Gelman]

Hi all. I was talking with Sebastian about fitting mixture of multivariate normals—it comes up for our hierarchical EP algorithm. Basically what we’re looking for is a little program that automatically does this, probably an R function that calls Stan. Here are the two challenges:
1. Avoiding degenerate solutions where one mixture component has zero variance.
2. Handling the aliasing, which requires identifying which mixture is which.

My idea is as follows:

1. Using a default weakly informative prior that constrains the relative values of the variances. I’m not sure how best to do this. A fairly lightweight (from a statistical perspective) approach would be to put a prior on the determinants of the cov matrix, but I’m not sure how to easily do this in Stan so that it wouldn’t be computationally expensive. Another approach would be to use one of our existing models such as inv-Wishart or LKJ with some hyperparameter that itself has a proper prior. Really there are three concerns here: first is the very basic issue of the prior constraining things to be away from the boundary; second is, if we are using a hierarchical model here, this could improve the fit by regularizing the cov matrices (under the reasonable assumption that they’re likely to be of comparable scale; third is that we want the whole thing to run fast so we want to avoid models that are expensive to evaluate.

2. How to handle the aliasing? In 1-d it’s simple enough because you can just constrain the means of the mixture components to be ordered: mU_1 < mu_2 < mu_3 etc. We could do the same in higher dimensions, using for example the first dimension and constraining those means to be increasing—but this could cause problems, for example, if there are 2 components whose means are close in dimension 1 but far apart in dimension 2. So here’s another idea: We first attempt to locate the K components using some crude algorithm such as nearest neighbor, so that we get K points in D-dimensional space. We then identify the model based on the distance of the modes from these points. I’m still not exactly sure how to do this: if K=2 it’s easy, we could just constrain the two modes based on the difference of the distances between the 2 points. (Actually, I’m not sure if it’s so easy to apply that constraint in Stan.) In higher dimensions it would be easier to use a soft constraint that penalizes each mode based on its distance from the crude mode. But that’s kinda sloppy because then we’re using the data twice and not quite fitting the model we want to fit.

The other issue is that we’d like the program to do full Bayes or optimize. I’m hoping that running stan on optimize setting will be easier in the new rstan that’s coming soon? I don’t like that “optimizing” function; I’d rather either call stan() with an optimize=TRUE setting or have a function called stan_optimize(). Also I would like the optimization to work with the stan() function that we have at that web-location, the version that has the make-like behavior.

Anyway, if we’re optimizing, things are simpler. We still have to deal with issue 1 above, but issue 2 is less of a concern because we can just pick any starting point and then it should find one of the modes.

See you
A

[Ben Goodrich]

On Thursday, May 7, 2015 at 9:38:59 AM UTC-4, Andrew Gelman wrote:

Here are the two challenges:
1.  Avoiding degenerate solutions where one mixture component has zero variance.

We need to be more specific about communicating this problem to users. It is pretty easy to ensure that a marginal variance is positive. The main difficulty is conditional variances going to zero, particularly when optimizing.

1.  Using a default weakly informative prior that constrains the relative values of the variances.  I’m not sure how best to do this.  A fairly lightweight (from a statistical perspective) approach would be to put a prior on the determinants of the cov matrix, but I’m not sure how to easily do this in Stan so that it wouldn’t be computationally expensive.  

Anything parameterized with (large) covariance matrices is expensive (and not so accurate numerically). Anything parameterized with Cholesky factors is cheap and pretty accurate. That includes (log) determinants.
 

Another approach would be to use one of our existing models such as inv-Wishart or LKJ with some hyperparameter that itself has a proper prior.

I think you want to avoid putting infinite density on singular matrices. This is kind of hard to do numerically with the inverse Wishart because the determinant term and the trace term pull in different directions as the covariance matrix approaches being singular. For LKJ it is simply a matter of constraining the shape parameter to be > 1. Also, the decomposition into standard deviations and a correlation matrix makes it easier to deal with the ordering below.

2.  How to handle the aliasing?  In 1-d it’s simple enough because you can just constrain the means of the mixture components to be ordered:  mU_1 < mu_2 < mu_3 etc.  We could do the same in higher dimensions, using for example the first dimension and constraining those means to be increasing—but this could cause problems, for example, if there are 2 components whose means are close in dimension 1 but far apart in dimension 2.  

It could be subject to the same problems, but I would try constraining the standard deviations to be increasing over the mixture components. You could do this for one variable or compose the standard deviations for each mixture as a simplex multiplied by a scale parameter and make the scale parameters to be a pos_ordered vector.

The other issue is that we’d like the program to do full Bayes or optimize.  I’m hoping that running stan on optimize setting will be easier in the new rstan that’s coming soon?  I don’t like that “optimizing” function; I’d rather either call stan() with an optimize=TRUE setting or have a function called stan_optimize().  Also I would like the optimization to work with the stan() function that we have at that web-location, the version that has the make-like behavior.

They take different settings, so it is a bit harder to add a optimize argument to stan(). But if having an alias for optimize() called stan_optimize() is worth enough for you to answer user's questions as to which one to use, I'll put it in.
 

Anyway, if we’re optimizing, things are simpler.  We still have to deal with issue 1 above, but issue 2 is less of a concern because we can just pick any starting point and then it should find one of the modes.

In sense that issue 2 is less of a problem, yes. But issue 1 is less simple under optimization because the optimum can correspond to a singular covariance matrix.

Ben

[Andrew Gelman]

On May 7, 2015, at 12:22 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Thursday, May 7, 2015 at 9:38:59 AM UTC-4, Andrew Gelman wrote:

Here are the two challenges:
1.  Avoiding degenerate solutions where one mixture component has zero variance.

We need to be more specific about communicating this problem to users. It is pretty easy to ensure that a marginal variance is positive. The main difficulty is conditional variances going to zero, particularly when optimizing.

1.  Using a default weakly informative prior that constrains the relative values of the variances.  I’m not sure how best to do this.  A fairly lightweight (from a statistical perspective) approach would be to put a prior on the determinants of the cov matrix, but I’m not sure how to easily do this in Stan so that it wouldn’t be computationally expensive.  

Anything parameterized with (large) covariance matrices is expensive (and not so accurate numerically). Anything parameterized with Cholesky factors is cheap and pretty accurate. That includes (log) determinants.

OK, then a hierarchical prior on the determinants (for example, lognormal) with a proper prior on the hyperparameters, that would do it.

 

Another approach would be to use one of our existing models such as inv-Wishart or LKJ with some hyperparameter that itself has a proper prior.

I think you want to avoid putting infinite density on singular matrices.

Of course!  The whole idea is to put a prior on that has zero density on singular matrices!  Maybe that wasn’t clear in my email.

This is kind of hard to do numerically with the inverse Wishart because the determinant term and the trace term pull in different directions as the covariance matrix approaches being singular. For LKJ it is simply a matter of constraining the shape parameter to be > 1. Also, the decomposition into standard deviations and a correlation matrix makes it easier to deal with the ordering below.

2.  How to handle the aliasing?  In 1-d it’s simple enough because you can just constrain the means of the mixture components to be ordered:  mU_1 < mu_2 < mu_3 etc.  We could do the same in higher dimensions, using for example the first dimension and constraining those means to be increasing—but this could cause problems, for example, if there are 2 components whose means are close in dimension 1 but far apart in dimension 2.  

It could be subject to the same problems, but I would try constraining the standard deviations to be increasing over the mixture components. You could do this for one variable or compose the standard deviations for each mixture as a simplex multiplied by a scale parameter and make the scale parameters to be a pos_ordered vector.

That’s a thought.  My gut is that organizing by location will be more stable than organizing by scale, but we could try it your way and see if it works.

The other issue is that we’d like the program to do full Bayes or optimize.  I’m hoping that running stan on optimize setting will be easier in the new rstan that’s coming soon?  I don’t like that “optimizing” function; I’d rather either call stan() with an optimize=TRUE setting or have a function called stan_optimize().  Also I would like the optimization to work with the stan() function that we have at that web-location, the version that has the make-like behavior.

They take different settings, so it is a bit harder to add a optimize argument to stan(). But if having an alias for optimize() called stan_optimize() is worth enough for you to answer user's questions as to which one to use, I'll put it in.

I think the problem is that optimize() is already an R function.  But, as Bob says, optimizing() is just a confusing name.

 

Anyway, if we’re optimizing, things are simpler.  We still have to deal with issue 1 above, but issue 2 is less of a concern because we can just pick any starting point and then it should find one of the modes.

In sense that issue 2 is less of a problem, yes. But issue 1 is less simple under optimization because the optimum can correspond to a singular covariance matrix.

That’s not a problem, we can use the boundary-avoiding Wishart (not inverse-Wishart) priors from this paper:

http://www.stat.columbia.edu/~gelman/research/published/chung_cov_matrices.pdf

A