From cf8cc886c2cb5de8c9cfc3fceeae9a8ddad8f42d Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Mon, 14 May 2018 12:09:25 +1000
Subject: [PATCH] update prob methods to allow more generic usage
---
tupak/prior.py | 58 ++++++++++++++++++++------------------------------
1 file changed, 23 insertions(+), 35 deletions(-)
diff --git a/tupak/prior.py b/tupak/prior.py
index 8eb15c937..bdce5908b 100644
--- a/tupak/prior.py
+++ b/tupak/prior.py
@@ -132,10 +132,8 @@ class Uniform(Prior):
def prob(self, val):
"""Return the prior probability of val"""
- if (self.minimum < val) & (val < self.maximum):
- return 1 / self.support
- else:
- return 0
+ in_prior = (val >= self.minimum) & (val <= self.maximum)
+ return 1 / self.support * in_prior
class DeltaFunction(Prior):
@@ -181,14 +179,12 @@ class PowerLaw(Prior):
def prob(self, val):
"""Return the prior probability of val"""
- if (val > self.minimum) and (val < self.maximum):
- if self.alpha == -1:
- return 1 / val / np.log(self.maximum / self.minimum)
- else:
- return val ** self.alpha * (1 + self.alpha) / (self.maximum ** (1 + self.alpha)
- - self.minimum ** (1 + self.alpha))
+ in_prior = (val >= self.minimum) & (val <= self.maximum)
+ if self.alpha == -1:
+ return np.nan_to_num(1 / val / np.log(self.maximum / self.minimum)) * in_prior
else:
- return 0
+ return np.nan_to_num(val ** self.alpha * (1 + self.alpha) / (self.maximum ** (1 + self.alpha)
+ - self.minimum ** (1 + self.alpha))) * in_prior
class Cosine(Prior):
@@ -205,13 +201,10 @@ class Cosine(Prior):
Prior.test_valid_for_rescaling(val)
return np.arcsin(-1 + val * 2)
- @staticmethod
- def prob(val):
+ def prob(self, val):
"""Return the prior probability of val, defined over [-pi/2, pi/2]"""
- if (val > np.minimum) and (val < np.maximum):
- return np.cos(val) / 2
- else:
- return 0
+ in_prior = (val >= self.minimum) & (val <= self.maximum)
+ return np.cos(val) / 2 * in_prior
class Sine(Prior):
@@ -226,15 +219,12 @@ class Sine(Prior):
This maps to the inverse CDF. This has been analytically solved for this case.
"""
Prior.test_valid_for_rescaling(val)
- return np.arccos(-1 + val * 2)
+ return np.arccos(1 - val * 2)
- @staticmethod
- def prob(val):
+ def prob(self, val):
"""Return the prior probability of val, defined over [0, pi]"""
- if (val > 0) and (val < np.pi):
- return np.sin(val) / 2
- else:
- return 0
+ in_prior = (val >= self.minimum) & (val <= self.maximum)
+ return np.sin(val) / 2 * in_prior
class Gaussian(Prior):
@@ -290,11 +280,9 @@ class TruncatedGaussian(Prior):
def prob(self, val):
"""Return the prior probability of val"""
- if (val > self.minimum) & (val < self.maximum):
- return np.exp(-(self.mu - val) ** 2 / (2 * self.sigma ** 2)) / (
- 2 * np.pi) ** 0.5 / self.sigma / self.normalisation
- else:
- return 0
+ in_prior = (val >= self.minimum) & (val <= self.maximum)
+ return np.exp(-(self.mu - val) ** 2 / (2 * self.sigma ** 2)) / (
+ 2 * np.pi) ** 0.5 / self.sigma / self.normalisation * in_prior
class Interped(Prior):
@@ -302,6 +290,7 @@ class Interped(Prior):
def __init__(self, xx, yy, minimum=None, maximum=None, name=None, latex_label=None):
"""Initialise object from arrays of x and y=p(x)"""
Prior.__init__(self, name, latex_label)
+ all_interpolated = interp1d(x=xx, y=yy, bounds_error=False, fill_value=0)
if minimum is None or minimum < min(xx):
self.minimum = min(xx)
else:
@@ -310,22 +299,21 @@ class Interped(Prior):
self.maximum = max(xx)
else:
self.maximum = maximum
- self.xx = xx[(xx > self.minimum) & (xx < self.maximum)]
- self.yy = yy[(xx > self.minimum) & (xx < self.maximum)]
+ self.xx = np.linspace(self.minimum, self.maximum, len(xx))
+ self.yy = all_interpolated(self.xx)
if np.trapz(self.yy, self.xx) != 1:
logging.info('Supplied PDF is not normalised, normalising.')
self.yy /= np.trapz(self.yy, self.xx)
self.YY = cumtrapz(self.yy, self.xx, initial=0)
+ # Need last element of cumulative distribution to be exactly one.
+ self.YY[-1] = 1
self.probability_density = interp1d(x=self.xx, y=self.yy, bounds_error=False, fill_value=0)
self.cumulative_distribution = interp1d(x=self.xx, y=self.YY, bounds_error=False, fill_value=0)
self.inverse_cumulative_distribution = interp1d(x=self.YY, y=self.xx, bounds_error=True)
def prob(self, val):
"""Return the prior probability of val"""
- if (val > self.minimum) & (val < self.maximum):
- return self.probability_density(val)
- else:
- return 0
+ return self.probability_density(val)
def rescale(self, val):
"""
--
GitLab