Merge branch 'prior_tests' into 'master'

Prior tests

See merge request Monash/tupak!35

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()

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):
         """