GPCoregionalizedRegresssion added

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

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

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
+from gp_coregionalized_regression import GPCoregionalizedRegression
+from sparse_gp_coregionalized_regression import SparseGPCoregionalizedRegression

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)