diff --git a/test/prior_tests.py b/test/prior_tests.py
index cc105c4d81618e6f6253d245fa09b309cba36e19..566dc3f2897a2a0adb4b3c89538eca8aedab2cb1 100644
--- a/test/prior_tests.py
+++ b/test/prior_tests.py
@@ -111,5 +111,75 @@ class TestFixMethod(unittest.TestCase):
self.assertRaises(ValueError, tupak.prior.fix, self.prior, np.nan)
+class TestPriorClasses(unittest.TestCase):
+
+ def setUp(self):
+
+ self.priors = [
+ tupak.prior.DeltaFunction(name='test', peak=1),
+ tupak.prior.Gaussian(name='test', mu=0, sigma=1),
+ tupak.prior.PowerLaw(name='test', alpha=0, minimum=0, maximum=1),
+ tupak.prior.PowerLaw(name='test', alpha=-1, minimum=1, maximum=1e2),
+ tupak.prior.Uniform(name='test', minimum=0, maximum=1),
+ tupak.prior.UniformComovingVolume(name='test', minimum=2e2, maximum=5e3),
+ tupak.prior.Sine(name='test'),
+ tupak.prior.Cosine(name='test'),
+ tupak.prior.Interped(name='test', xx=np.linspace(0, 10, 1000), yy=np.linspace(0, 10, 1000)**4,
+ minimum=3, maximum=5),
+ tupak.prior.TruncatedGaussian(name='test', mu=1, sigma=0.4, minimum=-1, maximum=1)
+ ]
+
+ def test_rescaling(self):
+ for prior in self.priors:
+ """Test the the rescaling works as expected."""
+ minimum_sample = prior.rescale(0)
+ self.assertAlmostEqual(minimum_sample, prior.minimum)
+ maximum_sample = prior.rescale(1)
+ self.assertAlmostEqual(maximum_sample, prior.maximum)
+ many_samples = prior.rescale(np.random.uniform(0, 1, 1000))
+ self.assertTrue(all((many_samples >= prior.minimum) & (many_samples <= prior.maximum)))
+ self.assertRaises(ValueError, lambda: prior.rescale(-1))
+
+ def test_sampling(self):
+ """Test that sampling from the prior always returns values within its domain."""
+ for prior in self.priors:
+ single_sample = prior.sample()
+ self.assertTrue((single_sample >= prior.minimum) & (single_sample <= prior.maximum))
+ many_samples = prior.sample(1000)
+ self.assertTrue(all((many_samples >= prior.minimum) & (many_samples <= prior.maximum)))
+
+ def test_prob(self):
+ """Test that the prior probability is non-negative in domain of validity and zero outside."""
+ for prior in self.priors:
+ # skip delta function prior in this case
+ if isinstance(prior, tupak.prior.DeltaFunction):
+ continue
+ if prior.maximum != np.inf:
+ outside_domain = np.linspace(prior.maximum + 1, prior.maximum + 1e4, 1000)
+ self.assertTrue(all(prior.prob(outside_domain) == 0))
+ if prior.minimum != -np.inf:
+ outside_domain = np.linspace(prior.minimum - 1e4, prior.minimum - 1, 1000)
+ self.assertTrue(all(prior.prob(outside_domain) == 0))
+ if prior.minimum == -np.inf:
+ prior.minimum = -1e5
+ if prior.maximum == np.inf:
+ prior.maximum = 1e5
+ domain = np.linspace(prior.minimum, prior.maximum, 1000)
+ self.assertTrue(all(prior.prob(domain) >= 0))
+ surround_domain = np.linspace(prior.minimum - 1, prior.maximum + 1, 1000)
+ prior.prob(surround_domain)
+
+ def test_normalized(self):
+ """Test that each of the priors are normalised, this needs care for delta function and Gaussian priors"""
+ for prior in self.priors:
+ if isinstance(prior, tupak.prior.DeltaFunction):
+ continue
+ elif isinstance(prior, tupak.prior.Gaussian):
+ domain = np.linspace(-1e2, 1e2, 1000)
+ else:
+ domain = np.linspace(prior.minimum, prior.maximum, 1000)
+ self.assertAlmostEqual(np.trapz(prior.prob(domain), domain), 1, 3)
+
+
if __name__ == '__main__':
unittest.main()
diff --git a/tupak/prior.py b/tupak/prior.py
index 8eb15c937131cfdb5d69ed1d47ebf6ac12097169..bdce5908b9217ee95778bcd15636e359bb56286f 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):
"""