joshspeagle / dynesty

Dynamic Nested Sampling package for computing Bayesian posteriors and evidences
https://dynesty.readthedocs.io/
MIT License
357 stars 77 forks source link

batch slice sampling stuck (No live points are above loglstar) #273

Closed segasai closed 3 years ago

segasai commented 3 years ago

Okay, I've still got an issue with dynamic slice sampling.

8971it [19:00,  7.87it/s, batch: 3 | bound: 283 | nc: 40 | ncall: 314184 | eff(%):  2.828 | loglstar: 64906.518 < 64911.523 < 64913.034 | logz: 64896.802 +/-  0.510 | stop:  1.482] 
6424024061129661440
Traceback (most recent call last):
 ,,,, line 421, in sampler_dyn
    dsampler.run_nested(
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py", line 1603, in run_nested
    passback = self.add_batch(nlive=nlive_batch,
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py", line 1723, in add_batch
    for results in self.sample_batch(nlive_new=nlive,
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py", line 1165, in sample_batch
    for it, results in enumerate(
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py", line 802, in sample
    u, v, logl, nc = self._new_point(loglstar_new, logvol)
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py", line 414, in _new_point
    u, v, logl, nc, blob = self._get_point_value(loglstar)
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py", line 397, in _get_point_value
    self._fill_queue(loglstar)
  File "/home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py", line 359, in _fill_queue
    raise RuntimeError('No live points are above loglstar. '
RuntimeError: No live points are above loglstar. Do you have likelihood plateau ? 
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
.....

--> 421         dsampler.run_nested(
    422             dlogz_init=.1,
    423             nlive_init=nlive,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in run_nested(self, nlive_init, maxiter_init, maxcall_init, dlogz_init, logl_max_init, n_effective_init, nlive_batch, wt_function, wt_kwargs, maxiter_batch, maxcall_batch, maxiter, maxcall, maxbatch, n_effective, stop_function, stop_kwargs, use_stop, save_bounds, print_progress, print_func, live_points)
   1601                     # Compute our sampling bounds using the provided
   1602                     # weight function.
-> 1603                     passback = self.add_batch(nlive=nlive_batch,
   1604                                               wt_function=wt_function,
   1605                                               wt_kwargs=wt_kwargs,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in add_batch(self, nlive, wt_function, wt_kwargs, maxiter, maxcall, logl_bounds, save_bounds, print_progress, print_func, stop_val)
   1721                 if logl_bounds is None:
   1722                     logl_bounds = wt_function(res, wt_kwargs)
-> 1723                 for results in self.sample_batch(nlive_new=nlive,
   1724                                                  logl_bounds=logl_bounds,
   1725                                                  maxiter=maxiter,

~/pyenv38/lib/python3.8/site-packages/dynesty/dynamicsampler.py in sample_batch(self, nlive_new, update_interval, logl_bounds, maxiter, maxcall, save_bounds)
   1163         # sample past the original bounds "for free".
   1164         for i in range(1):
-> 1165             for it, results in enumerate(
   1166                     self.sampler.sample(dlogz=0.,
   1167                                         logl_max=logl_max,

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in sample(self, maxiter, maxcall, dlogz, logl_max, n_effective, add_live, save_bounds, save_samples)
    800             # `logl > loglstar` using the bounding distribution and sampling
    801             # method from our sampler.
--> 802             u, v, logl, nc = self._new_point(loglstar_new, logvol)
    803             ncall += nc
    804             self.ncall += nc

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _new_point(self, loglstar, logvol)
    412         while True:
    413             # Get the next point from the queue
--> 414             u, v, logl, nc, blob = self._get_point_value(loglstar)
    415             ncall += nc
    416 

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _get_point_value(self, loglstar)
    395         # If the queue is empty, refill it.
    396         if self.nqueue <= 0:
--> 397             self._fill_queue(loglstar)
    398 
    399         # Grab the earliest entry.

~/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py in _fill_queue(self, loglstar)
    357             args = (np.nonzero(self.live_logl > loglstar)[0], )
    358             if len(args[0]) == 0:
--> 359                 raise RuntimeError('No live points are above loglstar. '
    360                                    'Do you have likelihood plateau ? ')
    361         else:

RuntimeError: No live points are above loglstar. Do you have likelihood plateau ? 
> /home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampler.py(359)_fill_queue()
    357             args = (np.nonzero(self.live_logl > loglstar)[0], )
    358             if len(args[0]) == 0:
--> 359                 raise RuntimeError('No live points are above loglstar. '
    360                                    'Do you have likelihood plateau ? ')
    361         else:

The diagnosis reveals the following problem. See the plot

xx

The black points show u in one dimension vs likelihood in saved_run, and red points show stuff from current new_run. Basically what seems to happen -- when the points are sampled for the batch and then the sampling continues from them, they miss the narrow posterior mode, so they are stuck in a lower mode, and therefore they can't reach the actual mode. Because the batch sampling is set up in such a way to not stop until we hit logl_max, it never finishes till all logl are constant and slice sampling can't go on.

I'm not sure there is an easy fix here. But in the same time, the failure doesn't sound so unlikely to me, if we have a wide mode and a bit higher narrow mode. So it's probably worth trying to do something about it.

segasai commented 3 years ago

One thought I had is to maybe not run this https://github.com/joshspeagle/dynesty/blob/2c70dbee7ac3868d251bdc7a3e99d565ae11506f/py/dynesty/dynamicsampler.py#L1166 with dlogz=0 but small dlogz. (1e-3,1e-5 or something). We can then warn a user if logl_max wasn't reached, but still proceed.

joshspeagle commented 3 years ago

That is an excellent idea! Great find. Yes, setting dlogz=1e-3 and then throwing a warning if the logl_max wasn't reached should be a perfect solution for this problem in general.