From 6ced5b124287c167c469d4d67776efd4eced2f11 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 12 Mar 2014 12:52:52 +0000 Subject: [PATCH] GPCoregionalizedRegresssion added --- GPy/core/gp.py | 9 ++-- GPy/likelihoods/mixed_noise.py | 58 ++++++++++++++++++++++ GPy/models/__init__.py | 6 +-- GPy/models/gp_coregionalized_regression.py | 44 ++++++++++++++++ 4 files changed, 110 insertions(+), 7 deletions(-) create mode 100644 GPy/likelihoods/mixed_noise.py create mode 100644 GPy/models/gp_coregionalized_regression.py diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 1add8268..b19f9ab2 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -27,7 +27,7 @@ class GP(Model): """ - def __init__(self, X, Y, kernel, likelihood, inference_method=None, Y_metadata=None, name='gp'): + def __init__(self, X, Y, kernel, likelihood, inference_method=None, name='gp', **Y_metadata): super(GP, self).__init__(name) assert X.ndim == 2 @@ -43,7 +43,7 @@ class GP(Model): _, self.output_dim = self.Y.shape if Y_metadata is not None: - self.Y_metadata = ObservableArray(Y_metadata) + self.Y_metadata = Y_metadata else: self.Y_metadata = None @@ -56,7 +56,7 @@ class GP(Model): #find a sensible inference method if inference_method is None: - if isinstance(likelihood, likelihoods.Gaussian): + if isinstance(likelihood, likelihoods.Gaussian) or isinstance(likelihood, likelihoods.MixedNoise): inference_method = exact_gaussian_inference.ExactGaussianInference() else: inference_method = expectation_propagation @@ -67,7 +67,8 @@ class GP(Model): self.add_parameter(self.likelihood) def parameters_changed(self): - self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) + self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, **self.Y_metadata) + self.likelihood.update_gradients(np.diag(grad_dict['dL_dK']), **self.Y_metadata) self.kern.update_gradients_full(grad_dict['dL_dK'], self.X) def log_likelihood(self): diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py new file mode 100644 index 00000000..b60f3adf --- /dev/null +++ b/GPy/likelihoods/mixed_noise.py @@ -0,0 +1,58 @@ +import numpy as np +from scipy import stats, special +from GPy.util.univariate_Gaussian import std_norm_pdf, std_norm_cdf +import link_functions +from likelihood import Likelihood +from ..core.parameterization import Param +from ..core.parameterization.transformations import Logexp +from ..core.parameterization import Parameterized +import itertools + +class MixedNoise(Likelihood): + def __init__(self, likelihoods_list, noise_index, variance = None, name='mixed_noise'): + + Nlike = len(likelihoods_list) + self.order = np.unique(noise_index) + + assert self.order.size == Nlike + + if variance is None: + variance = np.ones(Nlike) + else: + assert variance.size == Nlike + + super(Likelihood, self).__init__(name=name) + + self.add_parameters(*likelihoods_list) + self.likelihoods_list = likelihoods_list + self.noise_index = noise_index + self.log_concave = False + self.likelihoods_indices = [noise_index.flatten()==j for j in self.order] + + def covariance_matrix(self, Y, noise_index, **Y_metadata): + variance = np.zeros(Y.shape[0]) + for lik, ind in itertools.izip(self.likelihoods_list, self.likelihoods_indices): + variance[ind] = lik.variance + return np.diag(variance) + + def update_gradients(self, partial, noise_index, **Y_metadata): + [lik.update_gradients(partial[ind]) for lik,ind in itertools.izip(self.likelihoods_list, self.likelihoods_indices)] + + def predictive_values(self, mu, var, full_cov=False, noise_index=None, **Y_metadata): + _variance = np.array([ self.likelihoods_list[j].variance for j in noise_index ]) + 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 + + def predictive_variance(self, mu, sigma, noise_index, predictive_mean=None, **Y_metadata): + if isinstance(noise_index,int): + _variance = self.variance[noise_index] + else: + _variance = np.array([ self.variance[j] for j in noise_index ])[:,None] + return _variance + sigma**2 diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py index 83db4b8c..a253c63d 100644 --- a/GPy/models/__init__.py +++ b/GPy/models/__init__.py @@ -13,6 +13,6 @@ from warped_gp import WarpedGP from bayesian_gplvm import BayesianGPLVM from mrd import MRD from gradient_checker import GradientChecker -from gp_multioutput_regression import GPMultioutputRegression -from sparse_gp_multioutput_regression import SparseGPMultioutputRegression -from ss_gplvm import SSGPLVM \ No newline at end of file +from ss_gplvm import SSGPLVM +from gp_coregionalized_regression import GPCoregionalizedRegression +from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression diff --git a/GPy/models/gp_coregionalized_regression.py b/GPy/models/gp_coregionalized_regression.py new file mode 100644 index 00000000..313e09d4 --- /dev/null +++ b/GPy/models/gp_coregionalized_regression.py @@ -0,0 +1,44 @@ +# 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 GP +from .. import likelihoods +from .. import kern +from .. import util + +class GPCoregionalizedRegression(GP): + """ + Gaussian Process model for heteroscedastic multioutput regression + + This is a thin wrapper around the models.GP 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 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 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, kernel=None, likelihoods_list=None, name='GPCR',W_rank=1,kernel_name='X'): + + #Input and Output + X,Y,self.noise_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.noise_index,likelihoods_list) + + super(GPCoregionalizedRegression, self).__init__(X,Y,kernel,likelihood, noise_index=self.noise_index)