diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 2a60cb59..ef71ea2c 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -437,7 +437,7 @@ class GP(Model): mag[n] = np.sqrt(np.linalg.det(G[n, :, :])) return mag - def posterior_samples_f(self,X,size=10, full_cov=True): + def posterior_samples_f(self,X, size=10, full_cov=True, **predict_kwargs): """ Samples the posterior GP at the points X. @@ -448,20 +448,32 @@ class GP(Model): :param full_cov: whether to return the full covariance matrix, or just the diagonal. :type full_cov: bool. :returns: fsim: set of simulations - :rtype: np.ndarray (N x samples) + :rtype: np.ndarray (D x N x samples) (if D==1 we flatten out the first dimension) """ - m, v = self._raw_predict(X, full_cov=full_cov) + m, v = self._raw_predict(X, full_cov=full_cov, **predict_kwargs) if self.normalizer is not None: m, v = self.normalizer.inverse_mean(m), self.normalizer.inverse_variance(v) - v = v.reshape(m.size,-1) if len(v.shape)==3 else v - if not full_cov: - fsim = np.random.multivariate_normal(m.flatten(), np.diag(v.flatten()), size).T - else: - fsim = np.random.multivariate_normal(m.flatten(), v, size).T + + def sim_one_dim(m, v): + if not full_cov: + return np.random.multivariate_normal(m.flatten(), np.diag(v.flatten()), size).T + else: + return np.random.multivariate_normal(m.flatten(), v, size).T + if self.output_dim == 1: + return sim_one_dim(m, v) + else: + fsim = np.empty((self.output_dim, self.num_data, size)) + for d in range(self.output_dim): + if full_cov and v.ndim == 3: + fsim[d] = sim_one_dim(m[:, d], v[:, :, d]) + elif (not full_cov) and v.ndim == 2: + fsim[d] = sim_one_dim(m[:, d], v[:, d]) + else: + fsim[d] = sim_one_dim(m[:, d], v) return fsim - def posterior_samples(self, X, size=10, full_cov=False, Y_metadata=None): + def posterior_samples(self, X, size=10, full_cov=False, Y_metadata=None, likelihood=None, **predict_kwargs): """ Samples the posterior GP at the points X. @@ -473,11 +485,18 @@ class GP(Model): :type full_cov: bool. :param noise_model: for mixed noise likelihood, the noise model to use in the samples. :type noise_model: integer. - :returns: Ysim: set of simulations, a Numpy array (N x samples). + :returns: Ysim: set of simulations, + :rtype: np.ndarray (D x N x samples) (if D==1 we flatten out the first dimension) """ - fsim = self.posterior_samples_f(X, size, full_cov=full_cov) - Ysim = self.likelihood.samples(fsim, Y_metadata=Y_metadata) - return Ysim + fsim = self.posterior_samples_f(X, size, full_cov=full_cov, **predict_kwargs) + if likelihood is None: + likelihood = self.likelihood + if fsim.ndim == 3: + for d in range(fsim.shape[0]): + fsim[d] = likelihood.samples(fsim[d], Y_metadata=Y_metadata) + else: + fsim = likelihood.samples(fsim, Y_metadata=Y_metadata) + return fsim def input_sensitivity(self, summarize=True): """ diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py index 23bbbe6b..e3700b18 100644 --- a/GPy/plotting/__init__.py +++ b/GPy/plotting/__init__.py @@ -33,6 +33,7 @@ if config.get('plotting', 'library') is not 'none': 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 + GP.plot_samples = gpy_plot.gp_plots.plot_samples from ..core import SparseGP SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing diff --git a/GPy/plotting/gpy_plot/data_plots.py b/GPy/plotting/gpy_plot/data_plots.py index 1a94b4a4..f63aa742 100644 --- a/GPy/plotting/gpy_plot/data_plots.py +++ b/GPy/plotting/gpy_plot/data_plots.py @@ -179,8 +179,12 @@ def _plot_errorbars_trainset(self, canvas, 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), + if predict_kw is None: + predict_kw = {} + if 'Y_metadata' not in predict_kw: + predict_kw['Y_metadata'] = self.Y_metadata or {} + _, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, + apply_link, (0, 100), ycols, predict_kw) for d in ycols: plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], Y[rows,d], diff --git a/GPy/plotting/gpy_plot/gp_plots.py b/GPy/plotting/gpy_plot/gp_plots.py index f43ee265..7183834b 100644 --- a/GPy/plotting/gpy_plot/gp_plots.py +++ b/GPy/plotting/gpy_plot/gp_plots.py @@ -37,7 +37,7 @@ from .plot_util import helper_for_plot_data, update_not_existing_kwargs, \ def plot_mean(self, plot_limits=None, fixed_inputs=None, resolution=None, plot_raw=False, - Y_metadata=None, apply_link=False, + apply_link=False, which_data_ycols='all', levels=20, predict_kw=None, @@ -62,22 +62,21 @@ def plot_mean(self, plot_limits=None, fixed_inputs=None, """ canvas, kwargs = pl.get_new_canvas(kwargs) plots = _plot_mean(self, canvas, plot_limits, fixed_inputs, - resolution, plot_raw, Y_metadata, + resolution, plot_raw, apply_link, which_data_ycols, levels, predict_kw, **kwargs) return pl.show_canvas(canvas, plots) -@wraps(plot_mean) def _plot_mean(self, canvas, plot_limits=None, fixed_inputs=None, resolution=None, plot_raw=False, - Y_metadata=None, apply_link=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) if len(free_dims)<=2: - mu, _ = helper_predict_with_model(self, Xgrid, plot_raw, + mu, _, _ = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, None, get_which_data_ycols(self, which_data_ycols), predict_kw) @@ -141,7 +140,7 @@ def _plot_confidence(self, canvas, lower, upper, plot_limits=None, fixed_inputs= if len(free_dims)<=1: if len(free_dims)==1: - _, percs = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, + _, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, (lower, upper), ycols, predict_kw) @@ -155,11 +154,63 @@ def _plot_confidence(self, canvas, lower, upper, plot_limits=None, fixed_inputs= raise RuntimeError('Can only plot confidence interval in one input dimension') +def plot_samples(self, plot_limits=None, fixed_inputs=None, + resolution=None, plot_raw=True, + apply_link=False, + which_data_ycols='all', + samples=3, predict_kw=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 + :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 [defaults are 1D:200, 2D:50] + :param bool plot_raw: plot the latent function (usually denoted f) only? This is usually what you want! + :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 dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=) in here + :param int levels: for 2D plotting, the number of contour levels to use is + """ + canvas, kwargs = pl.get_new_canvas(kwargs) + plots = _plot_samples(self, canvas, plot_limits, fixed_inputs, + resolution, plot_raw, + apply_link, which_data_ycols, samples, + predict_kw, **kwargs) + return pl.show_canvas(canvas, plots) + +def _plot_samples(self, canvas, plot_limits=None, fixed_inputs=None, + resolution=None, plot_raw=False, + apply_link=False, + which_data_ycols=None, + samples=3, + predict_kw=None, **kwargs): + _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, fixed_inputs, resolution) + + if len(free_dims)<2: + + if len(free_dims)==1: + # 1D plotting: + _, _, samples = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, + None, get_which_data_ycols(self, which_data_ycols), predict_kw, samples) + update_not_existing_kwargs(kwargs, pl.defaults.samples_1d) + return dict(gpmean=[pl.plot(canvas, Xgrid[:, free_dims], samples, **kwargs)]) + else: + pass # Nothing to plot! + else: + raise RuntimeError('Cannot plot mean in more then 1 input dimensions') + + def plot_density(self, plot_limits=None, fixed_inputs=None, resolution=None, plot_raw=False, apply_link=False, which_data_ycols='all', - levels=20, + levels=35, predict_kw=None, **kwargs): """ @@ -178,7 +229,7 @@ def plot_density(self, plot_limits=None, fixed_inputs=None, :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 int levels: the number of levels in the density (number bigger then 1, where 35 is smooth and 1 is the same as plot_confidence). You can go higher then 50 if the result is not smooth enough for you. :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=) in here """ canvas, kwargs = pl.get_new_canvas(kwargs) @@ -193,7 +244,7 @@ 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, + levels=35, predict_kw=None, **kwargs): _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, fixed_inputs, resolution) @@ -203,7 +254,7 @@ def _plot_density(self, canvas, plot_limits=None, fixed_inputs=None, if len(free_dims)<=1: if len(free_dims)==1: - _, percs = helper_predict_with_model(self, Xgrid, plot_raw, + _, 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) diff --git a/GPy/plotting/gpy_plot/plot_util.py b/GPy/plotting/gpy_plot/plot_util.py index 04e4fa0e..8177ea63 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, percentiles, which_data_ycols, predict_kw): +def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, which_data_ycols, predict_kw, samples=0): """ Make the right decisions for prediction with a model based on the standard arguments of plotting. @@ -46,12 +46,12 @@ def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, wh if plot_raw: from ...likelihoods import Gaussian from ...likelihoods.link_functions import Identity - lik = Gaussian(Identity(), 0) # Make the likelihood not add any noise + lik = Gaussian(Identity(), 1e-9) # 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 {} + predict_kw['Y_metadata'] = {} if 'output_index' not in predict_kw['Y_metadata']: predict_kw['Y_metadata']['output_index'] = Xgrid[:,-1:].astype(np.int) @@ -61,6 +61,12 @@ def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, wh percentiles = self.predict_quantiles(Xgrid, quantiles=percentiles, **predict_kw) else: percentiles = [] + if samples > 0: + fsamples = self.posterior_samples(Xgrid, full_cov=True, size=samples, **predict_kw) + fsamples = fsamples[which_data_ycols] if fsamples.ndim == 3 else fsamples + else: + fsamples = None + # Filter out the ycolums which we want to plot: retmu = mu[:, which_data_ycols] percs = [p[:, which_data_ycols] for p in percentiles] @@ -70,8 +76,13 @@ def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, wh retmu[:, [i]] = self.likelihood.gp_link.transf(mu[:, [i]]) for perc in percs: perc[:, [i]] = self.likelihood.gp_link.transf(perc[:, [i]]) - - return retmu, percs + if fsamples is not None and fsamples.ndim == 3: + for s in range(fsamples.shape[-1]): + fsamples[i, :, s] = self.likelihood.gp_link.transf(fsamples[i, :, s]) + elif fsamples is not None: + for s in range(fsamples.shape[-1]): + fsamples[:, s] = self.likelihood.gp_link.transf(fsamples[:, s]) + return retmu, percs, fsamples def helper_for_plot_data(self, plot_limits, fixed_inputs, resolution): """ diff --git a/GPy/plotting/matplot_dep/defaults.py b/GPy/plotting/matplot_dep/defaults.py index a122783b..5effddec 100644 --- a/GPy/plotting/matplot_dep/defaults.py +++ b/GPy/plotting/matplot_dep/defaults.py @@ -52,5 +52,6 @@ yerrorbar = dict(ecolor=Tango.colorsHex['darkRed'], fmt='none', elinewidth=.5, a # GP plots meanplot_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=2) meanplot_2d = dict(cmap='hot', linewidth=.5) +samples_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=.3) 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 7ccf302a..9b510a97 100644 --- a/GPy/plotting/matplot_dep/plot_definitions.py +++ b/GPy/plotting/matplot_dep/plot_definitions.py @@ -100,7 +100,7 @@ class MatplotlibPlots(AbstractPlottingLibrary): ax = canvas plots = [] if not 'alpha' in kwargs.keys(): - kwargs['alpha'] = 1./(len(percentiles)) + kwargs['alpha'] = 3./max(4, (len(percentiles))) # pop where from kwargs where = kwargs.pop('where') if 'where' in kwargs else None