diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py
index 8f702669..23bbbe6b 100644
--- a/GPy/plotting/__init__.py
+++ b/GPy/plotting/__init__.py
@@ -29,15 +29,18 @@ if config.get('plotting', 'library') is not 'none':
     
     from ..core import GP
     GP.plot_data = gpy_plot.data_plots.plot_data
+    GP.plot_errorbars_trainset = gpy_plot.data_plots.plot_errorbars_trainset
     GP.plot_mean = gpy_plot.gp_plots.plot_mean
     GP.plot_confidence = gpy_plot.gp_plots.plot_confidence
-
+    GP.plot_density = gpy_plot.gp_plots.plot_density
+    
+    from ..core import SparseGP
+    SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing
+    
     from . import matplot_dep
     # Still to convert to new style:
     GP.plot = matplot_dep.models_plots.plot_fit
     GP.plot_f = matplot_dep.models_plots.plot_fit_f
-    GP.plot_density = matplot_dep.models_plots.plot_density
     
-    GP.plot_errorbars_trainset = matplot_dep.models_plots.plot_errorbars_trainset
     GP.plot_magnification = matplot_dep.dim_reduction_plots.plot_magnification
         
diff --git a/GPy/plotting/abstract_plotting_library.py b/GPy/plotting/abstract_plotting_library.py
index 7bfc662f..8911e8de 100644
--- a/GPy/plotting/abstract_plotting_library.py
+++ b/GPy/plotting/abstract_plotting_library.py
@@ -87,18 +87,35 @@ class AbstractPlottingLibrary(object):
         the kwargs are plotting library specific kwargs!
         """
         raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
-    
-    def scatter(self, canvas, X, Y, **kwargs):
+
+    def plot_axis_lines(self, ax, X, **kwargs):
+        """
+        Plot lines at the bottom of the axis at input location X.
+        
+        the kwargs are plotting library specific kwargs!
+        """
+        raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
+
+    def scatter(self, canvas, X, Y, c=None, vmin=None, vmax=None, **kwargs):
         """
         Make a scatter plot between X and Y on the canvas given.
         
         the kwargs are plotting library specific kwargs!
+        
+        :param canvas: the plotting librarys specific canvas to plot on.
+        :param array-like X: the inputs to plot.
+        :param array-like Y: the outputs to plot.
+        :param array-like c: the colorlevel for each point.
+        :param float vmin: minimum colorscale
+        :param float vmax: maximum colorscale
         """
         raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
 
     def xerrorbar(self, canvas, X, Y, error, **kwargs):
         """
         Make an errorbar along the xaxis for points at (X,Y) on the canvas.
+        if error is two dimensional, the lower error is error[:,0] and
+        the upper error is error[:,1]
         
         the kwargs are plotting library specific kwargs!
         """
@@ -107,6 +124,8 @@ class AbstractPlottingLibrary(object):
     def yerrorbar(self, canvas, X, Y, error, **kwargs):
         """
         Make errorbars along the yaxis on the canvas given.
