From 46eca3bbdd99d654411c095e656b20cbf42df3d6 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 7 Oct 2013 12:41:20 +0100 Subject: [PATCH] Plots tidied up. --- GPy/core/gp_base.py | 411 +++++++++++++++++++++++++++--------------- GPy/core/sparse_gp.py | 163 +++++++++++++---- 2 files changed, 391 insertions(+), 183 deletions(-) diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index bd0b877e..083f9980 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -3,13 +3,14 @@ from .. import kern from ..util.plot import gpplot, Tango, x_frame1D, x_frame2D import pylab as pb from GPy.core.model import Model +import warnings +from ..likelihoods import Gaussian, Gaussian_Mixed_Noise class GPBase(Model): """ Gaussian process base model for holding shared behaviour between sparse_GP and GP models. """ - def __init__(self, X, likelihood, kernel, normalize_X=False): self.X = X assert len(self.X.shape) == 2 @@ -57,7 +58,59 @@ class GPBase(Model): self.X = state.pop() Model.setstate(self, state) - def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None,output=None): + def posterior_samples_f(self,X,size=10,which_parts='all',full_cov=True): + """ + Samples the posterior GP at the points X. + + :param X: The points at which to take the samples. + :type X: np.ndarray, Nnew x self.input_dim. + :param size: the number of a posteriori samples to plot. + :type size: int. + :param which_parts: which of the kernel functions to plot (additively). + :type which_parts: 'all', or list of bools. + :param full_cov: whether to return the full covariance matrix, or just the diagonal. + :type full_cov: bool. + :returns: Ysim: set of simulations, a Numpy array (N x samples). + """ + m, v = self._raw_predict(X, which_parts=which_parts, full_cov=full_cov) + v = v.reshape(m.size,-1) if len(v.shape)==3 else v + if not full_cov: + Ysim = np.random.multivariate_normal(m.flatten(), np.diag(v.flatten()), size).T + else: + Ysim = np.random.multivariate_normal(m.flatten(), v, size).T + + return Ysim + + def posterior_samples(self,X,size=10,which_parts='all',full_cov=True,noise_model=None): + """ + Samples the posterior GP at the points X. + + :param X: the points at which to take the samples. + :type X: np.ndarray, Nnew x self.input_dim. + :param size: the number of a posteriori samples to plot. + :type size: int. + :param which_parts: which of the kernel functions to plot (additively). + :type which_parts: 'all', or list of bools. + :param full_cov: whether to return the full covariance matrix, or just the diagonal. + :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). + """ + Ysim = self.posterior_samples_f(X, size, which_parts=which_parts, full_cov=full_cov) + if isinstance(self.likelihood,Gaussian): + noise_std = np.sqrt(self.likelihood._get_params()) + Ysim += np.random.normal(0,noise_std,Ysim.shape) + elif isinstance(self.likelihood,Gaussian_Mixed_Noise): + assert noise_model is not None, "A noise model must be specified." + noise_std = np.sqrt(self.likelihood._get_params()[noise_model]) + Ysim += np.random.normal(0,noise_std,Ysim.shape) + else: + Ysim = self.likelihood.noise_model.samples(Ysim) + + return Ysim + + def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None): """ Plot the GP's view of the world, where the data is normalized and the - In one dimension, the function is plotted with a shaded region identifying two standard deviations. @@ -89,82 +142,41 @@ class GPBase(Model): fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - if not hasattr(self,'multioutput'): + if self.X.shape[1] == 1: + resolution = resolution or 200 + Xnew, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits) - if self.X.shape[1] == 1: - Xnew, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits) - if samples == 0: - m, v = self._raw_predict(Xnew, which_parts=which_parts) - gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v), axes=ax) - ax.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5) - else: - m, v = self._raw_predict(Xnew, which_parts=which_parts, full_cov=True) - v = v.reshape(m.size,-1) if len(v.shape)==3 else v - Ysim = np.random.multivariate_normal(m.flatten(), v, samples) - gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None, ], axes=ax) - for i in range(samples): - ax.plot(Xnew, Ysim[i, :], Tango.colorsHex['darkBlue'], linewidth=0.25) + m, v = self._raw_predict(Xnew, which_parts=which_parts) + if samples: + Ysim = self.posterior_samples_f(Xnew, samples, which_parts=which_parts, full_cov=True) + for yi in Ysim.T: + ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v), axes=ax) - ax.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5) - ax.set_xlim(xmin, xmax) - ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None]))) - ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) - ax.set_ylim(ymin, ymax) + ax.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5) + ax.set_xlim(xmin, xmax) + ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None]))) + ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) + ax.set_ylim(ymin, ymax) - if hasattr(self,'Z'): - Zu = self.Z * self._Xscale + self._Xoffset - ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) + elif self.X.shape[1] == 2: - elif self.X.shape[1] == 2: - resolution = resolution or 50 - Xnew, xmin, xmax, xx, yy = x_frame2D(self.X, plot_limits, resolution) - m, v = self._raw_predict(Xnew, which_parts=which_parts) - m = m.reshape(resolution, resolution).T - ax.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable - ax.scatter(self.X[:, 0], self.X[:, 1], 40, self.likelihood.Y, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max()) # @UndefinedVariable - ax.set_xlim(xmin[0], xmax[0]) - ax.set_ylim(xmin[1], xmax[1]) + resolution = resolution or 50 + Xnew, xmin, xmax, xx, yy = x_frame2D(self.X, plot_limits, resolution) + m, v = self._raw_predict(Xnew, which_parts=which_parts) + m = m.reshape(resolution, resolution).T + ax.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable + ax.scatter(self.X[:, 0], self.X[:, 1], 40, self.likelihood.Y, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max()) # @UndefinedVariable + ax.set_xlim(xmin[0], xmax[0]) + ax.set_ylim(xmin[1], xmax[1]) + + if samples: + warnings.warn("Samples only implemented for 1 dimensional inputs.") - else: - raise NotImplementedError, "Cannot define a frame with more than two input dimensions" else: - assert len(self.likelihood.noise_model_list) > output, 'The model has only %s outputs.' %self.num_outputs + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - if self.X.shape[1] == 2: - Xu = self.X[self.X[:,-1]==output ,0:1] - Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) - - if samples == 0: - m, v = self._raw_predict_single_output(Xnew, output=output, which_parts=which_parts) - gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v), axes=ax) - ax.plot(Xu[which_data], self.likelihood.Y[self.likelihood.index==output][:,None], 'kx', mew=1.5) - else: - m, v = self._raw_predict_single_output(Xnew, output=output, which_parts=which_parts, full_cov=True) - v = v.reshape(m.size,-1) if len(v.shape)==3 else v - Ysim = np.random.multivariate_normal(m.flatten(), v, samples) - gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None, ], axes=ax) - for i in range(samples): - ax.plot(Xnew, Ysim[i, :], Tango.colorsHex['darkBlue'], linewidth=0.25) - ax.set_xlim(xmin, xmax) - ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None]))) - ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) - ax.set_ylim(ymin, ymax) - - elif self.X.shape[1] == 3: - raise NotImplementedError, "Plots not implemented for multioutput models with 2D inputs...yet" - assert self.num_outputs >= output, 'The model has only %s outputs.' %self.num_outputs - - else: - raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - - if hasattr(self,'Z'): - Zu = self.Z[self.Z[:,-1]==output,:] - Zu = self.Z * self._Xscale + self._Xoffset - Zu = self.Z[self.Z[:,-1]==output ,0:1] #?? - ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) - - - def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, output=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): + def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): """ Plot the GP with noise where the likelihood is Gaussian. @@ -200,7 +212,6 @@ class GPBase(Model): :param fillcol: color of fill :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure """ - # TODO include samples if which_data == 'all': which_data = slice(None) @@ -208,98 +219,202 @@ class GPBase(Model): fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - if not hasattr(self,'multioutput'): + plotdims = self.input_dim - len(fixed_inputs) + if plotdims == 1: + resolution = resolution or 200 - plotdims = self.input_dim - len(fixed_inputs) - if plotdims == 1: - resolution = resolution or 200 + Xu = self.X * self._Xscale + self._Xoffset #NOTE self.X are the normalized values now - Xu = self.X * self._Xscale + self._Xoffset #NOTE self.X are the normalized values now + fixed_dims = np.array([i for i,v in fixed_inputs]) + freedim = np.setdiff1d(np.arange(self.input_dim),fixed_dims) - fixed_dims = np.array([i for i,v in fixed_inputs]) - freedim = np.setdiff1d(np.arange(self.input_dim),fixed_dims) + Xnew, xmin, xmax = x_frame1D(Xu[:,freedim], plot_limits=plot_limits) + Xgrid = np.empty((Xnew.shape[0],self.input_dim)) + Xgrid[:,freedim] = Xnew + for i,v in fixed_inputs: + Xgrid[:,i] = v - Xnew, xmin, xmax = x_frame1D(Xu[:,freedim], plot_limits=plot_limits) - Xgrid = np.empty((Xnew.shape[0],self.input_dim)) - Xgrid[:,freedim] = Xnew - for i,v in fixed_inputs: - Xgrid[:,i] = v + m, v, lower, upper = self.predict(Xgrid, which_parts=which_parts) - m, _, lower, upper = self.predict(Xgrid, which_parts=which_parts) - for d in range(m.shape[1]): - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) - ax.plot(Xu[which_data,freedim], self.likelihood.data[which_data, d], 'kx', mew=1.5) - ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) - ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) + if samples: #NOTE not tested with fixed_inputs + Ysim = self.posterior_samples(Xgrid, samples, which_parts=which_parts, full_cov=True) + for yi in Ysim.T: + ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. + for d in range(m.shape[1]): + gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) + ax.plot(Xu[which_data,freedim], self.likelihood.data[which_data, d], 'kx', mew=1.5) + ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) + ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + elif self.X.shape[1] == 2: - Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits,resolution=resolution) - m, _, lower, upper = self.predict(Xnew, which_parts=which_parts) - for d in range(m.shape[1]): - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax) - ax.plot(Xu[which_data], self.likelihood.data[which_data, d], 'kx', mew=1.5) - ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) - ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) + resolution = resolution or 50 + Xnew, _, _, xmin, xmax = x_frame2D(self.X, plot_limits, resolution) + x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution) + m, _, lower, upper = self.predict(Xnew, which_parts=which_parts) + m = m.reshape(resolution, resolution).T + ax.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable + Yf = self.likelihood.Y.flatten() + ax.scatter(self.X[:, 0], self.X[:, 1], 40, Yf, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) # @UndefinedVariable + ax.set_xlim(xmin[0], xmax[0]) + ax.set_ylim(xmin[1], xmax[1]) - elif self.X.shape[1] == 2: - resolution = resolution or 50 - Xnew, _, _, xmin, xmax = x_frame2D(self.X, plot_limits, resolution) - x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution) - m, _, lower, upper = self.predict(Xnew, which_parts=which_parts) - m = m.reshape(resolution, resolution).T - ax.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) # @UndefinedVariable - Yf = self.likelihood.Y.flatten() - ax.scatter(self.X[:, 0], self.X[:, 1], 40, Yf, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) # @UndefinedVariable - ax.set_xlim(xmin[0], xmax[0]) - ax.set_ylim(xmin[1], xmax[1]) - - else: - raise NotImplementedError, "Cannot define a frame with more than two input dimensions" + if samples: + warnings.warn("Samples only implemented for 1 dimensional inputs.") else: - assert len(self.likelihood.noise_model_list) > output, 'The model has only %s outputs.' %self.num_outputs - if self.X.shape[1] == 2: - resolution = resolution or 200 - Xu = self.X[self.X[:,-1]==output,:] #keep the output of interest - Xu = self.X * self._Xscale + self._Xoffset - Xu = self.X[self.X[:,-1]==output ,0:1] #get rid of the index column + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) - m, _, lower, upper = self.predict_single_output(Xnew, which_parts=which_parts,output=output) + def plot_single_output_f(self, output=None, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None): + """ + For a specific output, in a multioutput model, this function works just as plot_f on single output models. - for d in range(m.shape[1]): - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax) - ax.plot(Xu[which_data], self.likelihood.noise_model_list[output].data, 'kx', mew=1.5) - ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) - ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) - ax.set_xlim(xmin, xmax) - ax.set_ylim(ymin, ymax) + :param output: which output to plot (for multiple output models only) + :type output: integer (first output is 0) + :param samples: the number of a posteriori samples to plot + :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits + :param which_data: which if the training data to plot (default all) + :type which_data: 'all' or a slice object to slice self.X, self.Y + :param which_parts: which of the kernel functions to plot (additively) + :type which_parts: 'all', or list of bools + :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 full_cov: + :type full_cov: bool + :param fignum: figure to plot on. + :type fignum: figure number + :param ax: axes to plot on. + :type ax: axes handle + """ + assert output is not None, "An output must be specified." + assert len(self.likelihood.noise_model_list) > output, "The model has only %s outputs." %(self.output_dim + 1) - elif self.X.shape[1] == 3: - raise NotImplementedError, "Plots not yet implemented for multioutput models with 2D inputs" - resolution = resolution or 50 - - else: - raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - - """ - def samples_f(self,X,samples=10, which_data='all', which_parts='all',output=None): if which_data == 'all': which_data = slice(None) - if hasattr(self,'multioutput'): - np.hstack([X,np.ones((X.shape[0],1))*output]) + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) - m, v = self._raw_predict(X, which_parts=which_parts, full_cov=True) - v = v.reshape(m.size,-1) if len(v.shape)==3 else v - Ysim = np.random.multivariate_normal(m.flatten(), v, samples) - #gpplot(X, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None, ], axes=ax) - for i in range(samples): - ax.plot(X, Ysim[i, :], Tango.colorsHex['darkBlue'], linewidth=0.25) + if self.X.shape[1] == 2: + Xu = self.X[self.X[:,-1]==output ,0:1] + Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) + Xnew_indexed = self._add_output_index(Xnew,output) - """ + m, v = self._raw_predict(Xnew_indexed, which_parts=which_parts) + + if samples: + Ysim = self.posterior_samples_f(Xnew_indexed, samples, which_parts=which_parts, full_cov=True) + for yi in Ysim.T: + ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + + gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v), axes=ax) + ax.plot(Xu[which_data], self.likelihood.Y[self.likelihood.index==output][:,None], 'kx', mew=1.5) + ax.set_xlim(xmin, xmax) + ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None]))) + ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) + ax.set_ylim(ymin, ymax) + + elif self.X.shape[1] == 3: + raise NotImplementedError, "Plots not implemented for multioutput models with 2D inputs...yet" + #if samples: + # warnings.warn("Samples only implemented for 1 dimensional inputs.") + + else: + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" + + + def plot_single_output(self, output=None, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, fixed_inputs=[], linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue']): + """ + For a specific output, in a multioutput model, this function works just as plot_f on single output models. + + :param output: which output to plot (for multiple output models only) + :type output: integer (first output is 0) + :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 which_data: which if the training data to plot (default all) + :type which_data: 'all' or a slice object to slice self.X, self.Y + :param which_parts: which of the kernel functions to plot (additively) + :type which_parts: 'all', or list of bools + :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 levels: number of levels to plot in a contour plot. + :type levels: int + :param samples: the number of a posteriori samples to plot + :type samples: int + :param fignum: figure to plot on. + :type fignum: figure number + :param ax: axes to plot on. + :type ax: axes handle + :type output: integer (first output is 0) + :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 linecol: color of line to plot. + :type linecol: + :param fillcol: color of fill + :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure + """ + assert output is not None, "An output must be specified." + assert len(self.likelihood.noise_model_list) > output, "The model has only %s outputs." %(self.output_dim + 1) + if which_data == 'all': + which_data = slice(None) + + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) + + if self.X.shape[1] == 2: + resolution = resolution or 200 + + Xu = self.X[self.X[:,-1]==output,:] #keep the output of interest + Xu = self.X * self._Xscale + self._Xoffset + Xu = self.X[self.X[:,-1]==output ,0:1] #get rid of the index column + + Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) + Xnew_indexed = self._add_output_index(Xnew,output) + + + m, v, lower, upper = self.predict(Xnew_indexed, which_parts=which_parts,noise_model=output) + + if samples: #NOTE not tested with fixed_inputs + Ysim = self.posterior_samples(Xnew_indexed, samples, which_parts=which_parts, full_cov=True,noise_model=output) + for yi in Ysim.T: + ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + + for d in range(m.shape[1]): + gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) + ax.plot(Xu[which_data], self.likelihood.noise_model_list[output].data, 'kx', mew=1.5) + ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper)) + ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin) + ax.set_xlim(xmin, xmax) + ax.set_ylim(ymin, ymax) + + elif self.X.shape[1] == 3: + raise NotImplementedError, "Plots not implemented for multioutput models with 2D inputs...yet" + #if samples: + # warnings.warn("Samples only implemented for 1 dimensional inputs.") + + else: + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" + + + def _add_output_index(self,X,output): + """ + In a multioutput model, appends an index column to X to specify the output it is related to. + + :param X: Input data + :type X: np.ndarray, N x self.input_dim + :param output: output X is related to + :type output: integer in {0,..., output_dim-1} + + .. Note:: For multiple non-independent outputs models only. + """ + + assert hasattr(self,'multioutput'), 'This function is for multiple output models only.' + + index = np.ones((X.shape[0],1))*output + return np.hstack((X,index)) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 834bdd84..d4b33ed2 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -34,7 +34,6 @@ class SparseGP(GPBase): self.Z = Z self.num_inducing = Z.shape[0] -# self.likelihood = likelihood if X_variance is None: self.has_uncertain_inputs = False @@ -305,9 +304,8 @@ class SparseGP(GPBase): return mu, var[:, None] - def predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): + def predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False, **likelihood_args): """ - Predict the function(s) at the new point(s) Xnew. **Arguments** @@ -338,56 +336,90 @@ class SparseGP(GPBase): mu, var = self._raw_predict(Xnew, X_variance_new, full_cov=full_cov, which_parts=which_parts) # now push through likelihood - mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov) + mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, **likelihood_args) return mean, var, _025pm, _975pm - def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None, output=None): + + def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None): + """ + Plot the GP's view of the world, where the data is normalized and the + - In one dimension, the function is plotted with a shaded region identifying two standard deviations. + - In two dimsensions, a contour-plot shows the mean predicted function + - Not implemented in higher dimensions + + :param samples: the number of a posteriori samples to plot + :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits + :param which_data: which if the training data to plot (default all) + :type which_data: 'all' or a slice object to slice self.X, self.Y + :param which_parts: which of the kernel functions to plot (additively) + :type which_parts: 'all', or list of bools + :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 full_cov: + :type full_cov: bool + :param fignum: figure to plot on. + :type fignum: figure number + :param ax: axes to plot on. + :type ax: axes handle + + :param output: which output to plot (for multiple output models only) + :type output: integer (first output is 0) + """ if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) + if fignum is None and ax is None: + fignum = fig.num if which_data is 'all': which_data = slice(None) - GPBase.plot(self, samples=0, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=None, levels=20, ax=ax, output=output) + GPBase.plot_f(self, samples=samples, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=resolution, full_cov=full_cov, fignum=fignum, ax=ax) - if not hasattr(self,'multioutput'): + if self.X.shape[1] == 1: + if self.has_uncertain_inputs: + Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now + ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], + xerr=2 * np.sqrt(self.X_variance[which_data, 0]), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + Zu = self.Z * self._Xscale + self._Xoffset + ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) - if self.X.shape[1] == 1: - if self.has_uncertain_inputs: - Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now - ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], - xerr=2 * np.sqrt(self.X_variance[which_data, 0]), - ecolor='k', fmt=None, elinewidth=.5, alpha=.5) - Zu = self.Z * self._Xscale + self._Xoffset - ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) + elif self.X.shape[1] == 2: + Zu = self.Z * self._Xscale + self._Xoffset + ax.plot(Zu[:, 0], Zu[:, 1], 'wo') - elif self.X.shape[1] == 2: - Zu = self.Z * self._Xscale + self._Xoffset - ax.plot(Zu[:, 0], Zu[:, 1], 'wo') else: - if self.X.shape[1] == 2 and hasattr(self,'multioutput'): - """ - Xu = self.X[self.X[:,-1]==output,:] - if self.has_uncertain_inputs: - Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" - Xu = self.X[self.X[:,-1]==output ,0:1] #?? + def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None): + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) + if fignum is None and ax is None: + fignum = fig.num + if which_data is 'all': + which_data = slice(None) - ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], - xerr=2 * np.sqrt(self.X_variance[which_data, 0]), - ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + GPBase.plot(self, samples=samples, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=resolution, levels=20, fignum=fignum, ax=ax) - """ - Zu = self.Z[self.Z[:,-1]==output,:] - Zu = self.Z * self._Xscale + self._Xoffset - Zu = self.Z[self.Z[:,-1]==output ,0:1] #?? - ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) - #ax.set_ylim(ax.get_ylim()[0],) + if self.X.shape[1] == 1: + if self.has_uncertain_inputs: + Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now + ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], + xerr=2 * np.sqrt(self.X_variance[which_data, 0]), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + Zu = self.Z * self._Xscale + self._Xoffset + ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) - else: - raise NotImplementedError, "Cannot define a frame with more than two input dimensions" + elif self.X.shape[1] == 2: + Zu = self.Z * self._Xscale + self._Xoffset + ax.plot(Zu[:, 0], Zu[:, 1], 'wo') + + + else: + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" def predict_single_output(self, Xnew, output=0, which_parts='all', full_cov=False): """ @@ -470,3 +502,64 @@ class SparseGP(GPBase): var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) return mu, var[:, None] + + + def plot_single_output_f(self, output=None, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None): + + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) + if fignum is None and ax is None: + fignum = fig.num + if which_data is 'all': + which_data = slice(None) + + GPBase.plot_single_output_f(self, output=output, samples=samples, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=resolution, full_cov=full_cov, fignum=fignum, ax=ax) + + if self.X.shape[1] == 2: + if self.has_uncertain_inputs: + Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now + ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], + xerr=2 * np.sqrt(self.X_variance[which_data, 0]), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + Zu = self.Z * self._Xscale + self._Xoffset + Zu = Zu[Zu[:,1]==output,0:1] + ax.plot(Zu[:,0], np.zeros_like(Zu[:,0]) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) + + elif self.X.shape[1] == 2: + Zu = self.Z * self._Xscale + self._Xoffset + Zu = Zu[Zu[:,1]==output,0:2] + ax.plot(Zu[:, 0], Zu[:, 1], 'wo') + + + else: + raise NotImplementedError, "Cannot define a frame with more than two input dimensions" + + def plot_single_output(self, output=None, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None): + if ax is None: + fig = pb.figure(num=fignum) + ax = fig.add_subplot(111) + if fignum is None and ax is None: + fignum = fig.num + if which_data is 'all': + which_data = slice(None) + + GPBase.plot_single_output(self, samples=samples, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=resolution, levels=20, fignum=fignum, ax=ax, output=output) + + if self.X.shape[1] == 2: + if self.has_uncertain_inputs: + Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now + ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0], + xerr=2 * np.sqrt(self.X_variance[which_data, 0]), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + Zu = self.Z * self._Xscale + self._Xoffset + Zu = Zu[Zu[:,1]==output,0:1] + ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12) + + elif self.X.shape[1] == 3: + Zu = self.Z * self._Xscale + self._Xoffset + Zu = Zu[Zu[:,1]==output,0:1] + ax.plot(Zu[:, 0], Zu[:, 1], 'wo') + + else: + raise NotImplementedError, "Cannot define a frame with more than two input dimensions"