From d208e40a141b8128457c6e5ffadae14c34f1fd7e Mon Sep 17 00:00:00 2001 From: Shanika Date: Tue, 21 Jan 2020 18:35:25 -0800 Subject: [PATCH 1/7] initial commit --- gwpopulation/hyperpe.py | 61 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 58 insertions(+), 3 deletions(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index 91eadaed..31c0f924 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -41,14 +41,14 @@ class HyperparameterLikelihood(Likelihood): """ def __init__(self, posteriors, hyper_prior, sampling_prior=None, - ln_evidences=None, max_samples=1e100, - selection_function=lambda args: 1, + ln_evidences=None, max_samples=1e100, pastro=None, + fiducial_vt=None, selection_function=lambda args: 1, conversion_function=lambda args: (args, None), cupy=True): """ Parameters ---------- posteriors: list - An list of pandas data frames of samples sets of samples. + A list of pandas data frames of samples sets of samples. Each set may have a different size. These can contain a `prior` column containing the original prior values. @@ -61,6 +61,11 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, of the hyperparameter likelihood. If not provided, the original evidences will be set to 0. This produces a Bayes factor between the sampling power_prior and the hyperparameterised model. + pastro: array + An array of pastro values corresponding to the event posteriors + fiducial_vt: float + The visible spacetime volume described by the fiducial model + used to calculate pastro. selection_function: func Function which evaluates your population selection function. conversion_function: func @@ -109,6 +114,19 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, self.samples_factor =\ - self.n_posteriors * np.log(self.samples_per_posterior) + if pastro is not None: + if fiducial_vt is not None: + logger.info('Fiducial VT set to {}'.format(fiducial_vt)) + self.fiducial_vt = fiducial_vt + else: + logger.info('No fiducial VT provided, defaulting to 0.005') + self.fiducial_vt = 0.005 + if len(pastro) != len(posteriors): + raise ValueError( + 'Number of pastro values provided are not + equal to the number of posteriors') + self.pastro = pastro + def log_likelihood_ratio(self): self.parameters, added_keys = self.conversion_function(self.parameters) self.hyper_prior.parameters.update(self.parameters) @@ -276,3 +294,40 @@ def _get_selection_factor(self): def generate_rate_posterior_sample(self): pass + +class PastroLikelihood(HyperparameterLikelihood): + """ + A mixture model likelihood which utlises p_astro for the effective + noise likelihood. + This likelihood is marginalised over rate. + + See Eq. (40) of https://arxiv.org/abs/1912.09708 for a definition. + """ + def log_likelihood_ratio(self): + self.parameters, added_keys = self.conversion_function(self.parameters) + self.hyper_prior.parameters.update(self.parameters) + ln_l1 = self._compute_per_event_ln_bayes_factors() + ln_l1 += self._get_selection_factor + ln_l1 += self._get_fiducial_vt_factor() + + if added_keys is not None: + for key in added_keys: + self.parameters.pop(key) + + ln_l2 = self._get_pastro_factor + ln_l2 = xp.nan_to_num(ln_l2) + + ln_l = xp.sum(xp.logaddexp(ln_l1,ln_l2)) + ln_l += 0.5 * self._get_selection_factor() + + return float(xp.nan_to_num(ln_l)) + + def _get_selection_factor(self): + return - xp.log( + self.selection_function(self.parameters)) + + def _get_pastro_factor(self): + return xp.log() + + def _get_fiducial_vt_factor(self): + return xp.log(self.fiducial_vt) From 39a915c94bc9f6f24eb0c617af58ab8888ae7c94 Mon Sep 17 00:00:00 2001 From: Shanika Date: Tue, 21 Jan 2020 21:43:37 -0800 Subject: [PATCH 2/7] add pastro calculation --- gwpopulation/hyperpe.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index 70330ae5..eea32657 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -325,7 +325,8 @@ def _get_selection_factor(self): self.selection_function(self.parameters)) def _get_pastro_factor(self): - return xp.log() + pastro_factor = (1.0 - self.pastro)/self.pastro + return xp.log(pastro_factor) def _get_fiducial_vt_factor(self): return xp.log(self.fiducial_vt) From 35afeb9e5b16ffdfa6df13b05367e5dd7ea09b8a Mon Sep 17 00:00:00 2001 From: Shanika Date: Tue, 21 Jan 2020 22:08:47 -0800 Subject: [PATCH 3/7] fix code climate issues --- gwpopulation/hyperpe.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index eea32657..7ac0a659 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -115,14 +115,13 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, if pastro is not None: if fiducial_vt is not None: logger.info('Fiducial VT set to {}'.format(fiducial_vt)) - self.fiducial_vt = fiducial_vt + self.fiducial_vt = fiducial_vt else: logger.info('No fiducial VT provided, defaulting to 0.005') self.fiducial_vt = 0.005 if len(pastro) != len(posteriors): - raise ValueError( - 'Number of pastro values provided are not - equal to the number of posteriors') + raise ValueError('Number of pastro values provided are not ' + 'equal to the number of posteriors.') self.pastro = pastro def log_likelihood_ratio(self): @@ -293,6 +292,7 @@ def _get_selection_factor(self): def generate_rate_posterior_sample(self): pass + class PastroLikelihood(HyperparameterLikelihood): """ A mixture model likelihood which utlises p_astro for the effective @@ -315,7 +315,7 @@ def log_likelihood_ratio(self): ln_l2 = self._get_pastro_factor ln_l2 = xp.nan_to_num(ln_l2) - ln_l = xp.sum(xp.logaddexp(ln_l1,ln_l2)) + ln_l = xp.sum(xp.logaddexp(ln_l1, ln_l2)) ln_l += 0.5 * self._get_selection_factor() return float(xp.nan_to_num(ln_l)) From 2631e1405ecfe8c7599e9fe8c837738b294347a9 Mon Sep 17 00:00:00 2001 From: Shanika Date: Wed, 22 Jan 2020 15:08:58 -0800 Subject: [PATCH 4/7] fix typo --- gwpopulation/hyperpe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index 7ac0a659..10daabb2 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -305,7 +305,7 @@ def log_likelihood_ratio(self): self.parameters, added_keys = self.conversion_function(self.parameters) self.hyper_prior.parameters.update(self.parameters) ln_l1 = self._compute_per_event_ln_bayes_factors() - ln_l1 += self._get_selection_factor + ln_l1 += self._get_selection_factor() ln_l1 += self._get_fiducial_vt_factor() if added_keys is not None: From 043c3ff573120ebda1aa55a03a7e3185bd537adf Mon Sep 17 00:00:00 2001 From: Shanika Date: Wed, 22 Jan 2020 17:49:34 -0800 Subject: [PATCH 5/7] rename function --- gwpopulation/hyperpe.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index 10daabb2..f35dce79 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -42,7 +42,7 @@ class HyperparameterLikelihood(Likelihood): def __init__(self, posteriors, hyper_prior, sampling_prior=None, ln_evidences=None, max_samples=1e100, pastro=None, - fiducial_vt=None, selection_function=lambda args: 1, + fiducial_selection=None, selection_function=lambda args: 1, conversion_function=lambda args: (args, None), cupy=True): """ Parameters @@ -63,7 +63,7 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, the sampling power_prior and the hyperparameterised model. pastro: array An array of pastro values corresponding to the event posteriors - fiducial_vt: float + fiducial_selection: float The visible spacetime volume described by the fiducial model used to calculate pastro. selection_function: func @@ -113,12 +113,12 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, self.n_posteriors = len(posteriors) if pastro is not None: - if fiducial_vt is not None: - logger.info('Fiducial VT set to {}'.format(fiducial_vt)) - self.fiducial_vt = fiducial_vt + if fiducial_selection is not None: + logger.info('Fiducial VT set to {}'.format(fiducial_selection)) + self.fiducial_selection = fiducial_selection else: logger.info('No fiducial VT provided, defaulting to 0.005') - self.fiducial_vt = 0.005 + self.fiducial_selection = 0.005 if len(pastro) != len(posteriors): raise ValueError('Number of pastro values provided are not ' 'equal to the number of posteriors.') @@ -301,12 +301,18 @@ class PastroLikelihood(HyperparameterLikelihood): See Eq. (40) of https://arxiv.org/abs/1912.09708 for a definition. """ + def __init__(self, posteriors, hyper_prior, sampling_prior=None, + ln_evidences=None, max_samples=1e100, pastro=None, + fiducial_selection=None, selection_function=lambda args: 1, + conversion_function=lambda args: (args, None), cupy=True): + + def log_likelihood_ratio(self): self.parameters, added_keys = self.conversion_function(self.parameters) self.hyper_prior.parameters.update(self.parameters) ln_l1 = self._compute_per_event_ln_bayes_factors() ln_l1 += self._get_selection_factor() - ln_l1 += self._get_fiducial_vt_factor() + ln_l1 += self._get_fiducial_selection_factor() if added_keys is not None: for key in added_keys: @@ -328,5 +334,5 @@ def _get_pastro_factor(self): pastro_factor = (1.0 - self.pastro)/self.pastro return xp.log(pastro_factor) - def _get_fiducial_vt_factor(self): - return xp.log(self.fiducial_vt) + def _get_fiducial_selection_factor(self): + return xp.log(self.fiducial_selection) From 6ad67dcfbd414103de53ca5ef508817173c23533 Mon Sep 17 00:00:00 2001 From: Shanika Date: Wed, 22 Jan 2020 20:54:08 -0800 Subject: [PATCH 6/7] edit functions --- gwpopulation/hyperpe.py | 56 ++++++++++++++++++++++------------------- 1 file changed, 30 insertions(+), 26 deletions(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index f35dce79..f6b89bf8 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -41,8 +41,8 @@ class HyperparameterLikelihood(Likelihood): """ def __init__(self, posteriors, hyper_prior, sampling_prior=None, - ln_evidences=None, max_samples=1e100, pastro=None, - fiducial_selection=None, selection_function=lambda args: 1, + ln_evidences=None, max_samples=1e100, + selection_function=lambda args: 1, conversion_function=lambda args: (args, None), cupy=True): """ Parameters @@ -61,11 +61,9 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, of the hyperparameter likelihood. If not provided, the original evidences will be set to 0. This produces a Bayes factor between the sampling power_prior and the hyperparameterised model. - pastro: array - An array of pastro values corresponding to the event posteriors fiducial_selection: float The visible spacetime volume described by the fiducial model - used to calculate pastro. + determined by the sampling prior distribution. selection_function: func Function which evaluates your population selection function. conversion_function: func @@ -112,23 +110,11 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, self.n_posteriors = len(posteriors) - if pastro is not None: - if fiducial_selection is not None: - logger.info('Fiducial VT set to {}'.format(fiducial_selection)) - self.fiducial_selection = fiducial_selection - else: - logger.info('No fiducial VT provided, defaulting to 0.005') - self.fiducial_selection = 0.005 - if len(pastro) != len(posteriors): - raise ValueError('Number of pastro values provided are not ' - 'equal to the number of posteriors.') - self.pastro = pastro - def log_likelihood_ratio(self): self.parameters, added_keys = self.conversion_function(self.parameters) self.hyper_prior.parameters.update(self.parameters) ln_l = xp.sum(self._compute_per_event_ln_bayes_factors()) - ln_l += self._get_selection_factor() + ln_l += self.n_posteriors * self._get_selection_factor() if added_keys is not None: for key in added_keys: self.parameters.pop(key) @@ -146,7 +132,7 @@ def _compute_per_event_ln_bayes_factors(self): self.sampling_prior, axis=-1)) def _get_selection_factor(self): - return - self.n_posteriors * xp.log( + return - xp.log( self.selection_function(self.parameters)) def generate_extra_statistics(self, sample): @@ -303,9 +289,31 @@ class PastroLikelihood(HyperparameterLikelihood): """ def __init__(self, posteriors, hyper_prior, sampling_prior=None, ln_evidences=None, max_samples=1e100, pastro=None, - fiducial_selection=None, selection_function=lambda args: 1, - conversion_function=lambda args: (args, None), cupy=True): - + selection_function=lambda args: 1, + conversion_function=lambda args: (args, None), + fiducial_selection=None, cupy=True): + """ + Parameters + ---------- + pastro: array + An array of pastro values corresponding to the event posteriors. + fiducial_selection: float + The visible spacetime volume described by the fiducial model + determined by the sampling prior distribution. + """ + if pastro is not None: + if fiducial_selection is not None: + logger.info('Fiducial VT set to {}'.format(fiducial_selection)) + self.fiducial_selection = fiducial_selection + else: + logger.warning('No fiducial VT provided, defaulting to 1.0') + self.fiducial_selection = 1.0 + if len(pastro) != len(posteriors): + raise ValueError('Number of pastro values provided are not ' + 'equal to the number of posteriors.') + self.pastro = pastro + + super().__init__(self) def log_likelihood_ratio(self): self.parameters, added_keys = self.conversion_function(self.parameters) @@ -326,10 +334,6 @@ def log_likelihood_ratio(self): return float(xp.nan_to_num(ln_l)) - def _get_selection_factor(self): - return - xp.log( - self.selection_function(self.parameters)) - def _get_pastro_factor(self): pastro_factor = (1.0 - self.pastro)/self.pastro return xp.log(pastro_factor) From b1b30c987d1233275ed36b522986a3986aecbf9c Mon Sep 17 00:00:00 2001 From: Shanika Date: Wed, 22 Jan 2020 20:57:11 -0800 Subject: [PATCH 7/7] remove parameter comment string --- gwpopulation/hyperpe.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/gwpopulation/hyperpe.py b/gwpopulation/hyperpe.py index f6b89bf8..d2ce2c40 100644 --- a/gwpopulation/hyperpe.py +++ b/gwpopulation/hyperpe.py @@ -61,9 +61,6 @@ def __init__(self, posteriors, hyper_prior, sampling_prior=None, of the hyperparameter likelihood. If not provided, the original evidences will be set to 0. This produces a Bayes factor between the sampling power_prior and the hyperparameterised model. - fiducial_selection: float - The visible spacetime volume described by the fiducial model - determined by the sampling prior distribution. selection_function: func Function which evaluates your population selection function. conversion_function: func