diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py
index 37c424ad..6893011d 100644
--- a/GPy/kern/_src/add.py
+++ b/GPy/kern/_src/add.py
@@ -23,7 +23,7 @@ class Add(CombinationKernel):
                     subkerns.insert(i, part)
         super(Add, self).__init__(subkerns, name)
         self._exact_psicomp = self._check_exact_psicomp()
-        
+
     def _check_exact_psicomp(self):
         from .. import RBF,Linear,Bias,White
         n_kerns = len(self.parts)
@@ -164,11 +164,11 @@ class Add(CombinationKernel):
             # rbf X bias
             #elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)):
             elif isinstance(p1,  Bias) and isinstance(p2, (RBF, Linear)):
-                tmp = p2.psi1(Z, variational_posterior).sum(axis=0)
+                tmp = p2.psi1(Z, variational_posterior)
                 psi2 += p1.variance * (tmp[:, :, None] + tmp[:, None, :])
             #elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)):
             elif isinstance(p2, Bias) and isinstance(p1, (RBF, Linear)):
-                tmp = p1.psi1(Z, variational_posterior).sum(axis=0)
+                tmp = p1.psi1(Z, variational_posterior)
                 psi2 += p2.variance * (tmp[:, :, None] + tmp[:, None, :])
             elif isinstance(p2, (RBF, Linear)) and isinstance(p1, (RBF, Linear)):
                 assert np.intersect1d(p1.active_dims, p2.active_dims).size == 0, "only non overlapping kernel dimensions allowed so far"
diff --git a/GPy/kern/_src/static.py b/GPy/kern/_src/static.py
index c7f7b118..bcfec4a7 100644
--- a/GPy/kern/_src/static.py
+++ b/GPy/kern/_src/static.py
@@ -103,12 +103,17 @@ class Bias(Static):
         return ret
 
     def psi2n(self, Z, variational_posterior):
-        ret = np.empty((1, Z.shape[0], Z.shape[0]), dtype=np.float64)
+        ret = np.empty((variational_posterior.mean.shape[0], Z.shape[0], Z.shape[0]), dtype=np.float64)
         ret[:] = self.variance*self.variance
         return ret
 
     def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-        self.variance.gradient = dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum()*variational_posterior.shape[0]
+        if dL_dpsi2.ndim == 2:
+            self.variance.gradient = (dL_dpsi0.sum() + dL_dpsi1.sum()
+                                    + 2.*self.variance*dL_dpsi2.sum()*variational_posterior.shape[0])
+        else:
+            self.variance.gradient = (dL_dpsi0.sum() + dL_dpsi1.sum()
+                                    + 2.*self.variance*dL_dpsi2.sum())
 
 class Fixed(Static):
     def __init__(self, input_dim, covariance_matrix, variance=1., active_dims=None, name='fixed'):
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 5dcead8a..1fd92d98 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -503,6 +503,7 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase):
 
         from GPy.models import GradientChecker
         m = GradientChecker(f, df, kernel.param_array.copy())
+        m.checkgrad(verbose=1)
         self.assertTrue(m.checkgrad())
 
     def _test_Z(self, kernel, psi2n=False):