lots of fixes, including prediction being mean and variance only

GPy/likelihoods/gaussian.py

@@ -18,6 +18,7 @@ import link_functions
 from likelihood import Likelihood
 from ..core.parameterization import Param
 from ..core.parameterization.transformations import Logexp
+from scipy import stats
 
 class Gaussian(Likelihood):
     """
@@ -49,11 +50,14 @@ class Gaussian(Likelihood):
         if isinstance(gp_link, link_functions.Identity):
             self.log_concave = True
 
-    def covariance_matrix(self, Y, **Y_metadata):
+    def covariance_matrix(self, Y, Y_metadata=None):
         return np.eye(Y.shape[0]) * self.variance
 
-    def update_gradients(self, partial):
-        self.variance.gradient = np.sum(partial)
+    def update_gradients(self, grad):
+        self.variance.gradient = grad
+
+    def exact_inference_gradients(self, dL_dKdiag):
+        return dL_dKdiag.sum()
 
     def _preprocess_values(self, Y):
         """
@@ -76,16 +80,12 @@ class Gaussian(Likelihood):
         Z_hat = 1./np.sqrt(2.*np.pi*sum_var)*np.exp(-.5*(data_i - v_i/tau_i)**2./sum_var)
         return Z_hat, mu_hat, sigma2_hat
 
-    def predictive_values(self, mu, var, full_cov=False):
+    def predictive_values(self, mu, var, full_cov=False, Y_metadata=None):
         if full_cov:
             var += np.eye(var.shape[0])*self.variance
-            d = 2*np.sqrt(np.diag(var))
-            low, up = mu - d, mu + d
         else:
             var += self.variance
-            d = 2*np.sqrt(var)
-            low, up = mu - d, mu + d
-        return mu, var, low, up
+        return mu, var
 
     def predictive_mean(self, mu, sigma):
         return mu
@@ -93,6 +93,9 @@ class Gaussian(Likelihood):
     def predictive_variance(self, mu, sigma, predictive_mean=None):
         return self.variance + sigma**2
 
+    def predictive_quantiles(self, mu, var, quantiles, Y_metadata):
+        return  [stats.norm.ppf(q)*np.sqrt(var) + mu for q in quantiles]
+
     def pdf_link(self, link_f, y, extra_data=None):
         """
         Likelihood function given link(f)

GPy/likelihoods/likelihood.py

@@ -135,7 +135,7 @@ class Likelihood(Parameterized):
 
         return mean
 
-    def _predictive_variance(self,mu,variance,predictive_mean=None):
+    def _predictive_variance(self, mu,variance, predictive_mean=None):
         """
         Numerical approximation to the predictive variance: V(Y_star)
 
@@ -358,7 +358,7 @@ class Likelihood(Parameterized):
 
         return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta
 
-    def predictive_values(self, mu, var, full_cov=False, sampling=True, num_samples=10000):
+    def predictive_values(self, mu, var, full_cov=False, Y_metadata=None):
         """
         Compute  mean, variance and conficence interval (percentiles 5 and 95) of the  prediction.
 
@@ -366,14 +366,21 @@ class Likelihood(Parameterized):
         :param var: variance of the latent variable, f, of posterior
         :param full_cov: whether to use the full covariance or just the diagonal
         :type full_cov: Boolean
-        :param num_samples: number of samples to use in computing quantiles and
-                            possibly mean variance
-        :type num_samples: integer
-        :param sampling: Whether to use samples for mean and variances anyway
-        :type sampling: Boolean
-
         """
 
+        pred_mean = self.predictive_mean(mu, var, Y_metadata)
+        pred_var = self.predictive_variance(mu, var, pred_mean, Y_metadata)
+
+        return pred_mean, pred_var
+
+
+    def samples(self, gp):
+        """
+        Returns a set of samples of observations based on a given value of the latent variable.
+
+        :param gp: latent variable
+        """
+        raise NotImplementedError
         if sampling:
             #Get gp_samples f* using posterior mean and variance
             if not full_cov:
@@ -393,20 +400,4 @@ class Likelihood(Parameterized):
             q1 = np.percentile(samples, 2.5, axis=axis)[:,None]
             q3 = np.percentile(samples, 97.5, axis=axis)[:,None]
 
-        else:
 
-            pred_mean = self.predictive_mean(mu, var)
-            pred_var = self.predictive_variance(mu, var, pred_mean)
-            print "WARNING: Predictive quantiles are only computed when sampling."
-            q1 = np.repeat(np.nan,pred_mean.size)[:,None]
-            q3 = q1.copy()
-
-        return pred_mean, pred_var, q1, q3
-
-    def samples(self, gp):
-        """
-        Returns a set of samples of observations based on a given value of the latent variable.
-
-        :param gp: latent variable
-        """
-        raise NotImplementedError

GPy/likelihoods/mixed_noise.py

@@ -3,56 +3,57 @@ 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 gaussian import Gaussian
 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
+    def __init__(self, likelihoods_list, name='mixed_noise'):
 
         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, gradients):
+        self.gradient = gradients
 
-    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 exact_inference_gradients(self, dL_dKdiag, Y_metadata):
+        assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
+        ind = Y_metadata['output_index']
+        return np.array([dL_dKdiag[ind==i].sum() for i in range(len(self.likelihoods_list))])
 
-    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
+    def predictive_values(self, mu, var, full_cov=False, Y_metadata=None):
+        if all([isinstance(l, Gaussian) for l in self.likelihoods_list]):
+            ind = Y_metadata['output_index']
+            _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
         else:
-            var += _variance
-            d = 2*np.sqrt(var)
-            low, up = mu - d, mu + d
-        return mu, var, low, up
+            raise NotImplementedError
 
-    def predictive_variance(self, mu, sigma, noise_index, predictive_mean=None, **Y_metadata):
+    def predictive_variance(self, mu, sigma, **other_shit):
         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
+
+
+    def covariance_matrix(self, Y, Y_metadata):
+        assert all([isinstance(l, Gaussian) for l in self.likelihoods_list])
+        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)
+