Bounding fixes

[segasai]

Following up on #425

I've stumbled upon this https://github.com/joshspeagle/dynesty/blob/2d49cf9104e417a44ada7f2f1d013cbfee0acf01/py/dynesty/nestedsamplers.py#L498

and similar code for MultiEllipsoid. Basically the default initialisation for these bounds is zero centre and unit covariance matrix which is obviously incorrect for an ellipsoid covering the whole cube. This is pretty inconsequential for slice/walk samplers, but this has many issues in particular for uniform sampler.

• That will sample a significant fraction of the volume outside the cube
• The sampling will be missing a significant volume and will be biased towards zero, at least till the first update of the boundary.

The correct initialisation should be centre at [0.5,...0.5] and covar of diag([ndim/4,...,ndim/4])

While I have a patch ready #428, I'm creating this issue for traceability (and feedback)

Comments ?

[joshspeagle]

This definitely looks like a problem! Glad the fix is simple at least.

[segasai]

I think there is more to it in #425 though with bound updates that is not quite right, I'm currently working on a possible patch/refactoring to get rid of the _beyond_unit_bound which was always a bit obscure. with the aim of clearly separating the state of the sampler, i.e. are we using the bound or not , (and that should be just a scalar variable) from decision of possible updating the bound, which should be just in one single location.

The reason why I started doing that, is because I still do not understand why default ellipsoid centers/covariances are used at all for sampling. I thought they would always be updated before the run based on current points.

[segasai]

While we on this @joshspeagle

I'm looking at this line: https://github.com/joshspeagle/dynesty/blob/b237972175009dda0fb5efda32871330a8428ea8/py/dynesty/sampler.py#L399

This lines implies that the boundary update will be done based on local ncall variable of _new_point rather than the global ncall of the sampling run. That is not correct, is it ?

There is also a check here: https://github.com/joshspeagle/dynesty/blob/b237972175009dda0fb5efda32871330a8428ea8/py/dynesty/sampler.py#L803 But this relies on since_update which was set inside _new_point.

[joshspeagle]

So the logic for both of these was a bit hacky, but essentially if you initialize your multi-ellipsoidal bounds too early over a uniform distribution you run the risk of shredding the posterior and making sampling inefficient/creating tons of islands. So the default behaviour was to wait a bit before the first update by just sampling from the cube.

In terms of the _beyond_unit_bound functionality, that was also hacky to try and just repeat sampling if we somehow end up proposing outside of the priors because of our bounds (which can happen depending on the expansion factor, etc.).

I still do not understand why default ellipsoid centers/covariances are used at all for sampling. I thought they would always be updated before the run based on current points.

I also thought this is what should be happening, so I was a bit surprised to see it wasn't triggering properly.

This lines implies that the boundary update will be done based on local ncall variable of _new_point rather than the global ncall of the sampling run. That is not correct, is it ?

So IIRC, one of these is because in earlier iterations where the bounding was unstable, you could propose a terrible bound and then you kinda were stuck. So that might have been included to trigger an pre-emptive update to avoid this from happening.

Not sure about the second check...

[segasai]

Thanks for your comments @joshspeagle . I'm not sure that quite answered my question )) But I have now updated the #428 patch with a refactor (still somewhat preliminary) that gets rid of of the _beyound_unit. Now there is a variable inside the class that is describing the current state of the sampler (unit cube sampling or not ). And then there is update_bound_if_needed() function that performs the bound update if it is required. The key difference (and that's where the previous code was incorrect I believe is that the ncall inside the _new_point() was local, preventing updates if nqueue was small enough I believe. But anyway I think the current logic is much clearer. And I finally understand when/what could happen. Let's see if that solves the #431 as well