The next step is to optimize the noise models' parameters

GPy/likelihoods/ep.py

@@ -65,11 +65,17 @@ class EP(likelihood):
         return self.noise_model.predictive_values(mu,var)
 
     def _get_params(self):
-        return np.zeros(0)
+        #return np.zeros(0)
+        return self.noise_model._get_params()
+
     def _get_param_names(self):
-        return []
+        #return []
+        return self.noise_model._get_param_names()
+
     def _set_params(self,p):
-        pass # TODO: the EP likelihood might want to take some parameters...
+        #pass # TODO: the EP likelihood might want to take some parameters...
+        self.noise_model._set_params(p)
+
     def _gradients(self,partial):
         return np.zeros(0) # TODO: the EP likelihood might want to take some parameters...
 

GPy/likelihoods/noise_model_constructors.py

@@ -3,8 +3,6 @@
 
 import numpy as np
 import noise_models
-#from likelihood_functions import LikelihoodFunction
-#import gp_transformations
 
 def binomial(gp_link=None):
     """
@@ -12,20 +10,32 @@ def binomial(gp_link=None):
 
     :param gp_link: a GPy gp_link function
     """
-    #self.discrete = True
-    #self.support_limits = (0,1)
-
     if gp_link is None:
         gp_link = noise_models.gp_transformations.Probit()
-    else:
-        assert isinstance(gp_link,noise_models.gp_transformations.GPTransformation), 'gp_link function is not valid.'
+    #else:
+    #    assert isinstance(gp_link,noise_models.gp_transformations.GPTransformation), 'gp_link function is not valid.'
 
     if isinstance(gp_link,noise_models.gp_transformations.Probit):
-        analytical_moments = True
+        analytical_mean = True
     else:
-        analytical_moments = False
-    return noise_models.binomial_noise.Binomial(gp_link,analytical_moments)
+        analytical_mean = False
+    analytical_variance = False
+    return noise_models.binomial_noise.Binomial(gp_link,analytical_mean,analytical_variance)
 
+def gaussian(gp_link=None,variance=1.):
+    """
+    Construct a gaussian likelihood
+
+    :param gp_link: a GPy gp_link function
+    """
+    if gp_link is None:
+        gp_link = noise_models.gp_transformations.Identity()
+    #else:
+    #    assert isinstance(gp_link,noise_models.gp_transformations.GPTransformation), 'gp_link function is not valid.'
+
+    analytical_mean = True
+    analytical_variance = True
+    return noise_models.gaussian_noise.Gaussian(gp_link,analytical_mean,analytical_variance,variance)
 
 def poisson(gp_link=None):
     """
@@ -35,8 +45,8 @@ def poisson(gp_link=None):
     """
     if gp_link is None:
         gp_link = noise_models.gp_transformations.Log_ex_1()
-    else:
-        assert isinstance(gp_link,noise_models.gp_transformations.GPTransformation), 'gp_link function is not valid.'
-    #assert isinstance(gp_link,gp_transformations.GPTransformation), 'gp_link function is not valid.'
-    analytical_moments = False
-    return noise_models.poisson_noise.Poisson(gp_link,analytical_moments)
+    #else:
+    #    assert isinstance(gp_link,noise_models.gp_transformations.GPTransformation), 'gp_link function is not valid.'
+    analytical_mean = False
+    analytical_variance = False
+    return noise_models.poisson_noise.Poisson(gp_link,analytical_mean,analytical_variance)

GPy/likelihoods/noise_models/__init__.py

@@ -1,4 +1,5 @@
 import noise_distributions
 import binomial_noise
+import gaussian_noise
 import poisson_noise
 import gp_transformations

GPy/likelihoods/noise_models/binomial_noise.py

@@ -17,8 +17,8 @@ class Binomial(NoiseDistribution):
     L(x) = \\Phi (Y_i*f_i)
     $$
     """
-    def __init__(self,gp_link=None,analytical_moments=False):
-        super(Binomial, self).__init__(gp_link,analytical_moments)
+    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
+        super(Binomial, self).__init__(gp_link,analytical_mean,analytical_variance)
 
     def _preprocess_values(self,Y):
         """

GPy/likelihoods/noise_models/gaussian_noise.py

@@ -15,10 +15,18 @@ class Gaussian(NoiseDistribution):
     :param mean: mean value of the Gaussian distribution
     :param variance: mean value of the Gaussian distribution
     """
-    def __init__(self,gp_link=None,analytical_moments=False,mean=0,variance=1.):
-        self.mean = mean
+    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False,variance=1.):
         self.variance = variance
-        super(Gaussian, self).__init__(gp_link,analytical_moments)
+        super(Gaussian, self).__init__(gp_link,analytical_mean,analytical_variance)
+
+    def _get_params(self):
+        return self.variance
+
+    def _get_param_names(self):
+        return ['noise_model_variance']
+
+    def _set_params(self,p):
+        self.variance = p
 
     def _preprocess_values(self,Y):
         """
