From 59ae2f0e34879225b25192d8f9231450f417e967 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:00:23 +0000 Subject: [PATCH 1/7] coregionalization examples fixed --- GPy/examples/regression.py | 63 ++++++++++---------------------------- 1 file changed, 17 insertions(+), 46 deletions(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 190af93b..2a4b91b3 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -25,80 +25,51 @@ def olympic_marathon_men(optimize=True, plot=True): return m -def coregionalization_toy2(optimize=True, plot=True): +def coregionalization_toy(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions. """ #build a design matrix with a column of integers indicating the output X1 = np.random.rand(50, 1) * 8 X2 = np.random.rand(30, 1) * 5 - index = np.vstack((np.zeros_like(X1), np.ones_like(X2))) - X = np.hstack((np.vstack((X1, X2)), index)) #build a suitable set of observed variables Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. - Y = np.vstack((Y1, Y2)) - #build the kernel - k1 = GPy.kern.RBF(1) + GPy.kern.Bias(1) - k2 = GPy.kern.Coregionalize(2,1) - k = k1**k2 - m = GPy.models.GPRegression(X, Y, kernel=k) + m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) if optimize: m.optimize('bfgs', max_iters=100) if plot: - m.plot(fixed_inputs=[(1,0)]) - m.plot(fixed_inputs=[(1,1)], ax=pb.gca()) - + slices = GPy.util.multioutput.get_slices([X1,X2]) + m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) + m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) return m -#FIXME: Needs recovering once likelihoods are consolidated -#def coregionalization_toy(optimize=True, plot=True): -# """ -# A simple demonstration of coregionalization on two sinusoidal functions. -# """ -# X1 = np.random.rand(50, 1) * 8 -# X2 = np.random.rand(30, 1) * 5 -# X = np.vstack((X1, X2)) -# Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 -# Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05 -# Y = np.vstack((Y1, Y2)) -# -# k1 = GPy.kern.RBF(1) -# m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1]) -# m.constrain_fixed('.*rbf_var', 1.) -# m.optimize(max_iters=100) -# -# fig, axes = pb.subplots(2,1) -# m.plot(fixed_inputs=[(1,0)],ax=axes[0]) -# m.plot(fixed_inputs=[(1,1)],ax=axes[1]) -# axes[0].set_title('Output 0') -# axes[1].set_title('Output 1') -# return m - def coregionalization_sparse(optimize=True, plot=True): """ A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations. """ - #fetch the data from the non sparse examples - m = coregionalization_toy2(optimize=False, plot=False) - X, Y = m.X, m.Y + #build a design matrix with a column of integers indicating the output + X1 = np.random.rand(50, 1) * 8 + X2 = np.random.rand(30, 1) * 5 - k = GPy.kern.RBF(1)**GPy.kern.Coregionalize(2) + #build a suitable set of observed variables + Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05 + Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2. - #construct a model - m = GPy.models.SparseGPRegression(X,Y, num_inducing=25, kernel=k) - m.Z[:,1].fix() # don't optimize the inducing input indexes + m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2]) if optimize: - m.optimize('bfgs', max_iters=100, messages=1) + m.optimize('bfgs', max_iters=100) if plot: - m.plot(fixed_inputs=[(1,0)]) - m.plot(fixed_inputs=[(1,1)], ax=pb.gca()) + slices = GPy.util.multioutput.get_slices([X1,X2]) + m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0}) + m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca()) + pb.ylim(-3,) return m From ced160c0366ee4eb475e5c8290055dd1f7d11f7b Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:01:21 +0000 Subject: [PATCH 2/7] default None for Y_metadata in predictive_quantiles --- GPy/likelihoods/gaussian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 0c73e485..51b56d09 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -99,7 +99,7 @@ class Gaussian(Likelihood): def predictive_variance(self, mu, sigma, predictive_mean=None): return self.variance + sigma**2 - def predictive_quantiles(self, mu, var, quantiles, Y_metadata): + def predictive_quantiles(self, mu, var, quantiles, Y_metadata=None): return [stats.norm.ppf(q/100.)*np.sqrt(var) + mu for q in quantiles] def pdf_link(self, link_f, y, extra_data=None): From ef31b5f1c9d23565f566ae97c3abfb0dc00e7ef5 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:02:00 +0000 Subject: [PATCH 3/7] Lines not used deleted --- GPy/likelihoods/likelihood.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index e7bad74a..6872320f 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -407,10 +407,7 @@ class Likelihood(Parameterized): #compute the quantiles by sampling!!! N_samp = 1000 s = np.random.randn(mu.shape[0], N_samp)*np.sqrt(var) + mu - #ss_f = s.flatten() - #ss_y = self.samples(ss_f, Y_metadata) ss_y = self.samples(s, Y_metadata) - #ss_y = ss_y.reshape(mu.shape[0], N_samp) return [np.percentile(ss_y ,q, axis=1)[:,None] for q in quantiles] From b44fee93c4cf78ea5b8b3a3a9d8aeb00590bed85 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:02:36 +0000 Subject: [PATCH 4/7] function predictive_quantiles added --- GPy/likelihoods/mixed_noise.py | 42 +++++++++++++++------------------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py index b4960f3a..909d4002 100644 --- a/GPy/likelihoods/mixed_noise.py +++ b/GPy/likelihoods/mixed_noise.py @@ -11,7 +11,7 @@ import itertools class MixedNoise(Likelihood): def __init__(self, likelihoods_list, name='mixed_noise'): - + #NOTE at the moment this likelihood only works for using a list of gaussians super(Likelihood, self).__init__(name=name) self.add_parameters(*likelihoods_list) @@ -38,35 +38,32 @@ class MixedNoise(Likelihood): return np.array([dL_dKdiag[ind==i].sum() for i in range(len(self.likelihoods_list))]) def predictive_values(self, mu, var, full_cov=False, Y_metadata=None): - if all([isinstance(l, Gaussian) for l in self.likelihoods_list]): - ind = Y_metadata['output_index'].flatten() - _variance = np.array([self.likelihoods_list[j].variance for j in ind ]) - if full_cov: - var += np.eye(var.shape[0])*_variance - else: - var += _variance - return mu, var + ind = Y_metadata['output_index'].flatten() + _variance = np.array([self.likelihoods_list[j].variance for j in ind ]) + if full_cov: + var += np.eye(var.shape[0])*_variance else: - raise NotImplementedError + var += _variance + return mu, var - def predictive_variance(self, mu, sigma, **other_shit): - if isinstance(noise_index,int): - _variance = self.variance[noise_index] - else: - _variance = np.array([ self.variance[j] for j in noise_index ])[:,None] + def predictive_variance(self, mu, sigma, Y_metadata): + _variance = self.gaussian_variance(Y_metadata) return _variance + sigma**2 + def predictive_quantiles(self, mu, var, quantiles, Y_metadata): + ind = Y_metadata['output_index'].flatten() + outputs = np.unique(ind) + Q = np.zeros( (mu.size,len(quantiles)) ) + for j in outputs: + q = self.likelihoods_list[j].predictive_quantiles(mu[ind==j,:], + var[ind==j,:],quantiles,Y_metadata=None) + Q[ind==j,:] = np.hstack(q) + return [q[:,None] for q in Q.T] def covariance_matrix(self, Y, Y_metadata): - #assert all([isinstance(l, Gaussian) for l in self.likelihoods_list]) - #ind = Y_metadata['output_index'].flatten() - #variance = np.zeros(Y.shape[0]) - #for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))): - # variance[ind==j] = lik.variance - #return np.diag(variance) + #TODO make more general, to allow non-gaussian likelihoods return np.diag(self.gaussian_variance(Y_metadata).flatten()) - def samples(self, gp, Y_metadata): """ Returns a set of samples of observations based on a given value of the latent variable. @@ -84,4 +81,3 @@ class MixedNoise(Likelihood): _ysim = np.array([np.random.normal(lik.gp_link.transf(gpj), scale=np.sqrt(lik.variance), size=1) for gpj in gp_filtered.flatten()]) Ysim[flt,:] = _ysim.reshape(n1,N2) return Ysim - From 1a0e16a6f38a43d154b0b98114451a71f624199f Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:03:16 +0000 Subject: [PATCH 5/7] bug fixed --- GPy/models/gp_coregionalized_regression.py | 2 +- GPy/models/sparse_gp_coregionalized_regression.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/models/gp_coregionalized_regression.py b/GPy/models/gp_coregionalized_regression.py index 6d478fd9..5854fe63 100644 --- a/GPy/models/gp_coregionalized_regression.py +++ b/GPy/models/gp_coregionalized_regression.py @@ -36,7 +36,7 @@ class GPCoregionalizedRegression(GP): #Kernel if kernel is None: - kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name) + kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name) #Likelihood likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) diff --git a/GPy/models/sparse_gp_coregionalized_regression.py b/GPy/models/sparse_gp_coregionalized_regression.py index a97696d2..e6c2c4c9 100644 --- a/GPy/models/sparse_gp_coregionalized_regression.py +++ b/GPy/models/sparse_gp_coregionalized_regression.py @@ -43,7 +43,7 @@ class SparseGPCoregionalizedRegression(SparseGP): #Kernel if kernel is None: - kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name) + kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name) #Likelihood likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) From 31ad9b73bfc8376a24ec277ebd61ce4454af4cfd Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 18 Mar 2014 16:04:01 +0000 Subject: [PATCH 6/7] 1D inducing inputs modified for coregionalized models --- GPy/plotting/matplot_dep/models_plots.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index cbb213b1..b626758f 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -123,6 +123,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #add inducing inputs (if a sparse model is used) if hasattr(model,"Z"): #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] + if isinstance(model,SparseGPCoregionalizedRegression): + Z = Z[Z[:,-1] == Y_metadata['output_index'],:] Zu = Z[:,free_dims] z_height = ax.get_ylim()[0] plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) From dfb555bbf256a07a737fd2cddaaa5171b4dd9dc5 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Thu, 20 Mar 2014 16:50:17 +0000 Subject: [PATCH 7/7] missing docstrings --- .../latent_function_inference/expectation_propagation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py index 514a6dc7..ff60d2e3 100644 --- a/GPy/inference/latent_function_inference/expectation_propagation.py +++ b/GPy/inference/latent_function_inference/expectation_propagation.py @@ -11,9 +11,9 @@ class EP(object): :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) :type epsilon: float - :param eta: Power EP thing TODO: Ricardo: what, exactly? + :param eta: parameter for fractional EP updates. :type eta: float64 - :param delta: Power EP thing TODO: Ricardo: what, exactly? + :param delta: damping EP updates factor. :type delta: float64 """ self.epsilon, self.eta, self.delta = epsilon, eta, delta