joshspeagle / dynesty

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

Ellipsoid fix #219

Closed segasai closed 3 years ago

segasai commented 3 years ago

Hi,

This fixes https://github.com/joshspeagle/dynesty/issues/205 The first patch actually addresses the root cause. When you scale the ellipsoid to contain all the points, you inadequately correct for the roundoff errors.

See for example this where (similarly to bounding_ellipsoid()) function I determine the maximum norm to then scale the matrix such that the norm =1 The actual numbers deviate from 1 by up to 1e-5 for poorly conditioned matrices.

In [74]: delts=np.zeros(10000)
    ...: for i in range(len(delts)):
    ...:     ndim = 30
    ...:     df=30
    ...:     xs=np.random.normal(size=ndim)
    ...:     mat=scipy.stats.wishart(df=df,scale=np.ones(ndim)).rvs()
    ...:     val= xs@scipy.linalg.pinv(mat)@xs
    ...:     delts[i]=(xs@scipy.linalg.pinv(mat*val)@xs)
    ...: print (np.max(np.abs(delts-1)))
    ...: 
    ...: 
    ...: 
3.5509573712566578e-06

That in turn leads to points not belonging to ellipsoids and so forth.

I've addressed the true cause, but also made the code much more agrressive in error checking. Now the code will simply bail out if the points are outside an ellipsoid. I also removed a couple of blanket "try, catch"'es as that can just hide bugs.

Just in case this array below """pts=np.array([[0.74072701, 0.23989705, 0.33023426, 0.01215883, 0.51469837, 0.44890521, 0.77941745, 0.13952267, 0.35263657, 0.74846314, 0.70115809], [0.74165003, 0.23966594, 0.32937525, 0.01215881, 0.51136602, 0.44875993, 0.79001457, 0.15858188, 0.35211 , 0.74649555, 0.70140197], [0.74233832, 0.23939336, 0.32715769, 0.01202915, 0.50620443, 0.44851161, 0.77373995, 0.1382089 , 0.35128579, 0.74413639, 0.70130415], [0.74164991, 0.24003932, 0.32639726, 0.01209666, 0.50898113, 0.44849516, 0.78452441, 0.14235259, 0.34992684, 0.74262878, 0.70118466], [0.740734 , 0.23978596, 0.33174908, 0.01218793, 0.5178418 , 0.44902064, 0.76970458, 0.12575265, 0.35381555, 0.75058863, 0.70134179], [0.7413627 , 0.24014706, 0.32930086, 0.01216193, 0.50885556, 0.44870013, 0.76510886, 0.10653991, 0.35219345, 0.74414166, 0.70147609], [0.74136586, 0.24022307, 0.32878415, 0.01215922, 0.50828152, 0.4486667 , 0.76648194, 0.10463014, 0.35170491, 0.74305404, 0.70139347], [0.74246291, 0.23961114, 0.3239174 , 0.0119511 , 0.50187017, 0.44824451, 0.77924463, 0.12496938, 0.34893025, 0.73909268, 0.70081965], [0.74254531, 0.23948345, 0.32582755, 0.01199907, 0.5059507 , 0.44840498, 0.76854328, 0.11854014, 0.35026727, 0.74156595, 0.70119547], [0.74159331, 0.23977765, 0.32963856, 0.0121741 , 0.51518553, 0.44881403, 0.7744988 , 0.13457225, 0.35207124, 0.74506993, 0.70141684], [0.74126511, 0.24010425, 0.32773292, 0.01214011, 0.50901677, 0.44872399, 0.78655283, 0.1296933 , 0.35037293, 0.73994855, 0.70080575], [0.74119231, 0.23998482, 0.32951103, 0.01218396, 0.51324479, 0.44887061, 0.77898945, 0.13012461, 0.35164717, 0.74333866, 0.70109244], [0.74077727, 0.24024786, 0.33391715, 0.01233366, 0.51716158, 0.44911472, 0.76724817, 0.09321799, 0.35497593, 0.74826123, 0.7016588 ], [0.74166418, 0.23993247, 0.32590573, 0.0120518 , 0.50840416, 0.4484495 , 0.78806661, 0.1511295 , 0.34985441, 0.74231698, 0.70114253], [0.74241233, 0.2395396 , 0.32565818, 0.01200924, 0.50859475, 0.44838418, 0.77232206, 0.11893393, 0.35006803, 0.74193645, 0.70105261], [0.74100223, 0.24014051, 0.33376745, 0.01231775, 0.51387807, 0.44911211, 0.76502106, 0.08574744, 0.3549149 , 0.74741795, 0.70165675], [0.74099902, 0.24131237, 0.3114181 , 0.01171154, 0.47014757, 0.44745408, 0.86074656, 0.1370894 , 0.33862851, 0.71358677, 0.6980086 ], [0.74081512, 0.23939776, 0.33351629, 0.01219075, 0.51486742, 0.44926762, 0.77418365, 0.12992944, 0.35503526, 0.74890462, 0.70115509], [0.74231298, 0.23931151, 0.32800698, 0.01205665, 0.50939603, 0.44858468, 0.77498763, 0.1467109 , 0.35184026, 0.74587929, 0.70143051], [0.74265923, 0.23941667, 0.32541527, 0.01197518, 0.50639489, 0.44837442, 0.76457188, 0.11356411, 0.35005118, 0.74136728, 0.70117492], [0.74084963, 0.23977435, 0.33106627, 0.01216702, 0.514034 , 0.44896372, 0.77524194, 0.13096818, 0.35336353, 0.74891676, 0.70124602], [0.74236521, 0.23969941, 0.32358891, 0.01195424, 0.5059654 , 0.44821809, 0.7778404 , 0.12045847, 0.34859623, 0.73933573, 0.70078013], [0.74180653, 0.23964357, 0.32816407, 0.01210661, 0.51603759, 0.4486947 , 0.76440525, 0.11874132, 0.35119684, 0.74584656, 0.70125696], [0.74132644, 0.24009473, 0.32751781, 0.01213121, 0.51246674, 0.44864774, 0.78147269, 0.13673156, 0.35045549, 0.7426282 , 0.70112192]]) """ causes the unpatched dynesty dynesty.bounding.bounding_ellipsoid(points) to create ellipsoid that does not contain all the points.

joshspeagle commented 3 years ago

Fantastic patch and sleuthing for these bugfixes. Happy to merge this in.

segasai commented 3 years ago

Thanks for merging. It may be a good idea to add a few tests on this (ellipsoid boundaries) to the test suite. In the spirit of the script that I have there, generate random covar matrix, random points from it and then generate a bounding ellipsoid.

joshspeagle commented 3 years ago

Yea, that's a good idea. I'll open a new issue on this to flag it.