joshspeagle / dynesty

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

Issue with max. likelihood in corner #185

Closed JohannesBuchner closed 4 years ago

JohannesBuchner commented 4 years ago

Hi,

I set up the following toy test problem for a systematic evaluation of sampling methods across many types of problems I am making. The likelihood is a product of beta functions with random parameters. I am having issues with dynesty though, even when ndim=2. Is the issue that the ellipsoid wanders outside the cube? It stalls completely.

def transform(x):
    return x

rng = np.random.RandomState(ndim)
a_values = 10**rng.uniform(-1, 1, size=ndim)
b_values = 10**rng.uniform(-1, 1, size=ndim)

def loglike(theta):
    logprob = scipy.stats.beta.logpdf(theta, a=a_values, b=b_values)
    logprob[~(logprob > -1e300)] = -1e300
    return logprob.sum()

from dynesty import DynamicNestedSampler
sampler = DynamicNestedSampler(loglike, transform, ndim)
sampler.run_nested()

Here is the full output:

  self.expand_tot = self.vol_tot / vol_tot_orig
/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py:183: RuntimeWarning: invalid value encountered in double_scalars
  f = np.exp((np.log(vol) - np.log(self.vol)) / self.n)  # linear factor
/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py:401: RuntimeWarning: invalid value encountered in double_scalars
  self.expand_tot *= vol_tot / self.vol_tot
/home/user/.local/lib/python3.6/site-packages/dynesty/utils.py:43: RuntimeWarning: invalid value encountered in greater
  return (np.all(u[nonbounded] > 0.) and
10756it [00:40, 587.64it/s, batch: 0 | bound: 33 | nc: 4 | ncall: 43409 | eff(%): 24.496 | loglstar:   -inf < 16.300 <    inf | logz: -0.223 +/-  0.110 | dlogz:  0.643 >  0.010]

10756it [2:22:52,  1.25it/s, batch: 0 | bound: 33 | nc: 4 | ncall: 43409 | eff(%): 24.496 | loglstar:   -inf < 16.300 <    inf | logz: -0.223 +/-  0.110 | dlogz:  0.643 >  0.010]
Traceback (most recent call last):
  File "problems.py", line 233, in <module>
    run_sampler(paramnames, loglike, transform=transform, vectorized=vectorized)
  File "/mnt/data/daten/PostDoc2/home/Downloads/autoemcee/autosampler.py", line 74, in run_sampler
    sampler.run_nested(maxcall=max_ncalls)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/dynamicsampler.py", line 1625, in run_nested
    live_points=live_points):
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/dynamicsampler.py", line 841, in sample_initial
    n_effective=n_effective)):
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/sampler.py", line 782, in sample
    u, v, logl, nc = self._new_point(loglstar_new, logvol)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/sampler.py", line 380, in _new_point
    u, v, logl, nc, blob = self._get_point_value(loglstar)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/sampler.py", line 364, in _get_point_value
    self._fill_queue(loglstar)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/sampler.py", line 333, in _fill_queue
    point, axes = self.propose_point()
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/nestedsamplers.py", line 580, in propose_unif
    u, idx, q = self.mell.sample(rstate=self.rstate, return_q=True)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 454, in sample
    x = self.ells[0].sample(rstate=rstate)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 234, in sample
    return self.ctr + self.randoffset(rstate=rstate)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 218, in randoffset
    return np.dot(self.axes, randsphere(self.n, rstate=rstate))
  File "<__array_function__ internals>", line 6, in dot
KeyboardInterrupt

I pressed Ctrl-C after a while, the stack trace perhaps is telling where it is spending its time.

I also got a linalg error with ndim=10:

  warnings.warn("Random walk proposals appear to be "
/home/user/.local/lib/python3.6/site-packages/dynesty/sampling.py:216: UserWarning: Random number generation appears to be extremely inefficient. Adjusting the scale-factor accordingly.
  warnings.warn("Random number generation appears to be "
/home/user/.local/lib/python3.6/site-packages/dynesty/sampling.py:238: UserWarning: Random walk proposals appear to be extremely inefficient. Adjusting the scale-factor accordingly.
  warnings.warn("Random walk proposals appear to be "
23081it [12:25, 30.94it/s, batch: 0 | bound: 653 | nc: 25 | ncall: 1773987 | eff(%):  1.301 | loglstar:   -inf < 34.066 <    inf | logz: -3.775 +/-  0.242 | dlogz:  0.029 >  0.010]  
Traceback (most recent call last):
  File "problems.py", line 233, in <module>
    run_sampler(paramnames, loglike, transform=transform, vectorized=vectorized)
  File "/mnt/data/daten/PostDoc2/home/Downloads/autoemcee/autosampler.py", line 74, in run_sampler
    sampler.run_nested(maxcall=max_ncalls)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/dynamicsampler.py", line 1625, in run_nested
    live_points=live_points):
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/dynamicsampler.py", line 841, in sample_initial
    n_effective=n_effective)):
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/sampler.py", line 758, in sample
    bound = self.update(pointvol)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/nestedsamplers.py", line 566, in update
    pool=pool)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 585, in update
    firstell = bounding_ellipsoid(points, pointvol=pointvol)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 1363, in bounding_ellipsoid
    ell = Ellipsoid(ctr, covar)
  File "/home/user/.local/lib/python3.6/site-packages/dynesty/bounding.py", line 152, in __init__
    self.axes = lalg.cholesky(cov, lower=True)  # transformation axes
  File "/home/user/.local/lib/python3.6/site-packages/scipy/linalg/decomp_cholesky.py", line 91, in cholesky
    check_finite=check_finite)
  File "/home/user/.local/lib/python3.6/site-packages/scipy/linalg/decomp_cholesky.py", line 40, in _cholesky
    "definite" % info)
