diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 619d1687..d686064a 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -462,10 +462,8 @@ class kern(Parameterized):
                 pass
             # rbf X bias
             elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)):
-                target += 2 * p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :])
+                target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :])
             elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)):
-                tmp1 = p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :])
-                renorm = p1.variance*np.exp()
                 target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :])
             # linear X bias
             elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, Linear):
@@ -478,12 +476,21 @@ class kern(Parameterized):
                 target += p2.variance * (tmp[:, :, None] + tmp[:, None, :])
             # rbf X any
             elif isinstance(p1, (RBF, RBFInv)):
-                pass
+                psi11 = np.zeros((mu.shape[0], Z.shape[0]))
+                psi12 = np.zeros((mu.shape[0], Z.shape[0]))
+                p1.psi1(Z, mu, S, psi11)
+                p2.psi1(Z, mu, S, psi12)
+                
+                crossterms  = psi11[:, :, None] + psi12[:, None, :]
+                crossterms += psi12[:, :, None] + psi11[:, None, :]
+                
+                target += p1._crossterm_product_expectation(p2, Z, mu, S)
+                #import ipdb;ipdb.set_trace()
             elif isinstance(p2, (RBF, RBFInv)):
                 raise NotImplementedError # TODO
             else:
                 raise NotImplementedError, "psi2 cannot be computed for this kernel"
-        return target
+        return target        
 
     def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
         target = np.zeros(self.num_params)
diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py
index 585d687f..56a6b0eb 100644
--- a/GPy/kern/parts/rbf.py
+++ b/GPy/kern/parts/rbf.py
@@ -208,6 +208,16 @@ class RBF(Kernpart):
         self._psi_computations(Z, mu, S)
         target += self._psi2
 
+    def _crossterm_product_expectation(self, K, Z, mu, S):
+        # compute the crossterm expectation for K as the other kernel:
+        import ipdb;ipdb.set_trace()
+        Sigma = 1./self.lengthscale[None,:] + 1./S # is independent across M, 
+        M = (Z[None,:,:]/self.lengthscale[None,None,:] + (mu/S)[:,None,:]) / Sigma[:,None,:]
+        psi1_other = K.psi1()
+        self.variance
+        # return is [N x M x M]
+        return 
+
     def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
         """Shape N,num_inducing,num_inducing,Ntheta"""
         self._psi_computations(Z, mu, S)
diff --git a/GPy/testing/psi_stat_expectation_tests.py b/GPy/testing/psi_stat_expectation_tests.py
index 16904927..ae3d1022 100644
--- a/GPy/testing/psi_stat_expectation_tests.py
+++ b/GPy/testing/psi_stat_expectation_tests.py
@@ -27,7 +27,7 @@ def ard(p):
 @testing.deepTest(__test__())
 class Test(unittest.TestCase):
     input_dim = 9
-    num_inducing = 4
+    num_inducing = 13
     N = 30
     Nsamples = 9e6
 
@@ -51,13 +51,16 @@ class Test(unittest.TestCase):
 #                        GPy.kern.bias(self.input_dim) +
 #                        GPy.kern.white(self.input_dim)),
 #                     (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) +
-#                     GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) +
-#                     GPy.kern.linear(self.input_dim, np.random.rand(self.input_dim), ARD=True) +
+                    (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
+                     +GPy.kern.linear(self.input_dim, np.random.rand(self.input_dim), ARD=True)
 #                     GPy.kern.bias(self.input_dim) +
 #                     GPy.kern.white(self.input_dim)),
-        (GPy.kern.linear(self.input_dim, np.random.rand(self.input_dim), ARD=True) +
-                    GPy.kern.bias(self.input_dim, np.random.rand()) +
-                    GPy.kern.white(self.input_dim, np.random.rand())),
+                    ),
+        (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
+         +GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
+         #+GPy.kern.bias(self.input_dim, np.random.rand())
+         #+GPy.kern.white(self.input_dim, np.random.rand())),
+         ),
                 (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) +
                     GPy.kern.bias(self.input_dim, np.random.rand()) +
                     GPy.kern.white(self.input_dim, np.random.rand())),