diff --git a/GPy/examples/BGPLVM_demo.py b/GPy/examples/BGPLVM_demo.py
index 0fd58571..e92856ab 100644
--- a/GPy/examples/BGPLVM_demo.py
+++ b/GPy/examples/BGPLVM_demo.py
@@ -8,21 +8,19 @@ np.random.seed(123344)
 
 N = 10
 M = 3
-Q = 4
-D = 5
+Q = 2
+D = 4
 #generate GPLVM-like data
 X = np.random.rand(N, Q)
 k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
 K = k.K(X)
 Y = np.random.multivariate_normal(np.zeros(N),K,D).T
 
-<<<<<<< HEAD
-k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
-# k = GPy.kern.linear(Q, ARD = True)  + GPy.kern.white(Q, 0.00001)
-=======
+k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q)
+# k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
 # k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-k = GPy.kern.rbf(Q, ARD = False)  + GPy.kern.white(Q, 0.00001)
->>>>>>> master
+# k = GPy.kern.rbf(Q, ARD = False)  + GPy.kern.white(Q, 0.00001)
+
 m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
 m.constrain_positive('(rbf|bias|noise|white|S)')
 # m.constrain_fixed('S', 1)
diff --git a/GPy/examples/oil_flow_demo.py b/GPy/examples/oil_flow_demo.py
index 91fff4e4..71fb1bd3 100644
--- a/GPy/examples/oil_flow_demo.py
+++ b/GPy/examples/oil_flow_demo.py
@@ -32,7 +32,7 @@ Y -= Y.mean(axis=0)
 # Y /= Y.std(axis=0)
 
 Q = 5
-k = GPy.kern.linear(Q, ARD = False) + GPy.kern.white(Q)
+k = GPy.kern.linear(Q, ARD = True) + GPy.kern.white(Q)
 m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k, M = 20)
 m.constrain_positive('(rbf|bias|S|linear|white|noise)')
 
@@ -43,7 +43,7 @@ m.constrain_positive('(rbf|bias|S|linear|white|noise)')
 # plot_oil(m.X, np.array([1,1]), labels, 'PCA initialization')
 m.optimize(messages = True)
 # m.optimize('tnc', messages = True)
-plot_oil(m.X, m.kern.parts[0].lengthscale, labels, 'B-GPLVM')
+# plot_oil(m.X, m.kern.parts[0].lengthscale, labels, 'B-GPLVM')
 # # pb.figure()
 # m.plot()
 # pb.title('PCA initialisation')
diff --git a/GPy/kern/bias.py b/GPy/kern/bias.py
index d97f2f00..91594e4c 100644
--- a/GPy/kern/bias.py
+++ b/GPy/kern/bias.py
@@ -47,6 +47,10 @@ class bias(kernpart):
     def dKdiag_dX(self,partial,X,target):
         pass
 
+    #---------------------------------------#
+    #             PSI statistics            #
+    #---------------------------------------#
+
     def psi0(self, Z, mu, S, target):
         target += self.variance
 
@@ -59,27 +63,27 @@ class bias(kernpart):
     def dpsi0_dtheta(self, partial, Z, mu, S, target):
         target += partial.sum()
 
+    def dpsi1_dtheta(self, partial, Z, mu, S, target):
+        target += partial.sum()
+
+    def dpsi2_dtheta(self, partial, Z, mu, S, target):
+        target += 2.*self.variance*partial.sum()
+
+    
     def dpsi0_dZ(self, partial, Z, mu, S, target):
         pass
 
     def dpsi0_dmuS(self, partial, Z, mu, S, target_mu, target_S):
         pass
 
-    def dpsi1_dtheta(self, partial, Z, mu, S, target):
-        target += partial.sum()
-
     def dpsi1_dZ(self, partial, Z, mu, S, target):
         pass
 
     def dpsi1_dmuS(self, partial, Z, mu, S, target_mu, target_S):
         pass
- 
-    def dpsi2_dtheta(self, partial, Z, mu, S, target):
-        target += np.sum(2.*self.variance*partial)
 
     def dpsi2_dZ(self, partial, Z, mu, S, target):
         pass
 
     def dpsi2_dmuS(self, partial, Z, mu, S, target_mu, target_S):
         pass
