From ce573e6c366a7688af8b49ae1ba7ea8332e738ff Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Wed, 17 Oct 2018 23:11:21 +1100
Subject: [PATCH] add generation of distance posterior from marginalized
likelihood
---
bilby/gw/conversion.py | 47 +++++++++++++++++++++++++++++++++++++++++-
1 file changed, 46 insertions(+), 1 deletion(-)
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index ac40efe0c..1a0df0fe6 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -1,8 +1,9 @@
from __future__ import division
import numpy as np
+from pandas import DataFrame
from ..core.utils import logger, solar_mass
-from ..core.prior import DeltaFunction
+from ..core.prior import DeltaFunction, Interped
try:
from astropy.cosmology import z_at_value, Planck15
@@ -621,6 +622,11 @@ def _generate_all_cbc_parameters(sample, defaults, base_conversion,
output_sample, _ = base_conversion(output_sample)
output_sample = generate_mass_parameters(output_sample)
output_sample = generate_spin_parameters(output_sample)
+ if likelihood is not None:
+ if likelihood.distance_marginalization:
+ output_sample = \
+ generate_distance_samples_from_marginalized_likelihood(
+ output_sample, likelihood)
output_sample = generate_source_frame_parameters(output_sample)
compute_snrs(output_sample, likelihood)
return output_sample
@@ -947,3 +953,42 @@ def compute_snrs(sample, likelihood):
else:
logger.debug('Not computing SNRs.')
+
+
+def generate_distance_samples_from_marginalized_likelihood(sample, likelihood):
+ if isinstance(sample, dict):
+ pass
+ elif isinstance(sample, DataFrame):
+ for ii in range(len(sample)):
+ temp = _generate_distance_sample_from_marginalized_likelihood(
+ dict(sample.iloc[ii]), likelihood)
+ sample['luminosity_distance'][ii] = temp['luminosity_distance']
+ return sample
+
+
+def _generate_distance_sample_from_marginalized_likelihood(sample, likelihood):
+ signal_polarizations = \
+ likelihood.waveform_generator.frequency_domain_strain(sample)
+ rho_mf_sq = 0
+ rho_opt_sq = 0
+ for ifo in likelihood.interferometers:
+ signal = ifo.get_detector_response(signal_polarizations, sample)
+ rho_mf_sq += ifo.matched_filter_snr_squared(signal=signal)
+ rho_opt_sq += ifo.optimal_snr_squared(signal=signal)
+
+ rho_mf_sq_dist = \
+ rho_mf_sq * sample['luminosity_distance'] / \
+ likelihood._distance_array
+
+ rho_opt_sq_dist = \
+ rho_opt_sq * sample['luminosity_distance']**2 / \
+ likelihood._distance_array**2
+
+ distance_log_like = (rho_mf_sq_dist.real - rho_opt_sq_dist.real / 2)
+
+ distance_post = np.exp(distance_log_like - max(distance_log_like)) *\
+ likelihood.distance_prior_array
+
+ sample['luminosity_distance'] = Interped(
+ likelihood._distance_array, distance_post).sample()
+ return sample
--
GitLab