+        if error is two dimensional, the lower error is error[:,0] and
+        the upper error is error[:,1]
                 
         the kwargs are plotting library specific kwargs!
         """
diff --git a/GPy/plotting/gpy_plot/data_plots.py b/GPy/plotting/gpy_plot/data_plots.py
index eb0cdf46..1a94b4a4 100644
--- a/GPy/plotting/gpy_plot/data_plots.py
+++ b/GPy/plotting/gpy_plot/data_plots.py
@@ -32,7 +32,7 @@ from . import pl
 
 import numpy as np
 from .plot_util import get_x_y_var, get_free_dims, get_which_data_ycols,\
-    get_which_data_rows, update_not_existing_kwargs
+    get_which_data_rows, update_not_existing_kwargs, helper_predict_with_model
 
 def _plot_data(self, canvas, which_data_rows='all',
         which_data_ycols='all', visible_dims=None,
@@ -65,6 +65,8 @@ def _plot_data(self, canvas, which_data_rows='all',
             update_not_existing_kwargs(plot_kwargs, pl.defaults.data_2d)
             plots['dataplot'].append(pl.scatter(canvas, X[rows, free_dims[0]], X[rows, free_dims[1]], 
                                            c=Y[rows, d], vmin=Y.min(), vmax=Y.max(), **plot_kwargs))
+    elif len(free_dims) == 0:
+        pass #Nothing to plot!
     else:
         raise NotImplementedError("Cannot plot in more then two dimensions")
     return plots
@@ -93,3 +95,100 @@ def plot_data(self, which_data_rows='all',
     canvas, kwargs = pl.get_new_canvas(plot_kwargs)
     plots = _plot_data(self, canvas, which_data_rows, which_data_ycols, visible_dims, error_kwargs, **kwargs)
     return pl.show_canvas(canvas, plots)
+
+
+def plot_inducing(self, visible_dims=None, **plot_kwargs):
+    """
+    Plot the inducing inputs of a sparse gp model
+    
+    :param array-like visible_dims: an array specifying the input dimensions to plot (maximum two)
+    :param kwargs plot_kwargs: keyword arguments for the plotting library
+    """
+    canvas, kwargs = pl.get_new_canvas(plot_kwargs)
+    plots = _plot_inducing(self, canvas, visible_dims, **kwargs)
+    return pl.show_canvas(canvas, plots)
+
+def _plot_inducing(self, canvas, visible_dims, **plot_kwargs):
+    free_dims = get_free_dims(self, visible_dims, None)
+
+    Z = self.Z[:, free_dims]
+    plots = {}
+
+    #one dimensional plotting
+    if len(free_dims) == 1:
+        update_not_existing_kwargs(plot_kwargs, pl.defaults.inducing_1d)
+        plots['inducing'] = pl.plot_axis_lines(canvas, Z[:, free_dims], **plot_kwargs)
+    #2D plotting
+    elif len(free_dims) == 2:
+        update_not_existing_kwargs(plot_kwargs, pl.defaults.inducing_2d)
+        plots['inducing'] = pl.scatter(canvas, Z[:, free_dims[0]], Z[:, free_dims[1]], 
+                                       **plot_kwargs)
+    elif len(free_dims) == 0:
+        pass #Nothing to plot!
+    else:
+        raise NotImplementedError("Cannot plot in more then two dimensions")
+    return plots
+
+def plot_errorbars_trainset(self, which_data_rows='all',
+        which_data_ycols='all', fixed_inputs=None, 
+        plot_raw=False, apply_link=False,
+        predict_kw=None, **plot_kwargs):
+    """
+    Plot the errorbars of the GP likelihood on the training data.
+    These are the errorbars after the appropriate 
+    approximations according to the likelihood are done.
+    
+    This also works for heteroscedastic likelihoods.
+    
+    Give the Y_metadata in the predict_kw if you need it.
+    
+    :param which_data_rows: which of the training data to plot (default all)
+    :type which_data_rows: 'all' or a slice object to slice self.X, self.Y
+    :param which_data_ycols: when the data has several columns (independant outputs), only plot these
+    :param fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input dimension i should be set to value v.
+    :type fixed_inputs: a list of tuples
+    :param dict predict_kwargs: kwargs for the prediction used to predict the right quantiles.
+    :param kwargs plot_kwargs: kwargs for the data plot for the plotting library you are using
+    """
+    canvas, kwargs = pl.get_new_canvas(plot_kwargs)
+    plots = _plot_errorbars_trainset(self, canvas, which_data_rows, which_data_ycols, 
+                                     fixed_inputs, plot_raw, apply_link, predict_kw, **kwargs)
+    return pl.show_canvas(canvas, plots)
+
+def _plot_errorbars_trainset(self, canvas, 
+        which_data_rows='all', which_data_ycols='all', 
+        fixed_inputs=None,
+        plot_raw=False, apply_link=False,
+        predict_kw=None, **plot_kwargs):
+
+    ycols = get_which_data_ycols(self, which_data_ycols)
+    rows = get_which_data_rows(self, which_data_rows)
+
+    X, _, Y = get_x_y_var(self)
+    
+    if fixed_inputs is None:
+        fixed_inputs = []
+    free_dims = get_free_dims(self, None, fixed_inputs)    
+
+    Xgrid = X.copy()
+    for i, v in fixed_inputs:
+        Xgrid[:, i] = v
+
+    plots = []
+    
+    if len(free_dims)<2:
+        if len(free_dims)==1:
+            update_not_existing_kwargs(plot_kwargs, pl.defaults.yerrorbar)
+            _, percs = helper_predict_with_model(self, Xgrid, plot_raw, 
+                                              apply_link, (2.5, 97.5), 
+                                              ycols, predict_kw)
+            for d in ycols:
+                plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], Y[rows,d], 
+                                          np.vstack([Y[rows,d]-percs[0][rows,d], percs[1][rows,d]-Y[rows,d]]),
+                                          **plot_kwargs))
+            return dict(yerrorbars=plots)
+        else:
+            pass #Nothing to plot!
+    else:
+        raise NotImplementedError("Cannot plot in more then one dimension.")
+    return plots
\ No newline at end of file
diff --git a/GPy/plotting/gpy_plot/gp_plots.py b/GPy/plotting/gpy_plot/gp_plots.py
index 7a1f848a..f43ee265 100644
--- a/GPy/plotting/gpy_plot/gp_plots.py
+++ b/GPy/plotting/gpy_plot/gp_plots.py
@@ -32,40 +32,8 @@ import numpy as np
 from functools import wraps
 
 from . import pl
-from .plot_util import get_x_y_var, get_fixed_dims, get_free_dims, \
-    x_frame1D, x_frame2D, update_not_existing_kwargs, \
-    helper_predict_with_model
-
-def _helper_for_plots(self, plot_limits, fixed_inputs, resolution):
-    """
-    Figure out the data, free_dims and create an Xgrid for
-    the prediction. 
-    """
-    X, Xvar, Y = get_x_y_var(self)
-
-    #work out what the inputs are for plotting (1D or 2D)
-    fixed_dims = get_fixed_dims(self, fixed_inputs)
-    free_dims = get_free_dims(self, None, fixed_dims)
-    
-    if len(free_dims) == 1:
-        #define the frame on which to plot
-        resolution = resolution or 200
-        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits, resolution=resolution)
-        Xgrid = np.empty((Xnew.shape[0],self.input_dim))
-        Xgrid[:,free_dims] = Xnew
-        for i,v in fixed_dims:
-            Xgrid[:,i] = v
-        x = Xgrid
-        y = None
-    elif len(free_dims) == 2:
-        #define the frame for plotting on
-        resolution = resolution or 50
-        Xnew, x, y, xmin, xmax = x_frame2D(X[:,free_dims], plot_limits, resolution)
-        Xgrid = np.empty((Xnew.shape[0],self.input_dim))
-        Xgrid[:,free_dims] = Xnew
-        for i,v in fixed_dims:
-            Xgrid[:,i] = v    
-    return X, Xvar, Y, fixed_dims, free_dims, Xgrid, x, y, xmin, xmax, resolution
+from .plot_util import helper_for_plot_data, update_not_existing_kwargs, \
+    helper_predict_with_model, get_which_data_ycols
 
 def plot_mean(self, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=False,
@@ -76,6 +44,9 @@ def plot_mean(self, plot_limits=None, fixed_inputs=None,
               **kwargs):
     """
     Plot the mean of the GP.