-
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index c01ba815..d12febbb 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -3,10 +3,11 @@
 
 
 import numpy as np
+import pylab as pb
 from ..core.parameterised import parameterised
 from kernpart import kernpart
 import itertools
-from product_orthogonal import product_orthogonal 
+from product_orthogonal import product_orthogonal
 
 class kern(parameterised):
     def __init__(self,D,parts=[], input_slices=None):
@@ -386,3 +387,59 @@ class kern(parameterised):
 
         #TODO: there are some extra terms to compute here!
         return target_mu, target_S
+
+    def plot(self, x = None, plot_limits=None,which_functions='all',resolution=None,*args,**kwargs):
+        if which_functions=='all':
+            which_functions = [True]*self.Nparts
+        if self.D == 1:
+            if x is None:
+                x = np.zeros((1,1))
+            else:
+                x = np.asarray(x)
+                assert x.size == 1, "The size of the fixed variable x is not 1"
+                x = x.reshape((1,1))
+
+            if plot_limits == None:
+                xmin, xmax = (x-5).flatten(), (x+5).flatten()
+            elif len(plot_limits) == 2:
+                xmin, xmax = plot_limits
+            else:
+                raise ValueError, "Bad limits for plotting"
+
+            Xnew = np.linspace(xmin,xmax,resolution or 201)[:,None]
+            Kx = self.K(Xnew,x,slices2=which_functions)
+            pb.plot(Xnew,Kx,*args,**kwargs)
+            pb.xlim(xmin,xmax)
+            pb.xlabel("x")
+            pb.ylabel("k(x,%0.1f)" %x)
+
+        elif self.D == 2:
+            if x is None:
+                x = np.zeros((1,2))
+            else:
+                x = np.asarray(x)
+                assert x.size == 2, "The size of the fixed variable x is not 2"
+                x = x.reshape((1,2))
+
+            if plot_limits == None:
+                xmin, xmax = (x-5).flatten(), (x+5).flatten()
+            elif len(plot_limits) == 2:
+                xmin, xmax = plot_limits
+            else:
+                raise ValueError, "Bad limits for plotting"
+
+            resolution = resolution or 51
+            xx,yy = np.mgrid[xmin[0]:xmax[0]:1j*resolution,xmin[1]:xmax[1]:1j*resolution]
+            xg = np.linspace(xmin[0],xmax[0],resolution)
+            yg = np.linspace(xmin[1],xmax[1],resolution)
+            Xnew = np.vstack((xx.flatten(),yy.flatten())).T
+            Kx = self.K(Xnew,x,slices2=which_functions)
+            Kx = Kx.reshape(resolution,resolution).T
+            pb.contour(xg,yg,Kx,vmin=Kx.min(),vmax=Kx.max(),cmap=pb.cm.jet)
+            pb.xlim(xmin[0],xmax[0])
+            pb.ylim(xmin[1],xmax[1])
+            pb.xlabel("x1")
+            pb.ylabel("x2")
+            pb.title("k(x1,x2 ; %0.1f,%0.1f)" %(x[0,0],x[0,1]) )
+        else:
+            raise NotImplementedError, "Cannot plot a kernel with more than two input dimensions"
diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py
index d7869f0a..da4f79f4 100644
--- a/GPy/kern/linear.py
+++ b/GPy/kern/linear.py
@@ -90,22 +90,19 @@ class linear(kernpart):
 
     def psi0(self,Z,mu,S,target):
         self._psi_computations(Z,mu,S)
-        target += np.sum(self.variances*self.mu2_S)
+        target += np.sum(self.variances*self.mu2_S,1)
 
     def dpsi0_dtheta(self,partial,Z,mu,S,target):
         self._psi_computations(Z,mu,S)
-        tmp = (partial[:, None] * (np.sum(self.mu2_S,0)))
+        tmp = partial[:, None] * self.mu2_S
         if self.ARD:
             target += tmp.sum(0)
         else:
             target += tmp.sum()
 
     def dpsi0_dmuS(self,partial, Z,mu,S,target_mu,target_S):
