tests pass, small bugfixes

Author: bwengals

pymc3/gp/cov.py

@@ -4,14 +4,17 @@
 from functools import reduce
 from operator import mul, add
 
-__all__ = ['ExpQuad',
+__all__ = ['Constant',
+           'WhiteNoise',
+           'ExpQuad',
            'RatQuad',
            'Exponential',
            'Matern52',
            'Matern32',
            'Linear',
            'Polynomial',
            'Cosine',
+           'Periodic',
            'WarpedInput',
            'Gibbs']
 
@@ -55,6 +58,12 @@ def __call__(self, X, Xs=None, diag=False):
         else:
             return self.full(X, Xs)
 
+    def diag(self, X):
+        raise NotImplementedError
+
+    def full(self, X, Xs):
+        raise NotImplementedError
+
     def _slice(self, X, Xs):
         X = tt.as_tensor_variable(X[:, self.active_dims])
         if Xs is not None:
@@ -172,13 +181,13 @@ def __init__(self, sigma):
         self.sigma = sigma
 
     def diag(self, X):
-        return tt.square(self.sigma) * tt.ones(tt.stack([X.shape[0], ]))
+        return tt.alloc(tt.square(self.sigma), X.shape[0])
 
     def full(self, X, Xs=None):
         if Xs is None:
             return tt.square(self.sigma) * tt.eye(X.shape[0])
         else:
-            return tt.zeros(tt.stack([X.shape[0], Xs.shape[0]]))
+            return tt.alloc(0.0, X.shape[0], Xs.shape[0])
 
 
 class Stationary(Covariance):
@@ -187,9 +196,9 @@ class Stationary(Covariance):
 
     Parameters
     ----------
-    ls : If input_dim > 1, a list or array of scalars or PyMC3 random
+    ls : Lengthscale.  If input_dim > 1, a list or array of scalars or PyMC3 random
     variables.  If input_dim == 1, a scalar or PyMC3 random variable.
-    ls_inv : 1 / ls.  One of ls or ls_inv must be provided.
+    ls_inv : Inverse lengthscale.  1 / ls.  One of ls or ls_inv must be provided.
     """
 
     def __init__(self, input_dim, ls=None, ls_inv=None, active_dims=None):
@@ -222,13 +231,19 @@ def euclidean_dist(self, X, Xs):
 
     def diag(self, X):
         return tt.alloc(1.0, X.shape[0])
-        #return tt.ones(tt.stack([X.shape[0], ]))
 
     def full(self, X, Xs=None):
         raise NotImplementedError
 
 
 class Periodic(Stationary):
+    R"""
+    The Periodic kernel.
+
+    .. math::
+       k(x, x') = \mathrm{exp}\left( -\frac{2 \mathrm{sin}^2(\pi |x-x'| \frac{1}{T})}{\ell^2} \right)
+    """
+
     def __init__(self, input_dim, period, ls=None, ls_inv=None, active_dims=None):
         super(Periodic, self).__init__(input_dim, ls, ls_inv, active_dims)
         self.period = period
@@ -490,8 +505,7 @@ def full(self, X, Xs=None):
                 * tt.exp(-1.0 * r2 / (rx2 + rz2)))
 
     def diag(self, X):
-        return tt.ones(tt.stack([X.shape[0], ]))
-
+        return tt.alloc(1.0, X.shape[0])
 
 def handle_args(func, args):
     def f(x, args):

pymc3/gp/gp.py

@@ -24,7 +24,7 @@ def __init__(self, mean_func=Zero(), cov_func=Constant(0.0)):
 
     def __add__(self, other):
         same_attrs = set(self.__dict__.keys()) == set(other.__dict__.keys())
-        if not isinstance(self, type(other)) and not same_attrs:
+        if not isinstance(self, type(other)) or not same_attrs:
             raise ValueError("cant add different GP types")
         mean_total = self.mean_func + other.mean_func
         cov_total = self.cov_func + other.cov_func

pymc3/tests/test_gp.py

@@ -312,6 +312,52 @@ def test_2dard(self):
         Kd = theano.function([], cov(X, diag=True))()
         npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
 
+    def test_inv_lengthscale(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        with pm.Model() as model:
+            cov = pm.gp.cov.ExpQuad(1, ls_inv=10)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 0.53940, atol=1e-3)
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 0.53940, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+
+
+class TestWhiteNoise(object):
+    def test_1d(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        with pm.Model() as model:
+            cov = pm.gp.cov.WhiteNoise(sigma=0.5)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 0.0, atol=1e-3)
+        npt.assert_allclose(K[0, 0], 0.5**2, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+        # check predict
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 0.0, atol=1e-3)
+        # white noise predicting should return all zeros
+        npt.assert_allclose(K[0, 0], 0.0, atol=1e-3)
+
+
+class TestConstant(object):
+    def test_1d(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        with pm.Model() as model:
+            cov = pm.gp.cov.Constant(2.5)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 2.5, atol=1e-3)
+        npt.assert_allclose(K[0, 0], 2.5, atol=1e-3)
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 2.5, atol=1e-3)
+        npt.assert_allclose(K[0, 0], 2.5, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+
 
 class TestRatQuad(object):
     def test_1d(self):
@@ -383,6 +429,20 @@ def test_1d(self):
         npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
 
 
+class TestPeriodic(object):
+    def test_1d(self):
+        X = np.linspace(0, 1, 10)[:, None]
+        with pm.Model() as model:
+            cov = pm.gp.cov.Periodic(1, 0.1, 0.1)
+        K = theano.function([], cov(X))()
+        npt.assert_allclose(K[0, 1], 0.00288, atol=1e-3)
+        K = theano.function([], cov(X, X))()
+        npt.assert_allclose(K[0, 1], 0.00288, atol=1e-3)
+        # check diagonal
+        Kd = theano.function([], cov(X, diag=True))()
+        npt.assert_allclose(np.diag(K), Kd, atol=1e-5)
+
+
 class TestLinear(object):
     def test_1d(self):
         X = np.linspace(0, 1, 10)[:, None]
@@ -491,7 +551,7 @@ def setup_method(self):
             cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3])
             mean_func = pm.gp.mean.Constant(0.5)
             gp = pm.gp.Marginal(mean_func, cov_func)
-            f = gp.marginal_likelihood("f", X, y, noise=0.0)
+            f = gp.marginal_likelihood("f", X, y, noise=0.0, is_observed=False, observed=y)
             p = gp.conditional("p", Xnew)
         self.logp = model.logp({"p": pnew})
         self.X = X
@@ -558,27 +618,28 @@ def testApproximations(self, approx):
         approx_logp = model.logp({"f": self.y, "p": self.pnew})
         npt.assert_allclose(approx_logp, self.logp, atol=0, rtol=1e-2)
 
-    def testPredictCov(self):
+    @pytest.mark.parametrize('approx', ['FITC', 'VFE', 'DTC'])
+    def testPredictVar(self, approx):
         with pm.Model() as model:
             cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3])
             mean_func = pm.gp.mean.Constant(0.5)
-            gp = pm.gp.MarginalSparse(mean_func, cov_func, approx="DTC")
+            gp = pm.gp.MarginalSparse(mean_func, cov_func, approx=approx)
             f = gp.marginal_likelihood("f", self.X, self.X, self.y, self.sigma)
-        mu1, cov1 = self.gp.predict(self.Xnew, pred_noise=True)
-        mu2, cov2 = gp.predict(self.Xnew, pred_noise=True)
+        mu1, var1 = self.gp.predict(self.Xnew, diag=True)
+        mu2, var2 = gp.predict(self.Xnew, diag=True)
         npt.assert_allclose(mu1, mu2, atol=0, rtol=1e-3)
-        npt.assert_allclose(cov1, cov2, atol=0, rtol=1e-3)
+        npt.assert_allclose(var1, var2, atol=0, rtol=1e-3)
 
-    def testPredictVar(self):
+    def testPredictCov(self):
         with pm.Model() as model:
             cov_func = pm.gp.cov.ExpQuad(3, [0.1, 0.2, 0.3])
             mean_func = pm.gp.mean.Constant(0.5)
             gp = pm.gp.MarginalSparse(mean_func, cov_func, approx="DTC")
-            f = gp.marginal_likelihood("f", self.X, self.X, self.y, self.sigma)
-        mu1, var1 = self.gp.predict(self.Xnew, diag=True)
-        mu2, var2 = gp.predict(self.Xnew, diag=True)
+            f = gp.marginal_likelihood("f", self.X, self.X, self.y, self.sigma, is_observed=False)
+        mu1, cov1 = self.gp.predict(self.Xnew, pred_noise=True)
+        mu2, cov2 = gp.predict(self.Xnew, pred_noise=True)
         npt.assert_allclose(mu1, mu2, atol=0, rtol=1e-3)
-        npt.assert_allclose(var1, var2, atol=0, rtol=1e-3)
+        npt.assert_allclose(cov1, cov2, atol=0, rtol=1e-3)
 
 
 class TestGPAdditive(object):
@@ -621,7 +682,7 @@ def testAdditiveMarginal(self):
 
     @pytest.mark.parametrize('approx', ['FITC', 'VFE', 'DTC'])
     def testAdditiveMarginalSparse(self, approx):
-        Xu = np.random.randn(10, 1)
+        Xu = np.random.randn(10, 3)
         sigma = 0.1
         with pm.Model() as model1:
             gp1 = pm.gp.MarginalSparse(self.means[0], self.covs[0], approx=approx)
@@ -669,8 +730,8 @@ def testAdditiveLatent(self):
             fp2 = gptot.conditional("fp2", self.Xnew)
 
         fp = np.random.randn(self.Xnew.shape[0])
-        npt.assert_allclose(fp1.logp({"fp1": fp}), fp2.logp({"fp2": fp}), atol=0, rtol=1e-2)
-
+        npt.assert_allclose(fp1.logp({"fsum": self.y, "fp1": fp}),
+                            fp2.logp({"fsum": self.y, "fp2": fp}), atol=0, rtol=1e-2)
 
     def testAdditiveSparseRaises(self):
         # cant add different approximations