+
+    Give the Y_metadata in the predict_kw if you need it.
+   
     
     :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
     :type plot_limits: np.array
@@ -90,7 +61,10 @@ def plot_mean(self, plot_limits=None, fixed_inputs=None,
     :param int levels: for 2D plotting, the number of contour levels to use is 
     """
     canvas, kwargs = pl.get_new_canvas(kwargs)
-    plots = _plot_mean(self, canvas, plot_limits, fixed_inputs, resolution, plot_raw, Y_metadata, apply_link, which_data_ycols, levels, predict_kw, **kwargs)
+    plots = _plot_mean(self, canvas, plot_limits, fixed_inputs, 
+                       resolution, plot_raw, Y_metadata, 
+                       apply_link, which_data_ycols, levels, 
+                       predict_kw, **kwargs)
     return pl.show_canvas(canvas, plots)
 
 @wraps(plot_mean)
@@ -100,27 +74,30 @@ def _plot_mean(self, canvas, plot_limits=None, fixed_inputs=None,
               which_data_ycols=None,
               levels=20, 
               predict_kw=None, **kwargs):
-    if predict_kw is None:
-        predict_kw = {}
+    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, fixed_inputs, resolution)
 
-    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = _helper_for_plots(self, plot_limits, fixed_inputs, resolution)
-
-    if len(free_dims<=2):
-        which_data_ycols = get_which_data_ycols(self, which_data_ycols)
-        mu, _ = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, None, which_data_ycols, **predict_kw)
+    if len(free_dims)<=2:
+        mu, _ = helper_predict_with_model(self, Xgrid, plot_raw, 
+                                          apply_link, None, 
+                                          get_which_data_ycols(self, which_data_ycols), 
+                                          predict_kw)
         if len(free_dims)==1:
             # 1D plotting:
             update_not_existing_kwargs(kwargs, pl.defaults.meanplot_1d)
-            return dict(gpmean=[pl.plot(canvas, Xgrid, mu, **kwargs)])
+            return dict(gpmean=[pl.plot(canvas, Xgrid[:, free_dims], mu, **kwargs)])
         else:
             update_not_existing_kwargs(kwargs, pl.defaults.meanplot_2d)
             return dict(gpmean=[pl.contour(canvas, x, y, 
                                            mu.reshape(resolution, resolution), 
                                            levels=levels, **kwargs)])
+    elif len(free_dims)==0:
+        pass # Nothing to plot!
+    else:
+        raise RuntimeError('Cannot plot mean in more then 2 input dimensions')
 
 def plot_confidence(self, lower=2.5, upper=97.5, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=False,
-              Y_metadata=None, apply_link=False, 
+              apply_link=False, 
               which_data_ycols='all',
               predict_kw=None, 
               **kwargs):
@@ -128,6 +105,9 @@ def plot_confidence(self, lower=2.5, upper=97.5, plot_limits=None, fixed_inputs=
     Plot the confidence interval between the percentiles lower and upper.
     E.g. the 95% confidence interval is $2.5, 97.5$.
     Note: Only implemented for one dimension!
+
+    Give the Y_metadata in the predict_kw if you need it.
+   
     
     :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
     :type plot_limits: np.array
@@ -142,26 +122,99 @@ def plot_confidence(self, lower=2.5, upper=97.5, plot_limits=None, fixed_inputs=
     """
     canvas, kwargs = pl.get_new_canvas(kwargs)
     plots = _plot_confidence(self, canvas, lower, upper, plot_limits, 
-                             fixed_inputs, resolution, plot_raw, Y_metadata, 
+                             fixed_inputs, resolution, plot_raw, 
                              apply_link, which_data_ycols, 
                              predict_kw, **kwargs)
     return pl.show_canvas(canvas, plots)
 
 def _plot_confidence(self, canvas, lower, upper, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=False,
-              Y_metadata=None, apply_link=False, 
+              apply_link=False, 
               which_data_ycols=None, 
               predict_kw=None, 
               **kwargs):
-    if predict_kw is None:
-        predict_kw = {}
-    
-    _, _, _, _, _, Xgrid, _, _, _, _, _ = _helper_for_plots(self, plot_limits, fixed_inputs, resolution)
+    _, _, _, _, free_dims, Xgrid, _, _, _, _, _ = helper_for_plot_data(self, plot_limits, fixed_inputs, resolution)
 
+    ycols = get_which_data_ycols(self, which_data_ycols)
+    
     update_not_existing_kwargs(kwargs, pl.defaults.confidence_interval)
-    _, percs = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, (lower, upper), which_data_ycols, **predict_kw)
-
-    return dict(gpconfidence=pl.fill_between(canvas, Xgrid, percs[0], percs[1], **kwargs))
-
     
