Source code for pycbc.population.fgmc_laguerre

# Copyright (C) 2016 Jolien Creighton
#           (C) 2021 Jolien Creighton & Thomas Dent
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.

Based ultimately on code used for O1 rate calculations, see
and technical documentation at

import numpy
import scipy.stats as sst
import scipy.special as ssp
import scipy.integrate as sig
import scipy.optimize as sop

[docs]class augmented_rv_continuous(sst.rv_continuous): def __init__(self, unit='dimensionless', texunit=r'\mbox{dimensionless}', texsymb=r'x', **kwargs): ''' Parameters ---------- unit : string, optional units of independent variable texunit : string, optional units of independent variable, in tex format texsymb : string, optional symbol of independent variable, in tex format ''' super(augmented_rv_continuous, self).__init__(**kwargs) self._hpd_interval_vec = numpy.vectorize(self._hpd_interval_scalar) self.unit = unit self.texunit = texunit self.texsymb = texsymb def _hpd_interval_scalar(self, alpha): def width(a): return self.ppf(alpha + self.cdf(a)) - a a = self.ppf(1e-6) # a is displaced slightly from 0 b = self.ppf(alpha) if self.pdf(a) >= self.pdf(b): # upper limit return self.a, b a = sop.fminbound(width, a, self.ppf(1.0 - alpha)) b = self.ppf(alpha + self.cdf(a)) return a, b
[docs] def hpd_interval(self, alpha): ''' Confidence interval of highest probability density. Parameters ---------- alpha : array_like of float Probability that an rv will be drawn from the returned range. Each value should be in the range [0, 1]. Returns ------- a, b : ndarray of float end-points of range that contain ``100 * alpha %`` of the rv's possible values. ''' if isinstance(alpha, (float, numpy.number)): a, b = self._hpd_interval_scalar(alpha) else: a, b = self._hpd_interval_vec(alpha) return a, b
[docs]class count_posterior(augmented_rv_continuous): ''' Count posterior distribution. ''' def __init__(self, logbf, laguerre_n, Lambda0, prior=-0.5, name='count posterior', unit='signals/experiment', texunit=r'\mathrm{signals}/\mathrm{experiment}', texsymb=r'\Lambda_1'): ''' Parameters ---------- logbf : array_like logs of normalized foreground over background pdf ratios of events laguerre_n: int degree of generalized Laguerre polynomial for quadrature formula Lambda0 : float background rate (default=len(bayesfac)) prior : float or count_posterior, optional prior distribution power law of improper prior if float or count posterior distribution if count_posterior (default=-0.5: Jeffreys prior) ''' super(count_posterior, self).__init__(a=0.0, b=numpy.inf, name=name, unit=unit, texunit=texunit, texsymb=texsymb) self.Lambda0 = Lambda0 # weighted Bayes factor self.k = numpy.exp(numpy.array(logbf)) / self.Lambda0 # power-law priors self.alpha = prior if prior == 0: self.prior = lambda x: 1.0 elif prior > 0: self.prior = lambda x: x ** prior else: # regularize at x = 0 self.prior = lambda x: (x + self.xtol) ** prior # pre-compute Gaussian-Generalized-Laguerre quadrature # abscissas and weights, along with pdf at these abscissas self.x, w = ssp.la_roots(laguerre_n, self.alpha) self.p = numpy.array([ww * + self.k * xx) for xx, ww in zip(self.x, w)]) self.norm = 1.0 / sum(self.p) self.p *= self.norm def _pdf(self, x): # discourage underflows by evaluating ln L and using ln(1+x) function logL = -x + numpy.sum(numpy.log1p(self.k * x)) P = numpy.exp(logL) * self.prior(x) return self.norm * P def _cdf(self, x): return sig.quad(self._pdf, 0.0, x)
[docs] def expect(self, func): ''' Calculate expected value of a function with respect to the distribution. The expected value of a function ``f(x)`` with respect to a distribution ``dist`` is defined as:: E[x] = Integral(f(x) * dist.pdf(x)) Parameters ---------- func : callable Function for which integral is calculated. Takes only one argument. Returns ------- expect : float The calculated expected value. ''' # FIXME: not as feature rich as the expect method this overrides return sum(pp * func(xx) for xx, pp in zip(self.x, self.p))
def _munp(self, n): return self.expect(lambda x: x**n)
[docs] def p_bg(self, logbf): ''' Calculate the false alarm probabilities of the events. Parameters ---------- logbf : array_like Logs of foreground over background probability ratios of events. ''' # get weighted bayes factor k = numpy.exp(numpy.asarray(logbf)) / self.Lambda0 P0 = + numpy.outer(k, self.x)), self.p) if isinstance(k, (float, int, numpy.number)): return P0.item() if isinstance(k, numpy.ndarray) and k.ndim == 0: return P0.item() # except in special cases above, return array of values return P0
__all__ = ['augmented_rv_continuous', 'count_posterior']