@@ -36,32 +44,29 @@ class Gaussian(NoiseDistribution):
         :param v_i: mean/variance of the cavity distribution (float)
         """
         sigma2_hat = 1./(1./self.variance + tau_i)
-        mu_hat = sigma2_hat*(self.mean/self.variance + v_i)
-        Z_hat = np.sqrt(2*np.pi*sigma2_hat)*np.exp(-.5*(self.mean - v_i/tau_i)**2/(self.variance + 1./tau_i)) #TODO check
+        mu_hat = sigma2_hat*(data_i/self.variance + v_i)
+        sum_var = self.variance + 1./tau_i
+        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_mean_analytical(self,mu,sigma):
-        new_sigma2 = 1./(1./self.variance + 1./sigma**2)
-        return new_sigma2*(mu/sigma + self.mean/self.variance)
+        new_sigma2 = self.predictive_variance(mu,sigma)
+        return new_sigma2*(mu/sigma**2 + self.gp_link.transf(mu)/self.variance)
+
+    def _predictive_variance_analytical(self,mu,sigma,*args): #TODO *args?
+        return 1./(1./self.variance + 1./sigma**2)
 
     def _mass(self,gp,obs):
-        p = (self.gp_link.transf(gp)-self.mean)/np.sqrt(self.variance)
-        return std_norm_pdf(p)
+        return std_norm_pdf( (self.gp_link.transf(gp)-obs)/np.sqrt(self.variance) )
 
     def _nlog_mass(self,gp,obs):
-        p = (self.gp_link.transf(gp)-self.mean)/np.sqrt(self.variance)
-        return .5*np.log(2*np.pi*self.variance) + .5*(p-self.mean)**2/self.variance
+        return .5*((self.gp_link.transf(gp)-obs)**2/np.sqrt(self.variance) + np.log(2*np.pi*self.variance))
 
     def _dnlog_mass_dgp(self,gp,obs):
-        p = (self.gp_link.transf(gp)-self.mean)/np.sqrt(self.variance)
-        dp = self.gp_link.dtransf_df(gp)
-        return (p - self.mean)/self.variance * dp
+        return (self.gp_link.transf(gp)-obs)/np.sqrt(self.variance) * self.gp_link.dtransf_df(gp)
 
     def _d2nlog_mass_dgp2(self,gp,obs):
-        p = (self.gp_link.transf(gp)-self.mean)/np.sqrt(self.variance)
-        dp = self.gp_link.dtransf_df(gp)
-        d2p = self.gp_link.d2transf_df2(gp)
-        return dp**2/self.variance + (p - self.mean)/self.variance * d2p
+        return ((self.gp_link.transf(gp)-obs)*self.gp_link.d2transf_df2(gp) + self.gp_link.dtransf_df(gp)**2)/self.variance
 
     def _mean(self,gp):
         """
@@ -79,11 +84,10 @@ class Gaussian(NoiseDistribution):
         """
         Mass (or density) function
         """
-        p = self.gp_link.transf(gp)
-        return p*(1-p)
+        return self.variance
 
     def _dvariance_dgp(self,gp):
-        return self.gp_link.dtransf_df(gp)*(1. - 2.*self.gp_link.transf(gp))
+        return 0
 
     def _d2variance_dgp2(self,gp):
-        return self.gp_link.d2transf_df2(gp)*(1. - 2.*self.gp_link.transf(gp)) - 2*self.gp_link.dtransf_df(gp)**2
+        return 0

GPy/likelihoods/noise_models/noise_distributions.py

@@ -18,16 +18,30 @@ class NoiseDistribution(object):
     :param Y: observed output (Nx1 numpy.darray)
     ..Note:: Y values allowed depend on the LikelihoodFunction used
     """
-    def __init__(self,gp_link,analytical_moments=False):
+    def __init__(self,gp_link,analytical_mean=False,analytical_variance=False):
         #assert isinstance(gp_link,gp_transformations.GPTransformation), "gp_link is not a valid GPTransformation."#FIXME
         self.gp_link = gp_link
-        self.analytical_moments = analytical_moments
-        if self.analytical_moments:
+        self.analytical_mean = analytical_mean
+        self.analytical_variance = analytical_variance
+        if self.analytical_mean:
             self.moments_match = self._moments_match_analytical
             self.predictive_mean = self._predictive_mean_analytical
         else:
             self.moments_match = self._moments_match_numerical
             self.predictive_mean = self._predictive_mean_numerical
+        if self.analytical_variance:
+            self.predictive_variance = self._predictive_variance_analytical
+        else:
+            self.predictive_variance = self._predictive_variance_numerical
+
+    def _get_params(self):
+        return np.zeros(0)
+
+    def _get_param_names(self):
+        return []
+
+    def _set_params(self,p):
+        pass
 
     def _preprocess_values(self,Y):
         """
@@ -214,6 +228,12 @@ class NoiseDistribution(object):
         """
         pass
 
+    def _predictive_variance_analytical(self,mu,sigma):
+        """
+        If available, this function computes the predictive variance analytically.
+        """
+        pass
+
     def _predictive_mean_numerical(self,mu,sigma):
         """
         Laplace approximation to the predictive mean: E(Y_star) = E( E(Y_star|f_star) )
@@ -248,7 +268,7 @@ class NoiseDistribution(object):
         mean_squared = np.exp(-self._nlog_exp_conditional_mean_sq_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_mean_sq_dgp2(maximum,mu,sigma))*sigma)
         return mean_squared
 
-    def predictive_variance(self,mu,sigma,predictive_mean=None):
+    def _predictive_variance_numerical(self,mu,sigma,predictive_mean=None):
         """
         Laplace approximation to the predictive variance: V(Y_star) = E( V(Y_star|f_star) ) + V( E(Y_star|f_star) )
 

GPy/likelihoods/noise_models/poisson_noise.py

@@ -18,12 +18,12 @@ class Poisson(NoiseDistribution):
     L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
     $$
     """
-    def __init__(self,gp_link=None,analytical_moments=False):
+    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
         #self.discrete = True
         #self.support_limits = (0,np.inf)
 
-        #self.analytical_moments = False
-        super(Poisson, self).__init__(gp_link,analytical_moments)
+        #self.analytical_mean = False
+        super(Poisson, self).__init__(gp_link,analytical_mean,analytical_variance)
 
     def _preprocess_values(self,Y): #TODO
         #self.scale = .5*Y.max()