Skip to content
Snippets Groups Projects
Commit 4ff3dd42 authored by Christian Chapman-Bird's avatar Christian Chapman-Bird
Browse files

Added custom Bilby hyper-likelihood.

parent c9f62659
No related branches found
No related tags found
Loading
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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment