New Gaussian likelihood for multiple outputs

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):

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)

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])
 

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
 

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])
+
+
+
+

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