-    
\ No newline at end of file
+    if len(free_dims)<=1:
+        if len(free_dims)==1:
+            _, percs = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, 
+                                                 (lower, upper), 
+                                                 ycols, predict_kw)
+    
+            fills = []
+            for d in ycols:
+                fills.append(pl.fill_between(canvas, Xgrid[:,free_dims[0]], percs[0][:,d], percs[1][:,d], **kwargs))
+            return dict(gpconfidence=fills)
+        else:
+            pass #Nothing to plot!
+    else:
+        raise RuntimeError('Can only plot confidence interval in one input dimension')
+
+
+def plot_density(self, plot_limits=None, fixed_inputs=None,
+              resolution=None, plot_raw=False,
+              apply_link=False, 
+              which_data_ycols='all',
+              levels=20,
+              predict_kw=None, 
+              **kwargs):
+    """
+    Plot the confidence interval between the percentiles lower and upper.
+    E.g. the 95% confidence interval is $2.5, 97.5$.
+    Note: Only implemented for one dimension!
+
+    Give the Y_metadata in the predict_kw if you need it.
+   
+    :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
+    :type plot_limits: np.array
+    :param fixed_inputs: a list of tuple [(i,v), (i,v)...], specifying that input dimension i should be set to value v.
+    :type fixed_inputs: a list of tuples
+    :param int resolution: The resolution of the prediction [default:200]
+    :param bool plot_raw: plot the latent function (usually denoted f) only?
+    :param dict Y_metadata: the Y_metadata (for e.g. heteroscedastic GPs)
+    :param bool apply_link: whether to apply the link function of the GP to the raw prediction.
+    :param array-like which_data_ycols: which columns of y to plot (array-like or list of ints)
+    :param int levels: the number of levels in the density (number between 1 and 50, where 50 is very smooth and 1 is the same as plot_confidence) 
+    :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=<specific kernel>) in here
+    """
+    canvas, kwargs = pl.get_new_canvas(kwargs)
+    plots = _plot_density(self, canvas, plot_limits, 
+                             fixed_inputs, resolution, plot_raw,  
+                             apply_link, which_data_ycols, 
+                             levels,
+                             predict_kw, **kwargs)
+    return pl.show_canvas(canvas, plots)
+
+def _plot_density(self, canvas, plot_limits=None, fixed_inputs=None,
+              resolution=None, plot_raw=False,
+              apply_link=False, 
+              which_data_ycols=None,
+              levels=20, 
+              predict_kw=None, **kwargs):
+    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, fixed_inputs, resolution)
+
+    ycols = get_which_data_ycols(self, which_data_ycols)
+
+    update_not_existing_kwargs(kwargs, pl.defaults.density)
+
+    if len(free_dims)<=1:
+        if len(free_dims)==1:
+            _, percs = helper_predict_with_model(self, Xgrid, plot_raw, 
+                                          apply_link, np.linspace(2.5, 97.5, levels*2), 
+                                          get_which_data_ycols(self, which_data_ycols), 
+                                          predict_kw)
+            # 1D plotting:
+            fills = []
+            for d in ycols:
+                fills.append(pl.fill_gradient(canvas, Xgrid[:, free_dims[0]], [p[:,d] for p in percs], **kwargs))
+            return dict(gpdensity=fills)
+        else:
+            pass # Nothing to plot!
+    else:
+        raise RuntimeError('Can only plot density in one input dimension')
+
+            
\ No newline at end of file
diff --git a/GPy/plotting/gpy_plot/plot_util.py b/GPy/plotting/gpy_plot/plot_util.py
index 7f286399..04e4fa0e 100644
--- a/GPy/plotting/gpy_plot/plot_util.py
+++ b/GPy/plotting/gpy_plot/plot_util.py
@@ -31,7 +31,7 @@
 import numpy as np
 from scipy import sparse
 
