pdfs: make UniPdf multivariate; properly test multivariate case

pybayes/numpywrap.py

@@ -6,6 +6,6 @@
 
 # just import and flatten numpy types and functions
 
-from numpy import any as np_any, array, asarray, diag, dot, dot as dotvv, ndarray, zeros
+from numpy import any as np_any, array, asarray, diag, dot, dot as dotvv, ndarray, prod, zeros
 from numpy.linalg import cholesky, inv, slogdet
 from numpy.random import normal, uniform

pybayes/numpywrap.pyx

@@ -9,7 +9,7 @@
 # not callable from python etc.
 
 from numpy cimport import_array, int, npy_intp, NPY_DOUBLE, PyArray_EMPTY, PyArray_ISCARRAY_RO, PyArray_ISFARRAY_RO
-from numpy import any as np_any, array, asarray, diag, empty, zeros
+from numpy import any as np_any, array, asarray, diag, empty, prod, zeros
 from numpy.linalg import cholesky, slogdet
 from numpy.random import normal, uniform
 

pybayes/pdfs.pxd

@@ -25,7 +25,7 @@ cdef class Pdf(CPdf):
 
 
 cdef class UniPdf(Pdf):
-    cdef public double a, b
+    cdef public ndarray a, b
 
 
 cdef class GaussPdf(Pdf):

pybayes/pdfs.py

@@ -78,40 +78,43 @@ def sample(self):
 
 
 class UniPdf(Pdf):
-    """Simple uniform single-dimensional probability density function
+    """Simple uniform multivariate probability density function
 
     .. math: f(x|a, b) = \THETA {x-a} \THETA {b-x} {1} \over {b-a}  TODO
     """
 
     def __init__(self, a, b):
         """Initialise uniform distribution with left point a and right point b
 
-        a must be greater that b
+        a must be greater (in each dimension) than b
         """
-        if b <= a:
-            raise ValueError("b must be grater than a")
-        self.a = float(a)
-        self.b = float(b)
+        self.a = asarray(a)
+        self.b = asarray(b)
+        if a.ndim != 1 or b.ndim != 1:
+            raise ValueError("both a and b must be 1D numpy arrays (vectors)")
+        if a.shape[0] != b.shape[0]:
+            raise ValueError("a must have same shape as b")
+        if np_any(self.b <= self.a):
+            raise ValueError("b must be greater than a in each dimension")
 
     def shape(self):
-        return 1
+        return self.a.shape[0]
 
     def mean(self):
-        return array([(self.a+self.b)/2.])
+        return (self.a+self.b)/2.  # element-wise division
 
     def variance(self):
-        return array([((self.b-self.a)**2)/12.])
+        return ((self.b-self.a)**2)/12.  # element-wise power and division
 
     def eval_log(self, x):
         if x is None:  # cython-specific, but wont hurt in python
             raise ValueError("x must be numpy.ndarray")
-        x0 = x[0]
-        if x0 <= self.a or x0 >= self.b:
+        if np_any(x <= self.a) or np_any(x >= self.b):
             return float('-inf')
-        return -log(self.b-self.a)
+        return -log(prod(self.b-self.a))
 
     def sample(self):
-        return uniform(self.a, self.b, self.shape())
+        return uniform(-0.5, 0.5, self.shape()) * (self.b-self.a) + self.mean()
 
 
 class GaussPdf(Pdf):

pybayes/tests/test_pdfs.py

@@ -58,24 +58,30 @@ class TestUniPdf(ut.TestCase):
     """Test uniform pdf"""
 
     def setUp(self):
-        self.uni = pb.UniPdf(-10., 20.)
+        self.uni = pb.UniPdf(np.array([-10.]), np.array([20.]))
+        (self.a, self.b) = (np.array([0., -1., 2.]), np.array([1., 1., 4.]))
+        self.multiuni = pb.UniPdf(self.a, self.b)
 
     def test_init(self):
         self.assertEqual(type(self.uni), pb.UniPdf)
+        self.assertEqual(type(self.multiuni), pb.UniPdf)
 
     def test_invalid_init(self):
-        self.assertRaises(ValueError, pb.UniPdf, 1.0, 0.5)
+        self.assertRaises(ValueError, pb.UniPdf, np.zeros(5), np.array([1., 2., 3., -0.01, 3.]))  # b must be > a element-wise
 
     def test_shape(self):
         self.assertEqual(self.uni.shape(), 1)
+        self.assertEqual(self.multiuni.shape(), 3)
 
     def test_mean(self):
         self.assertTrue(np.all(self.uni.mean() == np.array([5.])))
         self.assertTrue(np.all(self.uni.cmean(None) == np.array([5.])))  # test also cond. variant
+        self.assertTrue(np.all(self.multiuni.mean() == np.array([0.5, 0., 3.])))
 
     def test_variance(self):
-        self.assertTrue(np.all(self.uni.variance() == np.array(75.)))
-        self.assertTrue(np.all(self.uni.cvariance(None) == np.array(75.)))  # cond variant
+        self.assertTrue(np.all(self.uni.variance() == np.array([75.])))
+        self.assertTrue(np.all(self.uni.cvariance(None) == np.array([75.])))  # cond variant
+        self.assertTrue(np.all(self.multiuni.variance() == np.array([1./12., 1./3., 1./3.])))
 
     def test_eval_log(self):
         self.assertEqual(self.uni.eval_log(np.array([-10.1])), float('-inf'))
@@ -85,13 +91,30 @@ def test_eval_log(self):
         self.assertEqual(self.uni.eval_log(np.array([-10.1])), float('-inf'))
         self.assertEqual(self.uni.ceval_log(np.array([-10.1]), None), float('-inf'))
 
+        (x, y, z) = (0.2, 0.8, 3.141592)  # this point belongs to multiuni's interval
+        one_under = np.array([[-0.1, y, z],  # one of the values is lower
+                             [x, -1.04, z],
+                             [x, y, -3.975]])
+        one_over = np.array([[1465.67, y, z],  # one of the values is greater
+                             [x, 1.000456, z],
+                             [x, y, 4.67]])
+        self.assertEqual(self.multiuni.eval_log(np.array([x, y, z])), log(1./(1.*2*2)))
+        for i in range(3):
+            self.assertEqual(self.multiuni.eval_log(one_under[i]), float('-inf'))
+            self.assertEqual(self.multiuni.eval_log(one_over[i]), float('-inf'))
+
     def test_sample(self):
         for i in range(0, 100):
             sample = self.uni.sample()
             csample = self.uni.csample(None)
             self.assertTrue(-10. <= sample[0] <= 20.)  # test sample is within bounds
             self.assertTrue(-10. <= csample[0] <= 20.)  # also for conditional variant
 
+        for i in range(0, 500):
+            sample = self.multiuni.sample()
+            self.assertTrue(np.all(self.a <= sample))
+            self.assertTrue(np.all(sample <= self.b))
+
 
 class TestGaussPdf(ut.TestCase):
     """Test Gaussian pdf"""