Walkers getting stuck outside uniform priors

[Nick Wogan]

Hi,

I'm using emcee to estimate atmospheric composition of a planet based on its spectrum. For these MCMC runs I use uniform priors for gas concentrations.In the prior I also say that sum(gas concentrations) > 1 has a zero probability (log(-infinity) probability). A sum larger than 1 would be unphysical.

In my experiments I find that many walkers often get stuck in part of the parameter space where sum(gas concentrations) > 1. I've attached an example of the walkers progress.

Looking for advice on how to prevent walkers from getting stuck in this unphysical region (sum(gas concetrations) > 1).

Thanks!

Nick

[Jo Bovy]

Hi Nick,

I think what you're doing is similar to sampling a simplex (the space of N numbers that sum to 1), but with the looser restriction that the sum has to be smaller than 1. I don't know whether there are standard ways of dealing with this, but there are standard changes of variables to sample the simplex. For example, I think stan uses the transformations presented in

https://arxiv.org/abs/1010.3436

I presented a similar transformation in appendix B of 

https://arxiv.org/abs/1704.03884

In your case, maybe you can re-parametrize your N concentration variables to be a) sum(concentrations) and b) the fractions concentration/sum(concentrations) and then transform a) such that it goes from -inf to inf (e.g., logit) and b) to be sampled as a simplex. Then you won't run into a hard boundary anywhere.

Best,

Jo

On Sun, Apr 10, 2022 at 9:20 PM Jo Bovy <bo...@astro.utoronto.ca> wrote:

Hi Nick,

I think what you're doing is similar to sampling a simplex (the space of N numbers that sum to 1), but with the looser restriction that the sum has to be smaller than 1. I don't know whether there are standard ways of dealing with this, but there are standard changes of variables to sample the simplex. For example, I think stan uses the transformations presented in

https://arxiv.org/abs/1010.3436

I presented a similar transformation in appendix B of 

https://arxiv.org/abs/1704.03884

In your case, maybe you can re-parametrize your N concentration variables to be a) sum(concentrations) and b) the fractions concentration/sum(concentrations) and then transform a) such that it goes from -inf to inf (e.g., logit) and b) to be sampled as a simplex. Then you won't run into a hard boundary anywhere.

Best,

Jo

--

Jo Bovy, University of Toronto

Web: https://astro.utoronto.ca/~bovy

Book: galaxiesbook.org

--
You received this message because you are subscribed to the Google Groups "emcee users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to emcee-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/ae5e3dfc-a4a0-48c4-a18b-321b76c3d07dn%40googlegroups.com.

[Dan Foreman-Mackey]

Great suggestions, Jo! In my experience, these kinds of transformations tend to be less important for emcee than they are for Stan, but still can be super helpful. 

Nick: if you truly are returning -np.inf from your logprob when the condition is not satisfied, then the walkers won't ever move there, so I expect the issue isn't actually what you think. I'll email some suggestions off the list so we don't add too much noise. 

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/CALeWRE4Q9PYfAkZVs_ZUAUyVy2FSH_%2BjFt26YVoe8LVH6aw_%3DA%40mail.gmail.com.

--

Dan Foreman-Mackey

Associate Research Scientist

Flatiron Institute

https://dfm.io

[Nick Wogan]

Not sure if my previous message to Jo isn't showing up. I might resend it here in a bit. 

Dan, I am returning -np.inf from my logprob. I guess the walkers are not ever moving outside of the uniform priors. But they get in a position where every move they attempt is outside of the priors. Thanks for any feedback you give via Email!

Nick

[Nick Wogan]

Response to Jo,

This is extremely helpful. Thanks so much!

I think what you are saying makes sense. I'll have to read Appendix B in your paper, and work through it a bit. But this overall seems like the right direction to take.

I'll report back to this thread once I have worked through the problem, and have a good understanding of it.

Thanks again,

Nick

[Nick Wogan]

Jo,

In your paper, Eq B2, says

y_i = logit(x_i) - logit(X_N)

Where X_N = [\frac{N-1}{N}, \frac{N-2}{N-1}, ..., \frac{1}{2}]

I am confused by B2, because X_N is a vector of length N-1, and x_i and y_i are scalars in vectors that are length N.

Can you clarify X_N in B2?

Nick

[Jo Bovy]

Hi Nick,

x_i in eq. (B1) is also a vector of length N-1 (that is, you only compute it for the first N-1 indices). My full code implementing everything in that appendix is here

https://github.com/jobovy/simple-m2m/blob/master/py/simplex.py

that should hopefully clarify everything (note that variables have different names).

I do agree with Dan that just returning -inf should allow the algorithm to work without the simplex transformation, I've done this many times. Given that the simplex is convex, for every proposal step, there should be many possibilities that return a valid set of concentrations. However, because the stretch move extends out to twice the distance between points, I guess that's where your chains go into territory that they shouldn't. Maybe you can try fiddling with the 'a' parameter in the stretch move, making it smaller so there is less overshooting into -inf territory? Others on this list might be able to comment on whether this is a good idea or not, but I think it's at least technically allowed.

Best,

Jo

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/1be9a030-9eb8-45fc-be9d-83cdf7ef9d1cn%40googlegroups.com.

[Nick Wogan]

Jo,

I tried several `a` parameters between 1 and 2. Resulted in similar behavior. MCMC started fine, but then "flatlined", where all moves led to sum(mixing ratios) > 1 regions of parameter space.

Your original suggestion of a variable transformation sounds promising. But would be happy to hear if there is an easier solution.

Nick

[Dan Foreman-Mackey]

Hi all,

For those following along, Nick and I chatted off-list a bit and it looks like there was a separate bug that was causing the prior not to be enforced, so these coordinate transforms may not actually be necessary after all.

Dan

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/d3cf1330-5c9f-4f56-93a0-149b18a997bfn%40googlegroups.com.