diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py
index 71921ab1..8bc7ca59 100644
--- a/GPy/core/parameterization/variational.py
+++ b/GPy/core/parameterization/variational.py
@@ -117,4 +117,4 @@ class SpikeAndSlabPosterior(VariationalPosterior):
         import sys
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ...plotting.matplot_dep import variational_plots
-        return variational_plots.plot(self,*args)
+        return variational_plots.plot_SpikeSlab(self,*args)
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 9ebb54a2..818dff69 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -515,3 +515,28 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
         lvm_visualizer.close()
 
     return m
+
+def ssgplvm_simulation_linear():
+    import numpy as np
+    import GPy
+    N, D, Q = 1000, 20, 5
+    pi = 0.2
+    
+    def sample_X(Q, pi):
+        x = np.empty(Q)
+        dies = np.random.rand(Q)
+        for q in xrange(Q):
+            if dies[q]<pi:
+                x[q] = np.random.randn()
+            else:
+                x[q] = 0.
+        return x
+    
+    Y = np.empty((N,D))
+    X = np.empty((N,Q))
+    # Generate data from random sampled weight matrices
+    for n in xrange(N):
+        X[n] = sample_X(Q,pi)
+        w = np.random.randn(D,Q)
+        Y[n] = np.dot(w,X[n])
+    
diff --git a/GPy/plotting/matplot_dep/variational_plots.py b/GPy/plotting/matplot_dep/variational_plots.py
index 72b857a6..cf00d8a2 100644
--- a/GPy/plotting/matplot_dep/variational_plots.py
+++ b/GPy/plotting/matplot_dep/variational_plots.py
@@ -44,3 +44,48 @@ def plot(parameterized, fignum=None, ax=None, colors=None):
     pb.draw()
     fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
     return fig
+
+def plot_SpikeSlab(parameterized, fignum=None, ax=None, colors=None):
+    """
+    Plot latent space X in 1D:
+
+        - if fig is given, create input_dim subplots in fig and plot in these
+        - if ax is given plot input_dim 1D latent space plots of X into each `axis`
+        - if neither fig nor ax is given create a figure with fignum and plot in there
+
+    colors:
+        colors of different latent space dimensions input_dim
+
+    """
+    if ax is None:
+        fig = pb.figure(num=fignum, figsize=(8, min(12, (2 * parameterized.mean.shape[1]))))
+    if colors is None:
+        colors = pb.gca()._get_lines.color_cycle
+        pb.clf()
+    else:
+        colors = iter(colors)
+    plots = []
+    means, variances, gamma = param_to_array(parameterized.mean, parameterized.variance, parameterized.binary_prob)
+    x = np.arange(means.shape[0])
+    for i in range(means.shape[1]):
+        # mean and variance plot
+        a = fig.add_subplot(means.shape[1]*2, 1, 2*i + 1)
+        a.plot(means, c='k', alpha=.3)
+        plots.extend(a.plot(x, means.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
+        a.fill_between(x,
+                        means.T[i] - 2 * np.sqrt(variances.T[i]),
+                        means.T[i] + 2 * np.sqrt(variances.T[i]),
+                        facecolor=plots[-1].get_color(),
+                        alpha=.3)
+        a.legend(borderaxespad=0.)
+        a.set_xlim(x.min(), x.max())
+        if i < means.shape[1] - 1:
+            a.set_xticklabels('')
+        # binary prob plot
+        a = fig.add_subplot(means.shape[1]*2, 1, 2*i + 2)
+        a.bar(x,gamma[:,i],bottom=0.,linewidth=0,align='center')
+        a.set_xlim(x.min(), x.max())
+        a.set_ylim([0.,1.])
+    pb.draw()
+    fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
+    return fig