From ef778f75d451e6339abbef5a58717c9c0b578a8c Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Thu, 10 Sep 2015 11:08:28 +0100 Subject: [PATCH 1/2] Updated sampling and plots to be correct shape, and changed plotting of sampling to be posterior samples p(f*|f), like it used to be, and samples_y to plot samples of p(y*|y) --- GPy/core/gp.py | 12 ++++--- GPy/likelihoods/exponential.py | 2 +- GPy/likelihoods/likelihood.py | 2 +- GPy/likelihoods/poisson.py | 8 +++-- GPy/plotting/matplot_dep/models_plots.py | 46 +++++++++++++----------- 5 files changed, 41 insertions(+), 29 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 9a199faa..903044b9 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -535,7 +535,7 @@ class GP(Model): which_data_ycols='all', fixed_inputs=[], levels=20, samples=0, fignum=None, ax=None, resolution=None, plot_raw=False, linecol=None,fillcol=None, Y_metadata=None, - data_symbol='kx', predict_kw=None, plot_training_data=True): + data_symbol='kx', predict_kw=None, plot_training_data=True, samples_y=0, apply_link=False): """ Plot the posterior of the GP. - In one dimension, the function is plotted with a shaded region identifying two standard deviations. @@ -558,7 +558,7 @@ class GP(Model): :param levels: number of levels to plot in a contour plot. :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure :type levels: int - :param samples: the number of a posteriori samples to plot + :param samples: the number of a posteriori samples to plot, p(f*|y) :type samples: int :param fignum: figure to plot on. :type fignum: figure number @@ -574,6 +574,10 @@ class GP(Model): :type data_symbol: color either as Tango.colorsHex object or character ('r' is red, 'g' is green) alongside marker type, as is standard in matplotlib. :param plot_training_data: whether or not to plot the training points :type plot_training_data: boolean + :param samples_y: the number of a posteriori samples to plot, p(y*|y) + :type samples_y: int + :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 """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import models_plots @@ -587,7 +591,7 @@ class GP(Model): levels, samples, fignum, ax, resolution, plot_raw=plot_raw, Y_metadata=Y_metadata, data_symbol=data_symbol, predict_kw=predict_kw, - plot_training_data=plot_training_data, **kw) + plot_training_data=plot_training_data, samples_y=samples_y, apply_link=apply_link, **kw) def plot_data(self, which_data_rows='all', @@ -613,7 +617,7 @@ class GP(Model): :param levels: number of levels to plot in a contour plot. :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure :type levels: int - :param samples: the number of a posteriori samples to plot + :param samples: the number of a posteriori samples to plot, p(f*|y) :type samples: int :param fignum: figure to plot on. :type fignum: figure number diff --git a/GPy/likelihoods/exponential.py b/GPy/likelihoods/exponential.py index 0a6c543d..ecf0977e 100644 --- a/GPy/likelihoods/exponential.py +++ b/GPy/likelihoods/exponential.py @@ -124,7 +124,7 @@ class Exponential(Likelihood): #d3lik_dlink3 = 6*y/(link_f**4) - 2./(link_f**3) return d3lik_dlink3 - def samples(self, gp): + def samples(self, gp, Y_metadata=None): """ Returns a set of samples of observations based on a given value of the latent variable. diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index e961dd1e..74c4c6fd 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -622,7 +622,7 @@ class Likelihood(Parameterized): Nf_samp = 300 Ny_samp = 1 s = np.random.randn(mu.shape[0], Nf_samp)*np.sqrt(var) + mu - ss_y = self.samples(s, Y_metadata, samples=Ny_samp) + ss_y = self.samples(s, Y_metadata)#, samples=Ny_samp) #ss_y = ss_y.reshape(mu.shape[0], mu.shape[1], Nf_samp*Ny_samp) pred_quantiles = [np.percentile(ss_y, q, axis=1)[:,None] for q in quantiles] diff --git a/GPy/likelihoods/poisson.py b/GPy/likelihoods/poisson.py index 1c3fec9e..cfe279bb 100644 --- a/GPy/likelihoods/poisson.py +++ b/GPy/likelihoods/poisson.py @@ -137,7 +137,7 @@ class Poisson(Likelihood): """ return self.gp_link.transf(gp) - def samples(self, gp, Y_metadata=None, samples=1): + def samples(self, gp, Y_metadata=None): """ Returns a set of samples of observations based on a given value of the latent variable. @@ -145,5 +145,7 @@ class Poisson(Likelihood): """ orig_shape = gp.shape gp = gp.flatten() - Ysim = np.random.poisson(self.gp_link.transf(gp), [samples, gp.size]).T - return Ysim.reshape(orig_shape+(samples,)) + # Ysim = np.random.poisson(self.gp_link.transf(gp), [samples, gp.size]).T + # return Ysim.reshape(orig_shape+(samples,)) + Ysim = np.random.poisson(self.gp_link.transf(gp)) + return Ysim.reshape(orig_shape) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 3a5a01d2..a15146cf 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -75,7 +75,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', levels=20, samples=0, fignum=None, ax=None, resolution=None, plot_raw=False, linecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'], Y_metadata=None, data_symbol='kx', - apply_link=False, samples_f=0, plot_uncertain_inputs=True, predict_kw=None, plot_training_data=True): + apply_link=False, samples_y=0, plot_uncertain_inputs=True, predict_kw=None, plot_training_data=True): """ Plot the posterior of the GP. - In one dimension, the function is plotted with a shaded region identifying two standard deviations. @@ -93,24 +93,30 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', :type which_data_rows: 'all' or a list of integers :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 levels: number of levels to plot in a contour plot. + :param levels: for 2D plotting, the number of contour levels to use is ax is None, create a new figure :type levels: int - :param samples: the number of a posteriori samples to plot p(y*|y) + :param samples: the number of a posteriori samples to plot p(f*|y) :type samples: int :param fignum: figure to plot on. :type fignum: figure number :param ax: axes to plot on. :type ax: axes handle + :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 plot_raw: Whether to plot the raw function p(f|y) + :type plot_raw: boolean :param linecol: color of line to plot. - :type linecol: + :type linecol: hex or color :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 - :param apply_link: apply the link function if plotting f (default false) + :type fillcol: hex or color + :param apply_link: apply the link function if plotting f (default false), as well as posterior samples if requested :type apply_link: boolean - :param samples_f: the number of posteriori f samples to plot p(f*|y) - :type samples_f: int + :param samples_y: the number of posteriori f samples to plot p(y*|y) + :type samples_y: int + :param plot_uncertain_inputs: plot the uncertainty of the inputs as error bars if they have uncertainty (BGPLVM etc.) + :type plot_uncertain_inputs: boolean + :param predict_kw: keyword args for _raw_predict and predict functions if required + :type predict_kw: dict :param plot_training_data: whether or not to plot the training points :type plot_training_data: boolean """ @@ -185,17 +191,17 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #optionally plot some samples if samples: #NOTE not tested with fixed_inputs - Ysim = model.posterior_samples(Xgrid, samples, Y_metadata=Y_metadata) - print(Ysim.shape) - print(Xnew.shape) - for yi in Ysim.T: - plots['posterior_samples'] = ax.plot(Xnew, yi[:,None], '#3300FF', linewidth=0.25) - #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. - - if samples_f: #NOTE not tested with fixed_inputs - Fsim = model.posterior_samples_f(Xgrid, samples_f) + Fsim = model.posterior_samples_f(Xgrid, samples) + if apply_link: + Fsim = model.likelihood.gp_link.transf(Fsim) for fi in Fsim.T: - plots['posterior_samples_f'] = ax.plot(Xnew, fi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) + plots['posterior_samples'] = ax.plot(Xnew, fi[:,None], '#3300FF', linewidth=0.25) + #ax.plot(Xnew, fi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. + + if samples_y: #NOTE not tested with fixed_inputs + Ysim = model.posterior_samples(Xgrid, samples_y, Y_metadata=Y_metadata) + for yi in Ysim.T: + plots['posterior_samples_y'] = ax.scatter(Xnew, yi[:,None], s=5, c=Tango.colorsHex['darkBlue'], marker='o', alpha=0.5) #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. From 5608dc8c67bc11c4080630d34a4c9ab1baabbd4a Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Thu, 10 Sep 2015 14:27:28 +0100 Subject: [PATCH 2/2] add ARD to MLP kernel --- GPy/kern/_src/mlp.py | 30 ++++++++++++++++++++++++++--- GPy/kern/_src/psi_comp/gaussherm.py | 9 ++++++++- GPy/testing/kernel_tests.py | 7 ++++++- 3 files changed, 41 insertions(+), 5 deletions(-) diff --git a/GPy/kern/_src/mlp.py b/GPy/kern/_src/mlp.py index c495b77b..40a8d2d7 100644 --- a/GPy/kern/_src/mlp.py +++ b/GPy/kern/_src/mlp.py @@ -33,9 +33,14 @@ class MLP(Kern): """ - def __init__(self, input_dim, variance=1., weight_variance=1., bias_variance=1., active_dims=None, name='mlp'): + def __init__(self, input_dim, variance=1., weight_variance=1., bias_variance=1., ARD=False, active_dims=None, name='mlp'): super(MLP, self).__init__(input_dim, active_dims, name) self.variance = Param('variance', variance, Logexp()) + self.ARD= ARD + if ARD: + wv = np.empty((input_dim,)) + wv[:] = weight_variance + weight_variance = wv self.weight_variance = Param('weight_variance', weight_variance, Logexp()) self.bias_variance = Param('bias_variance', bias_variance, Logexp()) self.link_parameters(self.variance, self.weight_variance, self.bias_variance) @@ -100,7 +105,22 @@ class MLP(Kern): X2_prod = self._comp_prod(X2) if X2 is not None else X_prod XTX = self._comp_prod(X,X2) if X2 is not None else self._comp_prod(X, X) common = var*four_over_tau/np.sqrt((X_prod[:,None]+1.)*(X2_prod[None,:]+1.)-np.square(XTX))*dL_dK - dw = (common*((XTX-b)/w-XTX*(((X_prod-b)/(w*(X_prod+1.)))[:,None]+((X2_prod-b)/(w*(X2_prod+1.)))[None,:])/2.)).sum() + if self.ARD: + if X2 is not None: + XX2 = X[:,None,:]*X2[None,:,:] if X2 is not None else X[:,None,:]*X[None,:,:] + XX = np.square(X) + X2X2 = np.square(X2) + Q = self.weight_variance.shape[0] + common_XTX = common*XTX + dw = np.dot(common.flat,XX2.reshape(-1,Q)) -( (common_XTX.sum(1)/(X_prod+1.)).T.dot(XX)+(common_XTX.sum(0)/(X2_prod+1.)).dot(X2X2))/2 + else: + XX2 = X[:,None,:]*X[None,:,:] + XX = np.square(X) + Q = self.weight_variance.shape[0] + common_XTX = common*XTX + dw = np.dot(common.flat,XX2.reshape(-1,Q)) - ((common_XTX.sum(0)+common_XTX.sum(1))/(X_prod+1.)).dot(XX)/2 + else: + dw = (common*((XTX-b)/w-XTX*(((X_prod-b)/(w*(X_prod+1.)))[:,None]+((X2_prod-b)/(w*(X2_prod+1.)))[None,:])/2.)).sum() db = (common*(1.-XTX*(1./(X_prod[:,None]+1.)+1./(X2_prod[None,:]+1.))/2.)).sum() if X2 is None: common = common+common.T @@ -118,7 +138,11 @@ class MLP(Kern): dvar = (dL_dKdiag*K).sum()/var X_prod = self._comp_prod(X) common = var*four_over_tau/(np.sqrt(1-np.square(X_prod/(X_prod+1)))*np.square(X_prod+1))*dL_dKdiag - dw = (common*(X_prod-b)).sum()/w + if self.ARD: + XX = np.square(X) + dw = np.dot(common,XX) + else: + dw = (common*(X_prod-b)).sum()/w db = common.sum() dX = common[:,None]*X*w*2 return dvar, dw, db, dX diff --git a/GPy/kern/_src/psi_comp/gaussherm.py b/GPy/kern/_src/psi_comp/gaussherm.py index 8e54e6a0..923b3eb0 100644 --- a/GPy/kern/_src/psi_comp/gaussherm.py +++ b/GPy/kern/_src/psi_comp/gaussherm.py @@ -7,11 +7,15 @@ An approximated psi-statistics implementation based on Gauss-Hermite Quadrature import numpy as np +from ....core.parameterization import Param from GPy.util.caching import Cache_this from ....util.linalg import tdot from . import PSICOMP class PSICOMP_GH(PSICOMP): + """ + TODO: support Psi2 with shape NxMxM + """ def __init__(self, degree=5, cache_K=True): self.degree = degree @@ -64,7 +68,10 @@ class PSICOMP_GH(PSICOMP): dtheta_old = kern.gradient.copy() dtheta = np.zeros_like(kern.gradient) - dZ = np.zeros_like(Z.values) + if isinstance(Z, Param): + dZ = np.zeros_like(Z.values) + else: + dZ = np.zeros_like(Z) dmu = np.zeros_like(mu) dS = np.zeros_like(S) for i in xrange(self.degree): diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 50a5aed8..cf102f99 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -252,6 +252,11 @@ class KernelGradientTestsContinuous(unittest.TestCase): continuous_kerns = ['RBF', 'Linear'] self.kernclasses = [getattr(GPy.kern, s) for s in continuous_kerns] + def test_MLP(self): + k = GPy.kern.MLP(self.D,ARD=True) + k.randomize() + self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose)) + def test_Matern32(self): k = GPy.kern.Matern32(self.D) k.randomize() @@ -464,7 +469,7 @@ class Kernel_Psi_statistics_GradientTests(unittest.TestCase): self.w3n = self.w3n+np.swapaxes(self.w3n, 1,2) def test_kernels(self): - from GPy.kern import RBF,Linear + from GPy.kern import RBF,Linear,MLP Q = self.Z.shape[1] kernels = [RBF(Q,ARD=True), Linear(Q,ARD=True)]