diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index d1350be5..a65c2aa3 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -455,10 +455,15 @@ class kern(parameterised):
                 p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1._psi1 * 2., Z, mu, S, target[ps2])
             # linear X bias
             elif p1.name == 'bias' and p2.name == 'linear':
-                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target)
+                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2])  # [ps1])
+                psi1 = np.zeros((mu.shape[0], Z.shape[0]))
+                p2.psi1(Z, mu, S, psi1)
+                p1.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps1])
             elif p2.name == 'bias' and p1.name == 'linear':
-                p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target)
-                pass
+                p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1])
+                psi1 = np.zeros((mu.shape[0], Z.shape[0]))
+                p1.psi1(Z, mu, S, psi1)
+                p2.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps2])
             # rbf X linear
             elif p1.name == 'linear' and p2.name == 'rbf':
                 raise NotImplementedError  # TODO
diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py
index 6d2a3e48..78a8732a 100644
--- a/GPy/kern/linear.py
+++ b/GPy/kern/linear.py
@@ -114,7 +114,7 @@ class linear(kernpart):
 
     def psi1(self,Z,mu,S,target):
         """the variance, it does nothing"""
-        self.K(mu,Z,target)
+        self._psi1 = self.K(mu, Z, target)
 
     def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target):
         """the variance, it does nothing"""
diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py
index b11b4532..5396e175 100644
--- a/GPy/testing/bgplvm_tests.py
+++ b/GPy/testing/bgplvm_tests.py
@@ -60,7 +60,7 @@ class BGPLVMTests(unittest.TestCase):
 
     #@unittest.skip('psi2 cross terms are NotImplemented for this combination')
     def test_linear_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
+        N, M, Q, D = 30, 5, 4, 30
         X = np.random.rand(N, Q)
         k = GPy.kern.linear(Q) +  GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
         K = k.K(X)
diff --git a/GPy/testing/psi_stat_tests.py b/GPy/testing/psi_stat_tests.py
index 6aeea60c..1a14e088 100644
--- a/GPy/testing/psi_stat_tests.py
+++ b/GPy/testing/psi_stat_tests.py
@@ -105,6 +105,18 @@ if __name__ == "__main__":
     import sys
     interactive = 'i' in sys.argv
     if interactive:
+        N, M, Q, D = 30, 5, 4, 30
+        X = numpy.random.rand(N, Q)
+        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        K = k.K(X)
+        Y = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T
+        Y -= Y.mean(axis=0)
+        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
+        m.ensure_default_constraints()
+        m.randomize()
+#         self.assertTrue(m.checkgrad())
+
         Q = 5
         N = 50
         M = 10
@@ -119,17 +131,17 @@ if __name__ == "__main__":
                GPy.kern.linear(Q) + GPy.kern.bias(Q),
                GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
 
-        for k in kernels:
-            m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
-                     M=M, kernel=k)
-            assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
+#         for k in kernels:
+#             m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
+#                      M=M, kernel=k)
+#             assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
 #
 #         m0 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
 #                          M=M, kernel=GPy.kern.linear(Q))
 #         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
 #                          M=M, kernel=kernel)
-        m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
-                         M=M, kernel=kernel)
+#         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
+#                          M=M, kernel=kernel)
         m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
                          M=M, kernel=GPy.kern.rbf(Q))
         m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,