From 4692880cef3d891d772803389d97e08d41009e8d Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:33:58 +0000 Subject: [PATCH 1/9] parameter missin in dL_dthetaL added --- .../latent_function_inference/exact_gaussian_inference.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py index e76575c6..c4ed0a0f 100644 --- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py +++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py @@ -49,6 +49,6 @@ class ExactGaussianInference(object): dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi) - dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK)) + dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK),Y_metadata) return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL} From d4735e54c83e0d0f7db5671ca3042cbc2f6544cc Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 10:34:08 +0000 Subject: [PATCH 2/9] minor changes --- GPy/examples/coreg_example.py | 37 +++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py index 967758c6..f0288f35 100644 --- a/GPy/examples/coreg_example.py +++ b/GPy/examples/coreg_example.py @@ -3,6 +3,42 @@ import pylab as pb import GPy pb.ion() +X1 = 100 * np.random.rand(3)[:,None] +X2 = 100 * np.random.rand(4)[:,None] +I1 = np.zeros_like(X1) +I2 = np.ones_like(X2) + +_X = np.vstack([ X1, X2 ]) +_I = np.vstack([ I1, I2 ]) + +X = np.hstack([ _X, _I ]) + +Bias = GPy.kern.Bias(1,active_dims=[0]) +Coreg = GPy.kern.Coregionalize(1,2,active_dims=[1]) +K = Bias.prod(Coreg,name='X') + +K.coregion.W = 0 +print K.coregion.W + +print Bias.K(_X,_X) +print K.K(X,X) + +pb.matshow(K.K(X,X)) + +stop + +Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")] +kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=12,kernels_list=Mlist,name='H') + + +m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern) +m.optimize() + + + + +""" + X1 = 100 * np.random.rand(100)[:,None] X2 = 100 * np.random.rand(100)[:,None] #X1.sort() @@ -28,3 +64,4 @@ slices = GPy.util.multioutput.get_slices([Y1,Y2]) m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0) m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1) +""" From fe82678d80a44d6195ca11a2b04ad464f58a62a0 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 18:57:34 +0000 Subject: [PATCH 3/9] Changes to allow heteroscedastic inference --- .../latent_function_inference/var_dtc.py | 37 ++++++++++--------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 7a0c14e8..f948ff38 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -2,7 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) from posterior import Posterior -from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify +from ...util.linalg import mdot, jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify from ...util import diag from ...core.parameterization.variational import VariationalPosterior import numpy as np @@ -48,7 +48,7 @@ class VarDTC(object): def get_VVTfactor(self, Y, prec): return Y * prec # TODO chache this, and make it effective - def inference(self, kern, X, Z, likelihood, Y): + def inference(self, kern, X, Z, likelihood, Y, Y_metadata=None): if isinstance(X, VariationalPosterior): uncertain_inputs = True psi0 = kern.psi0(Z, X) @@ -65,7 +65,9 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./np.fmax(likelihood.variance, 1e-6) + #beta = 1./np.fmax(likelihood.variance, 1e-6) + beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6) + # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) @@ -74,7 +76,7 @@ class VarDTC(object): trYYT = self.get_trYYT(Y) # do the inference: - het_noise = beta.size < 1 + het_noise = beta.size > 1 num_inducing = Z.shape[0] num_data = Y.shape[0] # kernel computations, using BGPLVM notation @@ -134,16 +136,16 @@ class VarDTC(object): # log marginal likelihood log_marginal = _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, - psi0, A, LB, trYYT, data_fit) + psi0, A, LB, trYYT, data_fit, Y) #put the gradients in the right places dL_dR = _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, - data_fit, num_data, output_dim, trYYT) + data_fit, num_data, output_dim, trYYT, Y) - dL_dthetaL = likelihood.exact_inference_gradients(dL_dR) + dL_dthetaL = likelihood.exact_inference_gradients(dL_dR,Y_metadata) if uncertain_inputs: grad_dict = {'dL_dKmm': dL_dKmm, @@ -385,7 +387,7 @@ def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, C return dL_dpsi0, dL_dpsi1, dL_dpsi2 -def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, data_fit, num_data, output_dim, trYYT): +def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, DBi_plus_BiPBi, Lm, A, psi0, psi1, beta, data_fit, num_data, output_dim, trYYT, Y): # the partial derivative vector for the likelihood if likelihood.size == 0: # save computation here. @@ -394,19 +396,20 @@ def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, if uncertain_inputs: raise NotImplementedError, "heteroscedatic derivates with uncertain inputs not implemented" else: - from ...util.linalg import chol_inv - LBi = chol_inv(LB) + #from ...util.linalg import chol_inv + #LBi = chol_inv(LB) + LBi, _ = dtrtrs(LB,np.eye(LB.shape[0])) + Lmi_psi1, nil = dtrtrs(Lm, psi1.T, lower=1, trans=0) _LBi_Lmi_psi1, _ = dtrtrs(LB, Lmi_psi1, lower=1, trans=0) - dL_dR = -0.5 * beta + 0.5 * likelihood.V**2 + dL_dR = -0.5 * beta + 0.5 * (beta*Y)**2 dL_dR += 0.5 * output_dim * (psi0 - np.sum(Lmi_psi1**2,0))[:,None] * beta**2 dL_dR += 0.5*np.sum(mdot(LBi.T,LBi,Lmi_psi1)*Lmi_psi1,0)[:,None]*beta**2 - dL_dR += -np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T * likelihood.Y * beta**2 + dL_dR += -np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T * Y * beta**2 dL_dR += 0.5*np.dot(_LBi_Lmi_psi1Vf.T,_LBi_Lmi_psi1).T**2 * beta**2 - else: # likelihood is not heteroscedatic dL_dR = -0.5 * num_data * output_dim * beta + 0.5 * trYYT * beta ** 2 @@ -414,11 +417,11 @@ def _compute_dL_dR(likelihood, het_noise, uncertain_inputs, LB, _LBi_Lmi_psi1Vf, dL_dR += beta * (0.5 * np.sum(A * DBi_plus_BiPBi) - data_fit) return dL_dR -def _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit): -#compute log marginal likelihood +def _compute_log_marginal_likelihood(likelihood, num_data, output_dim, beta, het_noise, psi0, A, LB, trYYT, data_fit,Y): + #compute log marginal likelihood if het_noise: - lik_1 = -0.5 * num_data * output_dim * np.log(2. * np.pi) + 0.5 * np.sum(np.log(beta)) - 0.5 * np.sum(likelihood.V * likelihood.Y) - lik_2 = -0.5 * output_dim * (np.sum(beta * psi0) - np.trace(A)) + lik_1 = -0.5 * num_data * output_dim * np.log(2. * np.pi) + 0.5 * np.sum(np.log(beta)) - 0.5 * np.sum(beta * Y**2) + lik_2 = -0.5 * output_dim * (np.sum(beta.flatten() * psi0) - np.trace(A)) else: lik_1 = -0.5 * num_data * output_dim * (np.log(2. * np.pi) - np.log(beta)) - 0.5 * beta * trYYT lik_2 = -0.5 * output_dim * (np.sum(beta * psi0) - np.trace(A)) From 7fdecf5e317a40461d1f5e31f3a9d912d74972c4 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 18:57:58 +0000 Subject: [PATCH 4/9] bug fixed --- GPy/core/gp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 70b7d695..d7f07a47 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -43,7 +43,7 @@ class GP(Model): _, self.output_dim = self.Y.shape if Y_metadata is None: - Y_metadata = {} + self.Y_metadata = {} else: self.Y_metadata = Y_metadata From 15901a48d43566eef9995eeeed7ee9df6eb3c146 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 18:58:31 +0000 Subject: [PATCH 5/9] Changes to allow compatibility with mixed noise likelihoods --- GPy/core/sparse_gp.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index d137ceff..a0b09564 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -31,7 +31,7 @@ class SparseGP(GP): """ - def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp'): + def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp', Y_metadata=None): #pick a sensible inference method if inference_method is None: @@ -45,7 +45,7 @@ class SparseGP(GP): self.Z = Param('inducing inputs', Z) self.num_inducing = Z.shape[0] - GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name) + GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata) self.add_parameter(self.Z, index=0) @@ -53,7 +53,7 @@ class SparseGP(GP): return isinstance(self.X, VariationalPosterior) def parameters_changed(self): - self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y) + self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y, self.Y_metadata) self.likelihood.update_gradients(self.grad_dict['dL_dthetaL']) if isinstance(self.X, VariationalPosterior): #gradients wrt kernel @@ -75,7 +75,6 @@ class SparseGP(GP): target += self.kern.gradient self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None) self.kern.gradient += target - #gradients wrt Z self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) From bea74ac5231e47cce34170273771fd61a39e6b73 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 18:59:37 +0000 Subject: [PATCH 6/9] new function added --- GPy/likelihoods/gaussian.py | 6 ++++++ GPy/likelihoods/mixed_noise.py | 30 +++++++++++++++++++----------- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 101aac4b..0c73e485 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -50,6 +50,12 @@ class Gaussian(Likelihood): if isinstance(gp_link, link_functions.Identity): self.log_concave = True + def betaY(self,Y,Y_metadata=None): + return Y/self.gaussian_variance(Y_metadata) + + def gaussian_variance(self, Y_metadata=None): + return self.variance + def covariance_matrix(self, Y, Y_metadata=None): return np.eye(Y.shape[0]) * self.variance diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py index bfcb5916..b4960f3a 100644 --- a/GPy/likelihoods/mixed_noise.py +++ b/GPy/likelihoods/mixed_noise.py @@ -18,6 +18,17 @@ class MixedNoise(Likelihood): self.likelihoods_list = likelihoods_list self.log_concave = False + def gaussian_variance(self, Y_metadata): + assert all([isinstance(l, Gaussian) for l in self.likelihoods_list]) + ind = Y_metadata['output_index'].flatten() + variance = np.zeros(ind.size) + for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))): + variance[ind==j] = lik.variance + return variance[:,None] + + def betaY(self,Y,Y_metadata): + return Y/self.gaussian_variance(Y_metadata=Y_metadata) + def update_gradients(self, gradients): self.gradient = gradients @@ -32,13 +43,9 @@ class MixedNoise(Likelihood): _variance = np.array([self.likelihoods_list[j].variance for j in ind ]) if full_cov: var += np.eye(var.shape[0])*_variance - #d = 2*np.sqrt(np.diag(var)) - #low, up = mu - d, mu + d else: var += _variance - #d = 2*np.sqrt(var) - #low, up = mu - d, mu + d - return mu, var#, low, up + return mu, var else: raise NotImplementedError @@ -51,12 +58,13 @@ class MixedNoise(Likelihood): 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) + #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) + return np.diag(self.gaussian_variance(Y_metadata).flatten()) def samples(self, gp, Y_metadata): From 385ce7d344ed00a79ca79b2b8ba792bb4256cc68 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 19:00:16 +0000 Subject: [PATCH 7/9] Changes to allow mixed noise likelihoods --- GPy/plotting/matplot_dep/models_plots.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index ae79569b..cbb213b1 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -7,6 +7,7 @@ import Tango from base_plots import gpplot, x_frame1D, x_frame2D from ...util.misc import param_to_array from ...models.gp_coregionalized_regression import GPCoregionalizedRegression +from ...models.sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression def plot_fit(model, plot_limits=None, which_data_rows='all', @@ -86,7 +87,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', lower = m - 2*np.sqrt(v) upper = m + 2*np.sqrt(v) else: - meta = {'output_index': Xgrid[:,-1:].astype(np.int)} if isinstance(model,GPCoregionalizedRegression) else None + if isinstance(model,GPCoregionalizedRegression) or isinstance(model,SparseGPCoregionalizedRegression): + meta = {'output_index': Xgrid[:,-1:].astype(np.int)} + else: + meta = None m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta) lower, upper = model.predict_quantiles(Xgrid, Y_metadata=meta) From d506604c58e6277044e083e2a0f14392ebfc02db Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 19:00:35 +0000 Subject: [PATCH 8/9] minor changes --- GPy/examples/coreg_example.py | 39 +++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 13 deletions(-) diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py index f0288f35..8f4cfcc6 100644 --- a/GPy/examples/coreg_example.py +++ b/GPy/examples/coreg_example.py @@ -2,9 +2,10 @@ import numpy as np import pylab as pb import GPy pb.ion() +pb.close('all') -X1 = 100 * np.random.rand(3)[:,None] -X2 = 100 * np.random.rand(4)[:,None] +X1 = np.arange(3)[:,None] +X2 = np.arange(4)[:,None] I1 = np.zeros_like(X1) I2 = np.ones_like(X2) @@ -13,27 +14,39 @@ _I = np.vstack([ I1, I2 ]) X = np.hstack([ _X, _I ]) +Y1 = np.sin(X1/8.) +Y2 = np.cos(X2/8.) + Bias = GPy.kern.Bias(1,active_dims=[0]) Coreg = GPy.kern.Coregionalize(1,2,active_dims=[1]) K = Bias.prod(Coreg,name='X') -K.coregion.W = 0 -print K.coregion.W +#K.coregion.W = 0 +#print K.coregion.W -print Bias.K(_X,_X) -print K.K(X,X) +#print Bias.K(_X,_X) +#print K.K(X,X) -pb.matshow(K.K(X,X)) +#pb.matshow(K.K(X,X)) -stop Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")] -kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=12,kernels_list=Mlist,name='H') - - -m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern) -m.optimize() +kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H') +kern.B.W = 0 +kern.B.kappa = 1. +#kern.B.W.fix() +#kern.B.kappa.fix() +#m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern) +m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], kernel=kern) +#m.optimize() +m.checkgrad(verbose=1) +fig = pb.figure() +ax0 = fig.add_subplot(211) +ax1 = fig.add_subplot(212) +slices = GPy.util.multioutput.get_slices([Y1,Y2]) +m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0) +#m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1) From 9cda9f991489332b51dc731bd2019deb24e6659c Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 17 Mar 2014 19:01:16 +0000 Subject: [PATCH 9/9] New model SparseGPCoregionalizedRegression --- GPy/models/__init__.py | 1 + .../sparse_gp_coregionalized_regression.py | 66 +++++++++++++++++++ 2 files changed, 67 insertions(+) create mode 100644 GPy/models/sparse_gp_coregionalized_regression.py diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py index 34e5a17e..d0988c9e 100644 --- a/GPy/models/__init__.py +++ b/GPy/models/__init__.py @@ -15,4 +15,5 @@ from mrd import MRD from gradient_checker import GradientChecker from ss_gplvm import SSGPLVM from gp_coregionalized_regression import GPCoregionalizedRegression +from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression #.py file not included!!! #from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression diff --git a/GPy/models/sparse_gp_coregionalized_regression.py b/GPy/models/sparse_gp_coregionalized_regression.py new file mode 100644 index 00000000..a97696d2 --- /dev/null +++ b/GPy/models/sparse_gp_coregionalized_regression.py @@ -0,0 +1,66 @@ +# Copyright (c) 2012 - 2014 the GPy Austhors (see AUTHORS.txt) +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +from ..core import SparseGP +from ..inference.latent_function_inference import VarDTC +from .. import likelihoods +from .. import kern +from .. import util + +class SparseGPCoregionalizedRegression(SparseGP): + """ + Sparse Gaussian Process model for heteroscedastic multioutput regression + + This is a thin wrapper around the SparseGP class, with a set of sensible defaults + + :param X_list: list of input observations corresponding to each output + :type X_list: list of numpy arrays + :param Y_list: list of observed values related to the different noise models + :type Y_list: list of numpy arrays + :param Z_list: list of inducing inputs (optional) + :type Z_list: empty list | list of numpy arrays + :param kernel: a GPy kernel, defaults to RBF ** Coregionalized + :type kernel: None | GPy.kernel defaults + :likelihoods_list: a list of likelihoods, defaults to list of Gaussian likelihoods + :type likelihoods_list: None | a list GPy.likelihoods + :param num_inducing: number of inducing inputs, defaults to 10 per output (ignored if Z_list is not empty) + :type num_inducing: integer | list of integers + + :param name: model name + :type name: string + :param W_rank: number tuples of the corregionalization parameters 'W' (see coregionalize kernel documentation) + :type W_rank: integer + :param kernel_name: name of the kernel + :type kernel_name: string + """ + + def __init__(self, X_list, Y_list, Z_list=[], kernel=None, likelihoods_list=None, num_inducing=10, X_variance=None, name='SGPCR',W_rank=1,kernel_name='X'): + + #Input and Output + X,Y,self.output_index = util.multioutput.build_XY(X_list,Y_list) + Ny = len(Y_list) + + #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) + + #Likelihood + likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list) + + #Inducing inputs list + if len(Z_list): + assert len(Z_list) == self.output_dim, 'Number of outputs do not match length of inducing inputs list.' + else: + if isinstance(num_inducing,np.int): + num_inducing = [num_inducing] * Ny + num_inducing = np.asarray(num_inducing) + assert num_inducing.size == Ny, 'Number of outputs do not match length of inducing inputs list.' + for ni,Xi in zip(num_inducing,X_list): + i = np.random.permutation(Xi.shape[0])[:ni] + Z_list.append(Xi[i].copy()) + + Z, _, Iz = util.multioutput.build_XY(Z_list) + + super(SparseGPCoregionalizedRegression, self).__init__(X, Y, Z, kernel, likelihood, inference_method=VarDTC(), Y_metadata={'output_index':self.output_index}) + self['.*inducing'][:,-1].fix()