diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py
index 464d7425..6382fbde 100644
--- a/GPy/models/Bayesian_GPLVM.py
+++ b/GPy/models/Bayesian_GPLVM.py
@@ -13,6 +13,7 @@ from numpy.linalg.linalg import LinAlgError
 import itertools
 from matplotlib.colors import colorConverter
 from matplotlib.figure import SubplotParams
+from GPy.inference.optimization import SCG
 
 class Bayesian_GPLVM(sparse_GP, GPLVM):
     """
@@ -184,16 +185,46 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         ax.plot(self.Z[:, input_1], self.Z[:, input_2], '^w')
         return ax
 
+    def do_test_latents(self, Y):
+        """
+        Compute the latent representation for a set of new points Y
+
+        Notes:
+        This will only work with a univariate Gaussian likelihood (for now)
+        """
+        assert not self.likelihood.is_heteroscedastic
+        N_test = Y.shape[0]
+        Q = self.Z.shape[1]
+        means = np.zeros((N_test,Q))
+        covars = np.zeros((N_test,Q))
+
+        dpsi0 = - 0.5 * self.D * self.likelihood.precision
+        dpsi2 = self.dL_dpsi2[0][None,:,:] # TODO: this may change if we ignore het. likelihoods
+        V = self.likelihood.precision*Y
+        dpsi1 = np.dot(self.Cpsi1V,V.T)
+
+        start = np.zeros(self.Q*2)
+
+        for n,dpsi1_n in enumerate(dpsi1.T[:,:,None]):
+            args = (self.kern,self.Z,dpsi0,dpsi1_n,dpsi2)
+            xopt,fopt,neval,status = SCG(f=latent_cost, gradf=latent_grad, x=start, optargs=args, display = False)
+
+            mu,log_S = xopt.reshape(2,1,-1)
+            means[n] = mu[0].copy()
+            covars[n] = np.exp(log_S[0]).copy()
+
+        return means, covars
+
+
     def plot_X_1d(self, fig=None, axes=None, fig_num="LVM mu S 1d", colors=None):
         """
         Plot latent space X in 1D:
-        
+
             -if fig is given, create Q subplots in fig and plot in these
             -if axes is given plot Q 1D latent space plots of X into each `axis`
             -if neither fig nor axes is given create a figure with fig_num and plot in there
-            
+
         colors:
-            
             colors of different latent space dimensions Q
         """
         import pylab
@@ -483,3 +514,63 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         cidd = figs[0].canvas.mpl_connect('motion_notify_event', motion)
 
         return ax1, ax2, ax3, ax4, ax5 # , ax6, ax7
+
+
+
+
+def latent_cost_and_grad(mu_S, kern,Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
+    """
+    objective function for fitting the latent variables for test points
+    (negative log-likelihood: should be minimised!)
+    """
+    mu,log_S = mu_S.reshape(2,1,-1)
+    S = np.exp(log_S)
+
+    psi0 = kern.psi0(Z,mu,S)
+    psi1 = kern.psi1(Z,mu,S)
+    psi2 = kern.psi2(Z,mu,S)
+
+    lik = dL_dpsi0*psi0 + np.dot(dL_dpsi1.flatten(),psi1.flatten()) + np.dot(dL_dpsi2.flatten(),psi2.flatten()) - 0.5*np.sum(np.square(mu) + S) + 0.5*np.sum(log_S)
+
+    mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0,Z,mu,S)
+    mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1,Z,mu,S)
+    mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2,Z,mu,S)
+
+    dmu = mu0 + mu1 + mu2 - mu
+    #dS = S0 + S1 + S2 -0.5 + .5/S
+    dlnS = S*(S0 + S1 + S2 -0.5) + .5
+    return -lik,-np.hstack((dmu.flatten(),dlnS.flatten()))
+
+def latent_cost(mu_S, kern,Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
+    """
+    objective function for fitting the latent variables (negative log-likelihood: should be minimised!)
+    This is the same as latent_cost_and_grad but only for the objective
+    """
+    mu,log_S = mu_S.reshape(2,1,-1)
+    S = np.exp(log_S)
+
+    psi0 = kern.psi0(Z,mu,S)
+    psi1 = kern.psi1(Z,mu,S)
+    psi2 = kern.psi2(Z,mu,S)
+
+    lik = dL_dpsi0*psi0 + np.dot(dL_dpsi1.flatten(),psi1.flatten()) + np.dot(dL_dpsi2.flatten(),psi2.flatten()) - 0.5*np.sum(np.square(mu) + S) + 0.5*np.sum(log_S)
+    return -float(lik)
+
+def latent_grad(mu_S, kern,Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
+    """
+    This is the same as latent_cost_and_grad but only for the grad
+    """
+    mu,log_S = mu_S.reshape(2,1,-1)
+    S = np.exp(log_S)
+
+    mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0,Z,mu,S)
+    mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1,Z,mu,S)
+    mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2,Z,mu,S)
+
+    dmu = mu0 + mu1 + mu2 - mu
+    #dS = S0 + S1 + S2 -0.5 + .5/S
+    dlnS = S*(S0 + S1 + S2 -0.5) + .5
+
+    return -np.hstack((dmu.flatten(),dlnS.flatten()))
+
+
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 259db8f2..0d9e5a3f 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -233,3 +233,5 @@ class sparse_GP(GP):
             var = Kxx - np.sum(Kx * np.dot(Kmmi_LmiBLmi, Kx), 0)
 
         return mu, var[:, None]
+
+