diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 2ed0f280..2eec11ee 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -595,6 +595,45 @@ class GP(Model):
                                      data_symbol=data_symbol, predict_kw=predict_kw,
                                      plot_training_data=plot_training_data, samples_y=samples_y, apply_link=apply_link, **kw)
 
+    def plot_density(self, levels=20, plot_limits=None, fignum=None, ax=None, 
+        fixed_inputs=[], plot_raw=False, edgecolor='none', facecolor='#3465a4', 
+        predict_kw=None,Y_metadata=None, 
+        apply_link=False, resolution=200, **patch_kw):
+        """
+        Plot the posterior density of the GP.
+          - In one dimension, the function is plotted with a shaded gradient, visualizing the density of the posterior.
+          - Only implemented for one dimension, for higher dimensions use `plot`.
+
+        :param levels: number of levels to plot in the density plot. This is a number between 1 and 100. 1 corresponds to the normal plot_fit.
+        :type levels: int
+        :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 index i should be set to value v.
+        :type fixed_inputs: a list of tuples
+        :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
+        :type resolution: int
+        :param edgecolor: color of line to plot [Tango.colorsHex['darkBlue']]
+        :type edgecolor: color either as Tango.colorsHex object or character ('r' is red, 'g' is green) as is standard in matplotlib
+        :param facecolor: color of fill [Tango.colorsHex['lightBlue']]
+        :type facecolor: color either as Tango.colorsHex object or character ('r' is red, 'g' is green) as is standard in matplotlib
+        :param Y_metadata: additional data associated with Y which may be needed
+        :type Y_metadata: dict
+        :param apply_link: if there is a link function of the likelihood, plot the link(f*) rather than f*, when plotting posterior samples f
+        :type apply_link: boolean
+        :param resolution: resolution of interpolation (how many points to interpolate of the posterior).
+        :type resolution: int
+        :param: patch_kw: the keyword arguments for the patchcollection fill.
+        """
+        assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
+        from ..plotting.matplot_dep import models_plots
+        return models_plots.plot_density(self, levels, plot_limits, fignum, ax,
+                                         fixed_inputs, plot_raw=plot_raw, 
+                                         Y_metadata=Y_metadata,
+                                         predict_kw=predict_kw,
+                                         apply_link=apply_link, 
+                                         edgecolor=edgecolor, facecolor=facecolor,
+                                         **patch_kw)
+        
 
     def plot_data(self, which_data_rows='all',
         which_data_ycols='all', visible_dims=None,
diff --git a/GPy/plotting/matplot_dep/base_plots.py b/GPy/plotting/matplot_dep/base_plots.py
index 5e513ec2..3746f79a 100644
--- a/GPy/plotting/matplot_dep/base_plots.py
+++ b/GPy/plotting/matplot_dep/base_plots.py
@@ -40,6 +40,88 @@ def gpplot(x, mu, lower, upper, edgecol='#3300FF', fillcol='#33CCFF', ax=None, f
 
     return plots
 
+def plot_gradient_fill(ax, x, percentiles, **kwargs):
+    plots = []
+
+    #here's the box
+    if 'linewidth' not in kwargs:
+        kwargs['linewidth'] = 0.5
+    if not 'alpha' in kwargs.keys():
+        kwargs['alpha'] = 1./(len(percentiles))
+    
+    # pop where from kwargs
+    where = kwargs.pop('where') if 'where' in kwargs else None
+    # pop interpolate, which we actually do not do here!
+    if 'interpolate' in kwargs: kwargs.pop('interpolate')
+    
+    def pairwise(inlist):
+        l = len(inlist)
+        for i in range(int(np.ceil(l/2.))):
+            yield inlist[:][i], inlist[:][(l-1)-i]
+    
+    polycol = []
+    for y1, y2 in pairwise(percentiles):
+        import matplotlib.mlab as mlab
+        # Handle united data, such as dates
+        ax._process_unit_info(xdata=x, ydata=y1)
+        ax._process_unit_info(ydata=y2)
+    
+        # Convert the arrays so we can work with them
+        from numpy import ma
+        x = ma.masked_invalid(ax.convert_xunits(x))
+        y1 = ma.masked_invalid(ax.convert_yunits(y1))
+        y2 = ma.masked_invalid(ax.convert_yunits(y2))
+    
+        if y1.ndim == 0:
+            y1 = np.ones_like(x) * y1
+        if y2.ndim == 0:
+            y2 = np.ones_like(x) * y2
+    
+        if where is None:
+            where = np.ones(len(x), np.bool)
+        else:
+            where = np.asarray(where, np.bool)
+    
+        if not (x.shape == y1.shape == y2.shape == where.shape):
+            raise ValueError("Argument dimensions are incompatible")
+    
+        mask = reduce(ma.mask_or, [ma.getmask(a) for a in (x, y1, y2)])
+        if mask is not ma.nomask:
+            where &= ~mask
+        
+        polys = []
+        for ind0, ind1 in mlab.contiguous_regions(where):
+            xslice = x[ind0:ind1]
+            y1slice = y1[ind0:ind1]
+            y2slice = y2[ind0:ind1]
+    
+            if not len(xslice):
+                continue
+    
+            N = len(xslice)
+            X = 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
+            # down to 0 even if none of the y1 sample points do
+            start = xslice[0], y2slice[0]
+            end = xslice[-1], y2slice[-1]
+    
+            X[0] = start
+            X[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]
+    
+            polys.append(X)
+        polycol.extend(polys)
+    from matplotlib.collections import PolyCollection
+    plots.append(PolyCollection(polycol, **kwargs))
+    ax.add_collection(plots[-1], autolim=True)
+    ax.autoscale_view()
+    return plots
 
 def gperrors(x, mu, lower, upper, edgecol=None, ax=None, fignum=None, **kwargs):
     _, axes = ax_default(fignum, ax)
diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py
index 208267c0..64d2d8fb 100644
--- a/GPy/plotting/matplot_dep/models_plots.py
+++ b/GPy/plotting/matplot_dep/models_plots.py
@@ -9,6 +9,7 @@ from ...models.sparse_gp_coregionalized_regression import SparseGPCoregionalized
 from scipy import sparse
 from ...core.parameterization.variational import VariationalPosterior
 from matplotlib import pyplot as plt
+from GPy.plotting.matplot_dep.base_plots import plot_gradient_fill
 
 
 def plot_data(model, which_data_rows='all',
@@ -299,6 +300,70 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
         raise NotImplementedError("Cannot define a frame with more than two input dimensions")
     return plots
 
+def plot_density(model, levels=20, plot_limits=None, fignum=None, ax=None, 
+        fixed_inputs=[], plot_raw=False, edgecolor='none', facecolor='#3465a4',  
+        predict_kw=None,Y_metadata=None, 
+        apply_link=False, resolution=200, **patch_kwargs):
+    #deal with optional arguments
+    if ax is None:
+        fig = plt.figure(num=fignum)
+        ax = fig.add_subplot(111)
+
+    if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs():
+        X = model.X.mean
+    else:
+        X = model.X
+    Y = model.Y
+    if sparse.issparse(Y): Y = Y.todense().view(np.ndarray)
+
+    if predict_kw is None:
+        predict_kw = {}
+
+    #work out what the inputs are for plotting (1D or 2D)
+    fixed_dims = np.array([i for i,v in fixed_inputs])
+    free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
+    plots = {}
+    #one dimensional plotting
+    if len(free_dims) == 1:
+        #define the frame on which to plot
+        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits, resolution=resolution)
+        Xgrid = np.empty((Xnew.shape[0],model.input_dim))
+        Xgrid[:,free_dims] = Xnew
+        for i,v in fixed_inputs:
+            Xgrid[:,i] = v
+
+        percs = np.linspace(2.5, 97.5, levels*2)
+
+        #make a prediction on the frame and plot it
+        if plot_raw:
+            raise NotImplementedError('What to do? What to do?')
+            #===================================================================
+            # m, v = model._raw_predict(Xgrid, **predict_kw)
+            # percentiles = model.predict_quantiles(Xgrid, )
+            # if apply_link:
+            #     lower = model.likelihood.gp_link.transf(m - 2*np.sqrt(v))
+            #     upper = model.likelihood.gp_link.transf(m + 2*np.sqrt(v))
+            #     #Once transformed this is now the median of the function
+            #     m = model.likelihood.gp_link.transf(m)
+            # else:
+            #     lower = m - 2*np.sqrt(v)
+            #     upper = m + 2*np.sqrt(v)
+            #===================================================================
+        else:
+            if isinstance(model,GPCoregionalizedRegression) or isinstance(model,SparseGPCoregionalizedRegression):
+                extra_data = Xgrid[:,-1:].astype(np.int)
+                if Y_metadata is None:
+                    Y_metadata = {'output_index': extra_data}
+                else:
+                    Y_metadata['output_index'] = extra_data
+            percentiles = [i[:, 0] for i in model.predict_quantiles(Xgrid, percs, Y_metadata=Y_metadata, **predict_kw)]
+        patch_kwargs['facecolor'] = facecolor
+        patch_kwargs['edgecolor'] = edgecolor
+        plots['density'] = plot_gradient_fill(ax, Xgrid[:, 0], percentiles, **patch_kwargs)
+    else:
+        raise NotImplementedError('Only 1D density plottable.')
+    return plots
+
 def plot_fit_f(model, *args, **kwargs):
     """
     Plot the GP's view of the world, where the data is normalized and before applying a likelihood.