numpy.linalg.LinAlgError: 10-th leading minor of the array is not positive definite
joshspeagle commented 4 years ago

What seems to be happening is the bounding ellipsoid becomes ill-defined. Here's what I get grabbing the final bound right before it crashes the same way:

bound = sampler.bound[-1]  # grab multi-ellipsoid

bound.ctrs, bound.covs  # get covariance

Out: 
array([[5.68727758e-01, 6.20765700e-10]]),
array([[[ 3.40732147e-01, -4.60818321e-10],
        [-4.60818321e-10,  2.58193405e-18]]])

Right before this happens, it has a bunch of ellipsoids to try and capture this structure:

sampler.bound[-2].ctrs, sampler.bound[-2].covs

Out:
array([[7.78994760e-01, 7.76837638e-10],
       [6.04067203e-01, 8.66522033e-10],
       [5.38103294e-01, 1.24924083e-09],
       [4.50995375e-01, 1.17161581e-09],
       [3.26508871e-01, 1.39851199e-09],
       [2.09271961e-01, 1.50004789e-09],
       [1.19143973e-01, 2.09419876e-09],
       [5.31522393e-02, 3.27813293e-09],
       [2.15845285e-02, 3.57887206e-09],
       [9.31883901e-03, 4.08545452e-09],
       [1.90845769e-03, 6.43768057e-09]]),
array([[[ 4.06616474e-02, -5.18704975e-11],
        [-5.18704975e-11,  1.04618778e-18]],

       [[ 1.71662201e-03,  5.10527174e-12],
        [ 5.10527174e-12,  1.91615551e-18]],

       [[ 9.08793200e-04, -1.22183574e-11],
        [-1.22183574e-11,  9.21219120e-19]],

       [[ 5.25601725e-03,  1.49128177e-11],
        [ 1.49128177e-11,  3.37010043e-18]],

       [[ 6.51900859e-03, -1.62753168e-11],
        [-1.62753168e-11,  3.46138065e-18]],

       [[ 4.30822456e-03,  2.91625038e-12],
        [ 2.91625038e-12,  4.61720812e-18]],

       [[ 4.95258689e-03,  2.00712150e-11],
        [ 2.00712150e-11,  7.56857005e-18]],

       [[ 6.83041254e-04, -3.15470053e-12],
        [-3.15470053e-12,  9.91778647e-18]],

       [[ 1.60135247e-04, -1.89601649e-11],
        [-1.89601649e-11,  2.85046645e-17]],

       [[ 2.51907790e-05, -1.16499007e-11],
        [-1.16499007e-11,  3.68064311e-17]],

       [[ 1.11738885e-05, -2.11335828e-11],
        [-2.11335828e-11,  1.65686054e-16]]])

So it looks like it's a failure in the way the bounds are being dealt with as the axes become very skewed. I experimented with adding in reflective or periodic boundary conditions, but it still fails. So looks like issues with the clustering/bounding strike again! (sigh...)

JohannesBuchner commented 4 years ago

https://johannesbuchner.github.io/UltraNest/ integrates this without problems. One trick I use there is to only accept an update when the new bounding geometry is smaller in volume and did not produce numerical errors. Perhaps such a check would work here too, if numerical issues are the cause.

joshspeagle commented 4 years ago

Yea, that's not too surprising -- the bounding distributions are definitely the least stable part of dynesty and have issues that have plagued the code for years. I have a bunch of checks I tried to throw in to make the bounds more robust but they don't really seem to have fixed these numerical issues. I probably should add some additional logic similar to yours to finally try and quash these. I'll mark it as a to-do item.

joshspeagle commented 4 years ago

Closing this since it should now be resolved per #195. Thanks @JohannesBuchner!!!