likelihood.py 3.62 KB
Newer Older
1 2 3 4 5 6 7 8 9
from __future__ import division, print_function

import logging
import numpy as np
from ..core.likelihood import Likelihood
from .model import Model


class HyperparameterLikelihood(Likelihood):
10
    """ A likelihood for inferring hyperparameter posterior distributions
11

12
    See Eq. (34) of https://arxiv.org/abs/1809.02293 for a definition.
13 14 15 16

    Parameters
    ----------
    posteriors: list
17 18
        An list of pandas data frames of samples sets of samples.
        Each set may have a different size.
Colm Talbot's avatar
Colm Talbot committed
19
    hyper_prior: `bilby.hyper.model.Model`
20
        The population model, this can alternatively be a function.
Colm Talbot's avatar
Colm Talbot committed
21
    sampling_prior: `bilby.hyper.model.Model`
22
        The sampling prior, this can alternatively be a function.
23 24 25 26 27
    log_evidences: list, optional
        Log evidences for single runs to ensure proper normalisation
        of the hyperparameter likelihood. If not provided, the original
        evidences will be set to 0. This produces a Bayes factor between
        the sampling prior and the hyperparameterised model.
28 29 30 31 32
    max_samples: int, optional
        Maximum number of samples to use from each set.

    """

33 34
    def __init__(self, posteriors, hyper_prior, sampling_prior,
                 log_evidences=None, max_samples=1e100):
35 36 37 38
        if not isinstance(hyper_prior, Model):
            hyper_prior = Model([hyper_prior])
        if not isinstance(sampling_prior, Model):
            sampling_prior = Model([sampling_prior])
39 40 41 42
        if log_evidences is not None:
            self.evidence_factor = np.sum(log_evidences)
        else:
            self.evidence_factor = np.nan
43 44 45 46 47 48 49 50 51
        self.posteriors = posteriors
        self.hyper_prior = hyper_prior
        self.sampling_prior = sampling_prior
        self.max_samples = max_samples
        Likelihood.__init__(self, hyper_prior.parameters)

        self.data = self.resample_posteriors()
        self.n_posteriors = len(self.posteriors)
        self.samples_per_posterior = self.max_samples
52 53
        self.samples_factor =\
            - self.n_posteriors * np.log(self.samples_per_posterior)
54

55
    def log_likelihood_ratio(self):
56
        self.hyper_prior.parameters.update(self.parameters)
Gregory Ashton's avatar
Gregory Ashton committed
57
        log_l = np.sum(np.log(np.sum(self.hyper_prior.prob(self.data) /
58 59
                       self.sampling_prior.prob(self.data), axis=-1)))
        log_l += self.samples_factor
60 61
        return np.nan_to_num(log_l)

62 63 64 65 66 67
    def noise_log_likelihood(self):
        return self.evidence_factor

    def log_likelihood(self):
        return self.noise_log_likelihood() + self.log_likelihood_ratio()

68
    def resample_posteriors(self, max_samples=None):
69 70 71 72 73 74 75 76 77 78 79 80 81 82
        """
        Convert list of pandas DataFrame object to dict of arrays.

        Parameters
        ----------
        max_samples: int, opt
            Maximum number of samples to take from each posterior,
            default is length of shortest posterior chain.
        Returns
        -------
        data: dict
            Dictionary containing arrays of size (n_posteriors, max_samples)
            There is a key for each shared key in self.posteriors.
        """
83 84 85 86 87
        if max_samples is not None:
            self.max_samples = max_samples
        for posterior in self.posteriors:
            self.max_samples = min(len(posterior), self.max_samples)
        data = {key: [] for key in self.posteriors[0]}
88 89
        logging.debug('Downsampling to {} samples per posterior.'.format(
            self.max_samples))
90 91 92 93 94 95 96
        for posterior in self.posteriors:
            temp = posterior.sample(self.max_samples)
            for key in data:
                data[key].append(temp[key])
        for key in data:
            data[key] = np.array(data[key])
        return data