search

joshspeagle / dynesty

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

Slice sampling scaling/updates #260

Closed segasai closed 3 years ago

segasai commented 3 years ago

While debugging some issue with slice sampling I've encountered this reference https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4063214/ (as well as 2002.06212 ) which suggests that the slice sampling step should be scaled by 2*nexpand/(nexpand+ncontract)

While the dynesty code https://github.com/joshspeagle/dynesty/blob/b3dbd52369a9b9172b20e76cb53ddd7f43bf9a1c/py/dynesty/nestedsamplers.py#L249 scales by nexpand/(2 ncontract)

Is there any reason for the nonstandard scaling ?

joshspeagle commented 3 years ago

I hadn't done a ton of digging into scaling heuristics at the time that I chose that (apparently nonstandard) strategy, so there isn't too much reason behind it. Happy to change to the more commonly used approach if it looks like it performs better.

segasai commented 3 years ago

I've tried making a change and testing on Rosenbrock or Funnel distributions without an obvious improvement in number of evaluations. So lets keep things as they are unless there is a proof it's better.

segasai commented 3 years ago

Since I'm constantly seeing warnings like this: /home/skoposov/pyenv38/lib/python3.8/site-packages/dynesty/sampling.py:446: UserWarning: The slice sampling interval is longer than the cube size

I looked further at this question, and I think the strong rationale to change the current scaling from nexpand/2ncontract to nexpand/(nexpand+ncontract) is that The latter is much less noisy estimate of the expand fraction if nexpand or ncontract is low I.e. 

In [49]: N0=20;N1=np.random.binomial(N0,0.5,size=100000);N2=N0-N1
    ...: 

In [50]: def mad(x):
    ...:     return np.nanmedian(np.abs(x-np.nanmedian(x)))
    ...: 

In [51]: estim1=np.log10(N1/(N1+N2))
<ipython-input-51-f10221a845f5>:1: RuntimeWarning: divide by zero encountered in log10
  estim1=np.log10(N1/(N1+N2))

In [52]: estim2=np.log10(N1/(2*N2))
<ipython-input-52-188308011057>:1: RuntimeWarning: divide by zero encountered in log10
  estim2=np.log10(N1/(2*N2))

In [54]: print(np.nanmedian(estim1),mad(estim1))
-0.3010299956639812 0.07918124604762483

In [55]: print(np.nanmedian(estim2),mad(estim2))
-0.3010299956639812 0.17609125905568127

Because of this we are much less likely to see really large or small expansions.

joshspeagle commented 3 years ago

Should now be resolved as of #319.