diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 61a28448..09cc9abe 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -378,10 +378,28 @@ def stick(kernel=None):
 
     return m
 
+def bcgplvm_linear_stick(kernel=None):
+    data = GPy.util.datasets.osu_run1()
+    # optimize
+    mapping = GPy.mappings.Linear(data['Y'].shape[1], 2)
+    m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
+    m.optimize(messages=1, max_f_eval=10000)
+    if GPy.util.visualize.visual_available:
+        plt.clf
+        ax = m.plot_latent()
+        y = m.likelihood.Y[0, :]
+        data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
+        lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
+        raw_input('Press enter to finish')
+
+    return m
+
 def bcgplvm_stick(kernel=None):
     data = GPy.util.datasets.osu_run1()
     # optimize
-    m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel)
+    back_kernel=GPy.kern.rbf(data['Y'].shape[1], lengthscale=10.)
+    mapping = GPy.mappings.Kernel(X=data['Y'], output_dim=2, kernel=back_kernel)
+    m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
     m.optimize(messages=1, max_f_eval=10000)
     if GPy.util.visualize.visual_available:
         plt.clf
diff --git a/GPy/models/bcgplvm.py b/GPy/models/bcgplvm.py
new file mode 100644
index 00000000..b6246c32
--- /dev/null
+++ b/GPy/models/bcgplvm.py
@@ -0,0 +1,50 @@
+# ## Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+import numpy as np
+import pylab as pb
+import sys, pdb
+from ..core import GP
+from ..models import GPLVM
+from ..mappings import *
+
+
+class BCGPLVM(GPLVM):
+    """
+    Back constrained Gaussian Process Latent Variable Model
+
+    :param Y: observed data
+    :type Y: np.ndarray
+    :param input_dim: latent dimensionality
+    :type input_dim: int
+    :param init: initialisation method for the latent space
+    :type init: 'PCA'|'random'
+    :param mapping: mapping for back constraint
+    :type mapping: GPy.core.Mapping object
+
+    """
+    def __init__(self, Y, input_dim, init='PCA', X=None, kernel=None, normalize_Y=False, mapping=None):
+        
+        if mapping is None:
+            mapping = Kernel(X=Y, output_dim=input_dim)
+        self.mapping = mapping
+        GPLVM.__init__(self, Y, input_dim, init, X, kernel, normalize_Y)
+        self.X = self.mapping.f(self.likelihood.Y)
+
+    def _get_param_names(self):
+        return self.mapping._get_param_names() + GP._get_param_names(self)
+
+    def _get_params(self):
+        return np.hstack((self.mapping._get_params(), GP._get_params(self)))
+
+    def _set_params(self, x):
+        self.mapping._set_params(x[:self.mapping.num_params])
+        self.X = self.mapping.f(self.likelihood.Y)
+        GP._set_params(self, x[self.mapping.num_params:])
+
+    def _log_likelihood_gradients(self):
+        dL_df = 2.*self.kern.dK_dX(self.dL_dK, self.X)
+        dL_dtheta = self.mapping.df_dtheta(dL_df, self.likelihood.Y)
+        return np.hstack((dL_dtheta.flatten(), GP._log_likelihood_gradients(self)))
+