From c5bfa29d547dc8ed804480f3452b3f34030c88d3 Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Mon, 21 Sep 2020 17:14:43 -0500
Subject: [PATCH] Fix bug in combining evidences
If `use_ratio` was None, the `log_evidence` was recalculated correctly,
but the log_bayes_factor was not updated. This meant the
log_bayes_factor from the first result file was stored causing
inconsistencies (the log_ev - log_ev_err != log_bf). This fixes this
issue and makes the interface more straight forward, always using the
log_evidence directly
---
bilby/core/result.py | 20 +++++++++++---------
1 file changed, 11 insertions(+), 9 deletions(-)
diff --git a/bilby/core/result.py b/bilby/core/result.py
index df91621d5..d604e52d1 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -1644,24 +1644,26 @@ class ResultList(list):
The result object with the combined evidences.
"""
self.check_nested_samples()
- if result.use_ratio:
- log_bayes_factors = np.array([res.log_bayes_factor for res in self])
- result.log_bayes_factor = logsumexp(log_bayes_factors, b=1. / len(self))
- result.log_evidence = result.log_bayes_factor + result.log_noise_evidence
- result_weights = np.exp(log_bayes_factors - np.max(log_bayes_factors))
- else:
- log_evidences = np.array([res.log_evidence for res in self])
- result.log_evidence = logsumexp(log_evidences, b=1. / len(self))
- result_weights = np.exp(log_evidences - np.max(log_evidences))
+
+ # Combine evidences
+ log_evidences = np.array([res.log_evidence for res in self])
+ result.log_evidence = logsumexp(log_evidences, b=1. / len(self))
+ result.log_bayes_factor = result.log_evidence - result.log_noise_evidence
+
+ # Propogate uncertainty in combined evidence
log_errs = [res.log_evidence_err for res in self if np.isfinite(res.log_evidence_err)]
if len(log_errs) > 0:
result.log_evidence_err = 0.5 * logsumexp(2 * np.array(log_errs), b=1. / len(self))
else:
result.log_evidence_err = np.nan
+
+ # Combined posteriors with a weighting
+ result_weights = np.exp(log_evidences - np.max(log_evidences))
posteriors = list()
for res, frac in zip(self, result_weights):
selected_samples = (np.random.uniform(size=len(res.posterior)) < frac)
posteriors.append(res.posterior[selected_samples])
+
# remove original nested_samples
result.nested_samples = None
result.sampler_kwargs = None
--
GitLab