joshspeagle / dynesty

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

Using dynesty with npdim #456

Closed MichaelDAlbrow closed 9 months ago

MichaelDAlbrow commented 10 months ago

Dynesty version dynesty 2.1.2 installed with pip

Describe the bug

This may be a bug or a fundamental misunderstanding on my part about the dimensions of the input and output of prior_transform. Please let me know if this is the case, and my advance apologies for raising this as a bug. (Michael Albrow)

I am trying to use dynesty with npdim < ndim. My prior transform accepts npdim parameters and returns ndim parameters.

_initialize_live_points in dynamicsampler.py has a call:

cur_live_v = M(prior_transform, np.asarray(cur_live_u))

or

cur_live_v = map(prior_transform, np.asarray(cur_live_u))

followed a little later by

live_v[ngoods:ngoods + nextra] = cur_live_v[cur_ind]

Earlier in the function, live_v is initialised with length npdim, but cur_live_v (being the output of the prior_transform) is of length ndim.

Dynesty output Traceback (most recent call last): File "/home/users/mda45/Microlensing/Galacticmodel_compress_10Nov21/ob170108.py", line 210, in sampler = dynesty.NestedSampler(loglike, prior_transform_total_samples, File "/home/users/mda45/miniconda3/lib/python3.10/site-packages/dynesty/dynesty.py", line 677, in new live_points, logvol_init, init_ncalls = _initialize_live_points( File "/home/users/mda45/miniconda3/lib/python3.10/site-packages/dynesty/dynamicsampler.py", line 458, in _initialize_live_p live_v[ngoods:ngoods + nextra] = cur_live_v[cur_ind] ValueError: could not broadcast input array from shape (227,10) into shape (227,4)

segasai commented 10 months ago

Thanks for the report.

I can say I never quite understood what npdim is supposed to do (what is the use-case) and there is not a single test written for that behaviour. So it is possible it got broken at some point.

Did this work in previous dynesty versions ?

This description " Number of parameters accepted by prior_transform. This might differ from ndim in the case where a parameter of loglikelihood is dependent upon multiple independently distributed parameters, some of which may be nuisance parameters." from dynesty.py doesn't really clarify much.

Can you comment on this @joshspeagle ?

MichaelDAlbrow commented 10 months ago

When I open dynamicsampler.py in github, and click on "blame", it seems that this code section may have been updated 10 months ago (line 422).

segasai commented 10 months ago

Yes, I know that code and it's recent, but I just have no clue if npdim stuff was working even before that change as there is zero testing of it.

MichaelDAlbrow commented 10 months ago

A quick and dirty solution that is working for me is to replace the line

live_v = np.zeros((nlive, npdim))

with

xx = prior_transform(live_u[0]) live_v = np.zeros((nlive,len(xx)))

joshspeagle commented 10 months ago

Oh geez, that's a real holdover.

Yes, that was originally implemented for a very particular edge case where you might simulate some parameters z ~ p(z) for the prior but then do some transformation/linear combination during/after the transform from u --> z --> x before feeding stuff into the likelihood L(x). (In all honesty, it can probably be deprecated at this point since the likelihood should probably implement that sort of transformation internally, keeping the number of dimensions fixed throughout.)

I can confirm this worked in very early versions (at least when dev testing), but I also have not tested this feature in a long time and would not at all be surprised if it is now broken given the drastic overhaul of much of the internals.

segasai commented 10 months ago

Ok. Thanks @joshspeagle . From your description it definitely sounds like this does need to be implemented inside dynesty. I will put then a deprecation warning when npdim=ndim and raise and error if npdim!=ndim.

segasai commented 10 months ago

Just an extra remark. It is certainly useful in some situations to return derived quantities from the prior and then save then with the posterior. At the moment that can be done with the 'blob' functionality and is possible for the likelihood function but not the prior itself. I think I could imagine implementing for the prior as well, but anyway that would need to be done differently from how npdim is treated.

segasai commented 9 months ago

I have merged the patch ( #457 ) that gets rid of broken npdim option. Thanks for the report. I'm closing this issue.