From 4ff3dd42543dacb8adbd8d3ccfbba3d639afcc97 Mon Sep 17 00:00:00 2001
From: Christian Chapman-Bird <c.chapmanbird@btinternet.com>
Date: Thu, 4 Nov 2021 16:04:08 +0000
Subject: [PATCH] Added custom Bilby hyper-likelihood.

---
 src/inference/likelihood.py | 50 +++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)
 create mode 100644 src/inference/likelihood.py

diff --git a/src/inference/likelihood.py b/src/inference/likelihood.py
new file mode 100644
index 0000000..0709ced
--- /dev/null
+++ b/src/inference/likelihood.py
@@ -0,0 +1,50 @@
+import bilby
+import numpy as np
+import os
+import pandas as pd
+
+class HyperLikelihood(bilby.Likelihood):
+    """
+    Hyperparameter Likelihood class. Calculates loglikelihood using samples from pandas dataframes in directory
+    eventdir (one per event).
+
+    If the population content is changed, this must be re-initialised.
+
+    pdet (optional): Detection probability function for the population parameters. If supplied, detection selection
+    effects will be included.
+    """
+
+    def __init__(self, population_container, eventdir, pdet=None):
+        self.eventdir = eventdir
+        self.populations = population_container
+
+        self.n_obs = len(os.listdir(eventdir))
+        self.events = [pd.read_csv(ev) for ev in os.listdir(eventdir) if '.csv' in ev]
+
+        self.pdetfun = pdet
+
+        super().__init__()
+
+    def _update_params(self):
+        self.populations.set_hypers(self.parameters)
+
+    def _get_beta(self):
+        if self.pdetfun is not None:
+            # draw samples from population, plug in to p_det(theta).
+            samples = self.populations.sample(int(1e5))
+            pdetvals = self.pdetfun(samples)
+            return np.mean(pdetvals)  # take the mean over our population
+        else:
+            return 1  # without selection effects, 100% of the population is assumed to be detectable
+
+    def log_likelihood(self):
+        self._update_params()
+        bet = np.log(self._get_beta())
+        log_lhood = bet * self.n_obs  # N_obs * log(beta(lambda))
+
+        for event in self.events:
+            sample_ppops = self.populations.get_pdf(dict(zip(event.index, event.values)))  # dictionary of {key: samples} pairs
+            sum_over_samples = np.sum(sample_ppops)  # take the sum of our samples
+            log_lhood += np.log(sum_over_samples) # take the log of the sum and add it to the log likelihood
+            log_lhood -= np.log(len(event)) # 1/S_i # normalise by subtracting the log of the number of samples
+        return log_lhood
\ No newline at end of file
-- 
GitLab