diff --git a/GPy/core/gp.py b/GPy/core/gp.py index fb5fe789..4edde5cd 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -185,7 +185,7 @@ class GP(GPBase): if isinstance(self.likelihood,EP_Mixed_Noise): mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output) else: - mean, var, _025pm, _975pm = self.likelihood_list[output].predictive_values(mu, var, full_cov) + mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output) return mean, var, _025pm, _975pm def _raw_predict_single_output(self, _Xnew, output=0, which_parts='all', full_cov=False,stop=False): diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index 51a9352c..fe297d6b 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -178,7 +178,7 @@ class GPBase(Model): 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[self.likelihood.index==output][:,None], 'kx', mew=1.5) - ax.plot(Xu[which_data], self.likelihood_list[output].data, 'kx', mew=1.5) + 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) diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 697f3554..e05d0688 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -250,6 +250,19 @@ def symmetric(k): return k_ def coregionalise(Nout,R=1, W=None, kappa=None): + """ + Construct coregionalisation kernel, based on the coregionlisation matrix B = np.dot(W,W.T) + np.eye(Nout)*kappa + + :param Nout: the number of outputs to corregionalise + :type Nout: int + :param R: the number of columns in the W matrix + :type R: int + :param W: W matrix + :type W: numpy array of dimensionality (Nout x R) + :param kappa: kappa vector + :type kappa: numpy array of dimensionality (Nout,) + + """ p = parts.coregionalise.Coregionalise(Nout,R,W,kappa) return kern(1,[p]) diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py index 55f437b1..0cb62eb0 100644 --- a/GPy/likelihoods/__init__.py +++ b/GPy/likelihoods/__init__.py @@ -1,6 +1,7 @@ from ep import EP from ep_mixed_noise import EP_Mixed_Noise from gaussian import Gaussian +from gaussian_mixed_noise import Gaussian_Mixed_Noise from noise_model_constructors import * # TODO: from Laplace import Laplace diff --git a/GPy/likelihoods/gaussian_mixed_noise.py b/GPy/likelihoods/gaussian_mixed_noise.py new file mode 100644 index 00000000..cf5e2ca8 --- /dev/null +++ b/GPy/likelihoods/gaussian_mixed_noise.py @@ -0,0 +1,99 @@ +# Copyright (c) 2013, Ricardo Andrade +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + +import numpy as np +from scipy import stats +from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs +from likelihood import likelihood +from . import Gaussian + + +class Gaussian_Mixed_Noise(likelihood): + def __init__(self, data_list, noise_params=None, normalize=True): + + if noise_params is None: + noise_params = [1.] * len(data_list) + + assert len(data_list) == len(noise_params) + + self.noise_model_list = [Gaussian(Y,variance=v,normalize = normalize) for Y,v in zip(data_list,noise_params)] + + + self.n_list = [data.size for data in data_list] + n_models = len(data_list) + self.n_params = [noise_model._get_params().size for noise_model in self.noise_model_list] + + self.index = np.vstack([np.repeat(i,n)[:,None] for i,n in zip(range(n_models),self.n_list)]) + + + self.data = np.vstack(data_list) + self.N, self.output_dim = self.data.shape + self._offset = np.zeros((1, self.output_dim)) + self._scale = np.ones((1, self.output_dim)) + + self.is_heteroscedastic = True #TODO check how to deal with this + self.Z = 0. # a correction factor which accounts for the approximation made + + self.set_data(data_list) + #self._variance = np.asarray(variance) + 1. + self._set_params(np.asarray(noise_params)) + + def set_data(self, data_list): + self.data = np.vstack(data_list) + self.N, D = self.data.shape + assert D == self.output_dim + self.Y = (self.data - self._offset) / self._scale + if D > self.N: + self.YYT = np.dot(self.Y, self.Y.T) + self.trYYT = np.trace(self.YYT) + self.YYT_factor = jitchol(self.YYT) + else: + self.YYT = None + self.trYYT = np.sum(np.square(self.Y)) + self.YYT_factor = self.Y + + def predictive_values(self,mu,var,full_cov,noise_model): + """ + Predicts the output given the GP + + :param mu: GP's mean + :param var: GP's variance + :param full_cov: whether to return the full covariance matrix, or just the diagonal + :type full_cov: False|True + :param noise_model: noise model to use + :type noise_model: integer + """ + if full_cov: + raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood" + return self.noise_model_list[noise_model].predictive_values(mu,var,full_cov) + + def _get_params(self): + return np.hstack([noise_model._get_params().flatten() for noise_model in self.noise_model_list]) + + def _get_param_names(self): + if len(self.noise_model_list) == 1: + names = self.noise_model_list[0]._get_param_names() + else: + names = [] + for noise_model,i in zip(self.noise_model_list,range(len(self.n_list))): + names.append(''.join(noise_model._get_param_names() + ['_%s' %i])) + return names + + def _set_params(self,p): + cs_params = np.cumsum([0]+self.n_params) + for i in range(len(self.n_params)): + self.noise_model_list[i]._set_params(p[cs_params[i]:cs_params[i+1]]) + self.precision = np.hstack([np.repeat(noise_model.precision,n) for noise_model,n in zip(self.noise_model_list,self.n_list)])[:,None] + self.V = (self.precision) * self.Y + self.VVT_factor = self.precision * self.YYT_factor + self.covariance_matrix = np.eye(self.N) * 1./self.precision + #self._variance = x + + def _gradients(self,partial): + #NOTE this is not tested + return np.hstack([noise_model._gradients(partial) for noise_model in self.noise_model_list]) + + + + diff --git a/GPy/models/gp_multioutput.py b/GPy/models/gp_multioutput.py index 96138bbf..bc696cbd 100644 --- a/GPy/models/gp_multioutput.py +++ b/GPy/models/gp_multioutput.py @@ -31,6 +31,7 @@ class GPMultioutput(GP): """ def __init__(self,X_list,Y_list,kernel_list=None,normalize_X=False,normalize_Y=False,W=1,mixed_noise_list=[]): #TODO W + #TODO: split into 2 models gp_mixed_noise and ep_mixed_noise assert len(X_list) == len(Y_list) index = [] @@ -41,7 +42,8 @@ class GPMultioutput(GP): i += 1 index = np.vstack(index) - self.likelihood_list = [] + """ + if mixed_noise_list == []: for Y in Y_list: self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y)) @@ -49,14 +51,18 @@ class GPMultioutput(GP): Y = np.vstack([l_.Y for l_ in self.likelihood_list]) likelihood = likelihoods.Gaussian(Y,normalize=False) likelihood.index = index - + """ + if mixed_noise_list == []: + likelihood = likelihoods.Gaussian_Mixed_Noise(Y_list,normalize=normalize_Y) + #TODO: allow passing the variance parameter into the model else: + self.likelihood_list = [] #TODO this is not needed assert len(Y_list) == len(mixed_noise_list) for noise,Y in zip(mixed_noise_list,Y_list): self.likelihood_list.append(likelihoods.EP(Y,noise)) + #TODO: allow normalization likelihood = likelihoods.EP_Mixed_Noise(Y_list, mixed_noise_list) - X = np.hstack([np.vstack(X_list),index]) original_dim = X.shape[1] - 1