Symmetry of stretch move proposal

[Hendrik Jennen]

Hi all,

I have been trying to understand the symmetry condition for the stretch move as stated in the Goodman and Weare paper from 2010, and implemented as such in de emcee library, see also my question here.

In summary, I would expect the stated condition of g(1/z) = z g(z) to be replaced by g(1/z) = z^2 g(z).

I would much appreciate any insights on my (mis)understandings here!

KR

Hendrik

[David W Hogg]

What’s your argument? It might help to understand that.

--
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/0993b8a0-0447-4026-ab90-f291de290739n%40googlegroups.com.

--

David W. Hogg — http://cosmo.nyu.edu/hogg/

[Dan Foreman-Mackey]

Looking at the stack exchange question, the main issue here is a misunderstanding of what the symmetry condition is. One qualitative way to think about it is: "the probability of proposing X_k starting from Y (but using the same proposal distribution) should be equal to the probability of proposing Y starting from X_k". What you call the "inverse proposal" in your stack exchange question never enters! Hopefully this sets you on the right track.

Best,

Dan

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/CAEZJrfo%2BVEDxJwkF-cVG2P3y4Xj925u7%3DT9uNrJK9nDPwQ_5KQ%40mail.gmail.com.

--

Dan Foreman-Mackey

Associate Research Scientist

Flatiron Institute

https://dfm.io

[Hendrik Jennen]

Hi David and Dan,

Thanks for your replies!

I have reconsidered it but still end up with the power of 2 instead of 1. I'll write out the derivation here; hopefully we find where my error slips in.

The proposal X_k -> Y = X_i + z (X_k - X_i) has a probability

\int_a^b dy g(y)

where [a, b] is some small interval around z. If we use the same proposal distribution g(y) to go back from Y to X_k = X_i + (1 / z) (Y - X_i) , this integral should equal

\int_{1 / b}^{1 / a} dy g(y),

since this implies that the probabilities to go from X_k to Y and vice versa are equal. The second integral we write as (use y -> 1/y)

\int_b^a (-) dy y^{-2} g(1 / y) = \int_a^b dy y^{-2} g(1 / y)

leading to

g(1 / y) = y^2 g(y).

Another way to see this is that an infinitesimal region [y, y+dy] with volume dy and probability g(y) dy transforms into a region [1 / (y+dy), 1/ y] with volume dy / y^2 and probability g(1 / y) dy / y^2.

Could you help me out where this goes wrong, or show what the correct derivation is to get to the g(1 / y) = y g(y) condition?

Thanks in advance,

Hendrik

[David W Hogg]

I don’t have a full answer, but there is an issue with this derivation, which is that the moves are in the X space not the y space, so the intervals to be compared are in X units not y units. That is, we are comparing densities in the X space. The reverse move has a different dX/dy than the forward move. I’ll try to derive something later. The thing that amazes me is that the dimension of the affine space doesn’t enter. It really feels like it should.

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/6eae36ab-5307-4baa-be50-2b6a43cf1cean%40googlegroups.com.

[Will M. Farr]

Riffing off David's message (and maybe saving him a derivation): at the point X_k and proposing a jump to Y via 

Y = X_i + z (X_k - X_i)

the proposal density in Y space (also X space!) is 

p(Y|X_k) dY = p(z) dz

but z is acting as a dimensionless radial coordinate for Y about the "origin" X_i, so p(Y|X_k) dY is 

p(Y | X_k) dY = p(Y | X_k) r^(n-1) dr dOmega = p(Y | X_k) z^(n-1) |X_k - X_i|^n dz dOmega

(where dOmega is the "angular" part of the Y volume element---the d(cos(theta))d(phi) in 3D, for example) whence 

p(Y|X_k) = p(z) / (z^(n-1) |X_k - X_i|^n)

If we are now at the point Y, and are proposing jumps to Y' via

Y' = X_i + z' (Y - X_i)

the density is again 

p(Y'|Y) = p(z') dz'/dY ~ p(z') / (z'^(n-1) |Y - X_i|^n)

In order for the jump to return to the same place, we must have Y' = X_k, which is achieved when z' = 1/z.  Note that |Y-X_i|^n = z^n |X_k - X_i|^n.  Putting all that together, the correct factor that goes into the Metropolis-Hastings acceptance probability for the jump proposal ratio when a jump is proposed from X_k -> Y is 

p(Y' = X_k | Y) / p(Y | X_k) = p(z'=1/z)/p(z) z^(n-2)

Clearly you get the formula from the emcee paper if you choose p(z'=1/z)/p(z) = z, or p(z'=1/z) = z p(z), which is their "symmetric" condition.  I think what they meant by "symmetric" (this language also appears in Goodman and Weare) is that *if* you were in dimension 1, then the proposal ratio appearing in the M-H acceptance probability would be unity; their condition achieves this.  I think what *you* mean by "symmetric" is that there are equal densities in z' at z' = 1/z and in z at the point z; this is satisfied by a "flat in log" proposal, with p(z) = 1/z, or p(z'=1/z)/p(z) = z^2 => p(z'=1/z) = z^2 p(z), as you write.  In this case, the factor that appears in the M-H acceptance probability is z^n (not ^(n-1)!).

As long as you put the right factors of z^whatever into the M-H acceptance probability, you can do whatever you want.  Personally, I've observed a bit better performance with your version of symmetric; it is the unique density on z that, assuming z \in [1/a, a] with a > 1 (i.e. the default emcee a = 2), gives you an equal chance of moving *toward* X_i and *away* from X_i, which is a nice symmetry to have.  (For a while, perhaps now, this was the default proposal in `ptemcee`, for example.). But the performance change is not super large.

In any case, I guess that symmetric is in the eye of the beholder.  I hope this made things clear!

Cheers,

Will

P.S.---With apologies to Hogg!  I sent this just to him ("Reply") a moment ago, so now he gets two copies!  (Thanks for noticing, David!)

To view this discussion on the web visit https://groups.google.com/d/msgid/emcee-users/CAEZJrfrwoct%2BVC6AgX3xv4vnPcM-paqkX7VOcuN7%2BOBTP6ymug%40mail.gmail.com.

[Will M. Farr]

I suppose there is also the "symmetric" proposal that makes the proposal ratio in the M-H acceptance probability unity which satisfies p(z' = 1/z) = z^(2-n) p(z).  This gives p(z) = z^((n-2)/2) (compare to the other "symmetric" proposal in n = 1 dimension). Maybe that would be a good proposal to use---I don't know.  But it will very strongly prefer jumps *away* from X_i as the dimension gets larger.

Will

[Hendrik Jennen]

Hi Will,

Thanks a lot for your answer, and especially for showing how to derive the correct ratio of proposals in X-space; that explains a lot!

I indeed interpreted symmetric in 1 fixed way, only looking at the 1D proposal distribution, and did not consider the interplay with the acceptance ratio.

Cheers,

Hendrik