Skip to content

Resolve "Add histogram distribution as option in the hierarchical analysis"

Matthew Pitkin requested to merge hist_hierarch into master

Resolves #45 (closed).

This adds a non-parametric histogram-style distribution to the hierarchical analysis options.

Currently the tests pass, but trying to use this for an analysis fails (e.g., using test result in '/home/matt/cwinpy_test/pe/outdir') with:

from bilby.core.result import ResultList
from cwinpy.hierarchical import MassQuadrupoleDistribution

data = ResultList(["gaussian_example_result.json"])

data[0].posterior["Q22"] = data[0].posterior["x"]
data[0].priors["Q22"] = data[0].priors["x"]

distkwargs = {"low": data[0].posterior["Q22"].min(), "high": data[0].posterior["Q22"].max(), "nbins": 10}
intmethod = "numerical"
sampler_kwargs = {
    "sample": "unif",
    "nlive": 500,
    "gzip": True,
    "outdir": "exponential_distribution",
    "label": intmethod,
    "check_point_plot": False,
}

mqd = MassQuadrupoleDistribution(
    data=data,
    distribution="histogram",
    distkwargs=distkwargs,
    sampler_kwargs=sampler_kwargs,
    integration_method=intmethod,
    nsamples=500,  # number of Q22 samples to store/use
)

res = mqd.sample()

gives:

19:06 bilby INFO    : Running for label 'numerical', output will be saved to 'exponential_distribution'
19:06 bilby WARNING : Parameter weight has no default prior and is set to None, this will not be sampled and may cause an error.
19:06 bilby INFO    : Search parameters:
19:06 bilby INFO    :   weight0 = Dirichlet(order=0, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight1 = Dirichlet(order=1, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight2 = Dirichlet(order=2, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight3 = Dirichlet(order=3, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight4 = Dirichlet(order=4, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight5 = Dirichlet(order=5, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight6 = Dirichlet(order=6, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight7 = Dirichlet(order=7, n_dimensions=10, label='weight')
19:06 bilby INFO    :   weight8 = Dirichlet(order=8, n_dimensions=10, label='weight')
19:06 bilby INFO    : Single likelihood evaluation took 3.556e-04 s
0it [00:00, ?it/s]19:06 bilby INFO    : Using sampler Dynesty with kwargs {'bound': 'multi', 'sample': 'unif', 'verbose': True, 'periodic': None, 'reflective': None, 'check_point_delta_t': 600, 'nlive': 500, 'first_update': None, 'walks': 100, 'npdim': None, 'rstate': None, 'queue_size': 1, 'pool': None, 'use_pool': None, 'live_points': None, 'logl_args': None, 'logl_kwargs': None, 'ptform_args': None, 'ptform_kwargs': None, 'enlarge': 1.5, 'bootstrap': None, 'vol_dec': 0.5, 'vol_check': 8.0, 'facc': 0.2, 'slices': 5, 'update_interval': 300, 'print_func': <bound method Dynesty._print_func of <bilby.core.sampler.dynesty.Dynesty object at 0x7fad2e9e8c40>>, 'dlogz': 0.1, 'maxiter': None, 'maxcall': None, 'logl_max': inf, 'add_live': True, 'print_progress': True, 'save_bounds': False, 'n_effective': None, 'maxmcmc': 5000, 'nact': 5}
19:06 bilby INFO    : Checkpoint every check_point_delta_t = 600s
19:06 bilby INFO    : Using dynesty version 1.0.1
19:06 bilby INFO    : Resume file exponential_distribution/numerical_resume.pickle does not exist.
19:06 bilby INFO    : Generating initial points from the prior
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-99-62c69578641e> in <module>
----> 1 res = mqd.sample()

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/cwinpy/hierarchical.py in sample(self, **run_kwargs)
   1673             )
   1674         else:
-> 1675             self._result = bilby.run_sampler(
   1676                 likelihood=self._likelihood,
   1677                 priors=self._prior,

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/sampler/__init__.py in run_sampler(likelihood, priors, label, outdir, sampler, use_ratio, injection_parameters, conversion_function, plot, default_priors_file, clean, meta_data, save, gzip, result_class, npool, **kwargs)
    182         result = sampler._run_test()
    183     else:
--> 184         result = sampler.run_sampler()
    185     end_time = datetime.datetime.now()
    186 

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/sampler/dynesty.py in run_sampler(self)
    354             if self.kwargs['live_points'] is None:
    355                 self.kwargs['live_points'] = (
--> 356                     self.get_initial_points_from_prior(self.kwargs['nlive'])
    357                 )
    358             self.sampler = dynesty.NestedSampler(

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/sampler/base_sampler.py in get_initial_points_from_prior(self, npoints)
    440         while len(unit_cube) < npoints:
    441             unit = np.random.rand(self.ndim)
--> 442             theta = self.prior_transform(unit)
    443             if self.check_draw(theta, warning=False):
    444                 unit_cube.append(unit)

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/sampler/dynesty.py in prior_transform(self, theta)
    689 
    690         """
--> 691         return self.priors.rescale(self._search_parameter_keys, theta)
    692 
    693 

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/prior/dict.py in rescale(self, keys, theta)
    688             required_variables = {k: result[k] for k in getattr(self[key], 'required_variables', [])}
    689             result[key] = self[key].rescale(theta[index], **required_variables)
--> 690         return [result[key] for key in keys]
    691 
    692     def _update_rescale_keys(self, keys):

~/miniconda3/envs/cwinpy/lib/python3.8/site-packages/bilby/core/prior/dict.py in <listcomp>(.0)
    688             required_variables = {k: result[k] for k in getattr(self[key], 'required_variables', [])}
    689             result[key] = self[key].rescale(theta[index], **required_variables)
--> 690         return [result[key] for key in keys]
    691 
    692     def _update_rescale_keys(self, keys):

I'm not sure if this is due to the implementation in CWInPy, or just an issue with bilby.

Edited by Matthew Pitkin

Merge request reports