-def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, which_data_ycols, percentiles, **predict_kw):
+def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, which_data_ycols, predict_kw):
     """
     Make the right decisions for prediction with a model 
     based on the standard arguments of plotting.
@@ -39,29 +39,73 @@ def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, which_data_ycol
     This is quite complex and will take a while to understand,
     so do not change anything in here lightly!!! 
     """
+    # Put some standards into the predict_kw so that prediction is done automatically:
+    if predict_kw is None:
+        predict_kw = {}
     if 'likelihood' not in predict_kw:
         if plot_raw:
             from ...likelihoods import Gaussian
-            lik = Gaussian(0) # Make the likelihood not add any noise
+            from ...likelihoods.link_functions import Identity
+            lik = Gaussian(Identity(), 0) # Make the likelihood not add any noise
         else:
             lik = None
         predict_kw['likelihood'] = lik
+    if 'Y_metadata' not in predict_kw:
+        predict_kw['Y_metadata'] = self.Y_metadata or {}
+    if 'output_index' not in predict_kw['Y_metadata']:
+        predict_kw['Y_metadata']['output_index'] = Xgrid[:,-1:].astype(np.int)
 
     mu, _ = self.predict(Xgrid, **predict_kw)
     
     if percentiles is not None:
         percentiles = self.predict_quantiles(Xgrid, quantiles=percentiles, **predict_kw)