-        target_mu += np.sum(partial[:, None],0) * (2.0*mu*self.variances)
-        target_S += np.sum(partial[:, None] * self.variances, 0)
-
-    def dpsi0_dZ(self,Z,mu,S,target):
-        pass
+        target_mu += partial[:, None] * (2.0*mu*self.variances)
+        target_S += partial[:, None] * self.variances
 
     def psi1(self,Z,mu,S,target):
         """the variance, it does nothing"""
@@ -149,7 +146,7 @@ class linear(kernpart):
     def dpsi2_dZ(self,partial,Z,mu,S,target):
         self._psi_computations(Z,mu,S)
         mu2_S = np.sum(self.mu2_S,0)# Q,
-        target += (partial[:,:,:,None]* (Z * mu2_S * np.square(self.variances))).sum(0).sum(1)
+        target += (partial[:,:,:,None] * (self.mu2_S[:,None,None,:]*(Z*np.square(self.variances)[None,:])[None,None,:,:])).sum(0).sum(1)
 
     #---------------------------------------#
     #            Precomputations            #
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index 5babfa4f..16eda459 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -155,21 +155,20 @@ class rbf(kernpart):
         self._psi_computations(Z,mu,S)
         d_var = 2.*self._psi2/self.variance
         d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom)
-        d_length = d_length.sum(0)
+
         target[0] += np.sum(partial*d_var)
         dpsi2_dlength = d_length*partial[:,:,:,None]
         if not self.ARD:
             target[1] += dpsi2_dlength.sum()
         else:
             target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)
-
+            
     def dpsi2_dZ(self,partial,Z,mu,S,target):
         self._psi_computations(Z,mu,S)
         term1 = 0.5*self._psi2_Zdist/self.lengthscale2 # M, M, Q
         term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, Q
         dZ = self._psi2[:,:,:,None] * (term1[None] + term2)
-        target += (partial[:,:,:,None]*dZ).sum(0).sum(0) # <----------------- TODO not sure about the first ':' here, should be a None (WAS a none in the debug branch)
-
+        target += (partial[:,:,:,None]*dZ).sum(0).sum(0)
 
     def dpsi2_dmuS(self,partial,Z,mu,S,target_mu,target_S):
         """Think N,M,M,Q """
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 66d80d11..73d9416a 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -34,7 +34,7 @@ class sparse_GP(GP):
     """
 
     def __init__(self, X, likelihood, kernel, Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False):
-        self.scale_factor = 1.0# a scaling factor to help keep the algorithm stable
+        self.scale_factor = 100.0# a scaling factor to help keep the algorithm stable
 
         self.Z = Z
         self.Zslices = Zslices
diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py
new file mode 100644
index 00000000..c49bdfda
--- /dev/null
+++ b/GPy/testing/bgplvm_tests.py
@@ -0,0 +1,48 @@
+# Copyright (c) 2012, Nicolo Fusi
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import unittest
+import numpy as np
+import GPy
+
+class BGPLVMTests(unittest.TestCase):
+    def test_bias_kern(self):
+        N, M, Q, D = 10, 3, 2, 4
+        X = np.random.rand(N, Q)
+        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        K = k.K(X)
+        Y = np.random.multivariate_normal(np.zeros(N),K,D).T
+        k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        m.constrain_positive('(rbf|bias|noise|white|S)')
+        m.randomize()
+        self.assertTrue(m.checkgrad())
+
+    def test_linear_kern(self):
+        N, M, Q, D = 10, 3, 2, 4
+        X = np.random.rand(N, Q)
+        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        K = k.K(X)
+        Y = np.random.multivariate_normal(np.zeros(N),K,D).T
+        k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        m.constrain_positive('(linear|bias|noise|white|S)')
+        m.randomize()
+        self.assertTrue(m.checkgrad())
+
+    def test_rbf_kern(self):
+        N, M, Q, D = 10, 3, 2, 4
+        X = np.random.rand(N, Q)
+        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        K = k.K(X)
+        Y = np.random.multivariate_normal(np.zeros(N),K,D).T
+        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        m.constrain_positive('(rbf|bias|noise|white|S)')
+        m.randomize()
+        self.assertTrue(m.checkgrad())
+
+        
+if __name__ == "__main__":
+    print "Running unit tests, please be (very) patient..."
+    unittest.main()