-    else: percentiles = {}
+    else: percentiles = []
 
-    retmu = np.empty((Xgrid.shape[0], len(ycols)))
+    # Filter out the ycolums which we want to plot:
+    retmu = mu[:, which_data_ycols]
+    percs = [p[:, which_data_ycols] for p in percentiles]
     
     if plot_raw and apply_link:
-        for i, d in enumerate(ycols):
-            retmu = self.likelihood.gp_link.transf(mu[:, [i]])
-            for perc in percentiles:
+        for i in range(len(which_data_ycols)):
+            retmu[:, [i]] = self.likelihood.gp_link.transf(mu[:, [i]])
+            for perc in percs:
                 perc[:, [i]] = self.likelihood.gp_link.transf(perc[:, [i]])
 
-    return mu, percentiles
+    return retmu, percs
+
+def helper_for_plot_data(self, plot_limits, fixed_inputs, resolution):
+    """
+    Figure out the data, free_dims and create an Xgrid for
+    the prediction. 
+    """
+    X, Xvar, Y = get_x_y_var(self)
+
+    #work out what the inputs are for plotting (1D or 2D)
+    if fixed_inputs is None:
+        fixed_inputs = []
+    fixed_dims = get_fixed_dims(self, fixed_inputs)
+    free_dims = get_free_dims(self, None, fixed_dims)
+    
+    if len(free_dims) == 1:
+        #define the frame on which to plot
+        resolution = resolution or 200
+        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits, resolution=resolution)
+        Xgrid = np.empty((Xnew.shape[0],self.input_dim))
+        Xgrid[:,free_dims] = Xnew
+        for i,v in fixed_inputs:
+            Xgrid[:,i] = v
+        x = Xgrid
+        y = None
+    elif len(free_dims) == 2:
+        #define the frame for plotting on
+        resolution = resolution or 50
+        Xnew, x, y, xmin, xmax = x_frame2D(X[:,free_dims], plot_limits, resolution)
+        Xgrid = np.empty((Xnew.shape[0],self.input_dim))
+        Xgrid[:,free_dims] = Xnew
+        for i,v in fixed_inputs:
+            Xgrid[:,i] = v    
+    return X, Xvar, Y, fixed_dims, free_dims, Xgrid, x, y, xmin, xmax, resolution
+
 
 def update_not_existing_kwargs(to_update, update_from):
     """
@@ -112,8 +156,6 @@ def get_fixed_dims(model, fixed_inputs):
     """
     Work out the fixed dimensions from the fixed_inputs list of tuples.
     """
-    if fixed_inputs is None:
-        fixed_inputs = []
     return np.array([i for i,_ in fixed_inputs])
 
 def get_which_data_ycols(model, which_data_ycols):
diff --git a/GPy/plotting/matplot_dep/defaults.py b/GPy/plotting/matplot_dep/defaults.py
index c8db4b19..a122783b 100644
--- a/GPy/plotting/matplot_dep/defaults.py
+++ b/GPy/plotting/matplot_dep/defaults.py
@@ -43,11 +43,14 @@ from . import Tango
 
 # Data:
 data_1d = dict(lw=1.5, marker='x', edgecolor='k')
-data_2d = dict(s=35, edgecolors='none', linewidth=0., cmap=cm.get_cmap('hot'))
+data_2d = dict(s=35, edgecolors='none', linewidth=0., cmap=cm.get_cmap('hot'), alpha=.5)
+inducing_1d = dict(lw=0, s=500, facecolors=Tango.colorsHex['darkRed'])
+inducing_2d = dict(s=14, edgecolors='k', linewidth=.4, facecolors='white', alpha=.5)
 xerrorbar = dict(ecolor='k', fmt='none', elinewidth=.5, alpha=.5)
-yerrorbar = dict(ecolor=Tango.colorsHex['darkBlue'], fmt='none', elinewidth=.5, alpha=.5)
+yerrorbar = dict(ecolor=Tango.colorsHex['darkRed'], fmt='none', elinewidth=.5, alpha=.5)
 
 # GP plots
 meanplot_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=2)
 meanplot_2d = dict(cmap='hot', linewidth=.5)
-confidence_interval = dict(linecolor=Tango.colorsHex['darkBlue'],fillcolor=Tango.colorsHex['lightBlue'])
\ No newline at end of file
+confidence_interval = dict(edgecolor=Tango.colorsHex['darkBlue'],linewidth=.5,facecolor=Tango.colorsHex['lightBlue'],alpha=.3)
+density = dict(facecolor=Tango.colorsHex['mediumBlue'],edgecolors='none')
\ No newline at end of file
diff --git a/GPy/plotting/matplot_dep/plot_definitions.py b/GPy/plotting/matplot_dep/plot_definitions.py
index e256ef38..7ccf302a 100644
--- a/GPy/plotting/matplot_dep/plot_definitions.py
+++ b/GPy/plotting/matplot_dep/plot_definitions.py
@@ -59,14 +59,24 @@ class MatplotlibPlots(AbstractPlottingLibrary):
             ax.figure.tight_layout()
         except:
             pass
-        return ax
+        return plots
     
     def scatter(self, ax, X, Y, **kwargs):
         return ax.scatter(X, Y, **kwargs)
     
     def plot(self, ax, X, Y, **kwargs):
         return ax.plot(X, Y, **kwargs)
-    
+
+    def plot_axis_lines(self, ax, X, **kwargs):
+        from matplotlib import transforms
+        from matplotlib.path import Path
+        if 'transform' not in kwargs:
+            kwargs['transform'] = transforms.blended_transform_factory(ax.transData, ax.transAxes)
+        if 'marker' not in kwargs:
+            kwargs['marker'] = Path([[-.2,0.],    [-.2,.5],    [0.,1.],    [.2,.5],     [.2,0.],     [-.2,0.]],
+                                    [Path.MOVETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.LINETO, Path.CLOSEPOLY])
+        return ax.scatter(X, np.zeros_like(X), **kwargs)
+
     def xerrorbar(self, ax, X, Y, error, **kwargs):
         if not('linestyle' in kwargs or 'ls' in kwargs):
             kwargs['ls'] = 'none' 
@@ -74,7 +84,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
     
     def yerrorbar(self, ax, X, Y, error, **kwargs):
         if not('linestyle' in kwargs or 'ls' in kwargs):
-            kwargs['ls'] = 'none' 
+            kwargs['ls'] = 'none'
         return ax.errorbar(X, Y, yerr=error, **kwargs)
     
     def imshow(self, ax, X, **kwargs):
@@ -84,7 +94,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
         return ax.contour(X, Y, C, levels=np.linspace(C.min(), C.max(), levels), **kwargs)
 
     def fill_between(self, ax, X, lower, upper, **kwargs):
-        return ax.fill_between(X.flatten(), lower.flatten(), upper.flatten(), **kwargs)
+        return ax.fill_between(X, lower, upper, **kwargs)
 
     def fill_gradient(self, canvas, X, percentiles, **kwargs):
         ax = canvas
@@ -141,7 +151,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
                     continue
         
                 N = len(xslice)
-                X = np.zeros((2 * N + 2, 2), np.float)
+                p = np.zeros((2 * N + 2, 2), np.float)
         
                 # the purpose of the next two lines is for when y2 is a
                 # scalar like 0 and we want the fill to go all the way
@@ -149,15 +159,15 @@ class MatplotlibPlots(AbstractPlottingLibrary):
                 start = xslice[0], y2slice[0]
                 end = xslice[-1], y2slice[-1]
         
-                X[0] = start
-                X[N + 1] = end
+                p[0] = start
+                p[N + 1] = end
         
-                X[1:N + 1, 0] = xslice
-                X[1:N + 1, 1] = y1slice
-                X[N + 2:, 0] = xslice[::-1]
-                X[N + 2:, 1] = y2slice[::-1]
+                p[1:N + 1, 0] = xslice
+                p[1:N + 1, 1] = y1slice
+                p[N + 2:, 0] = xslice[::-1]
+                p[N + 2:, 1] = y2slice[::-1]
         
-                polys.append(X)
+                polys.append(p)
             polycol.extend(polys)
         from matplotlib.collections import PolyCollection
         plots.append(PolyCollection(polycol, **kwargs))