gamma noise added

GPy/likelihoods/noise_model_constructors.py

@@ -27,14 +27,15 @@ def gaussian(gp_link=None,variance=1.):
     Construct a gaussian likelihood
 
     :param gp_link: a GPy gp_link function
+    :param variance: scalar
     """
     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_mean = False
-    analytical_variance = True
+    analytical_variance = False
     return noise_models.gaussian_noise.Gaussian(gp_link,analytical_mean,analytical_variance,variance)
 
 def poisson(gp_link=None):
@@ -50,3 +51,18 @@ def poisson(gp_link=None):
     analytical_mean = False
     analytical_variance = False
     return noise_models.poisson_noise.Poisson(gp_link,analytical_mean,analytical_variance)
+
+def gamma(gp_link=None,beta=1.):
+    """
+    Construct a Gamma likelihood
+
+    :param gp_link: a GPy gp_link function
+    :param beta: scalar
+    """
+    if gp_link is None:
+        gp_link = noise_models.gp_transformations.Log_ex_1()
+    analytical_mean = False
+    analytical_variance = False
+    return noise_models.gamma_noise.Gamma(gp_link,analytical_mean,analytical_variance,beta)
+
+

GPy/likelihoods/noise_models/__init__.py

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

GPy/likelihoods/noise_models/binomial_noise.py

@@ -87,7 +87,7 @@ class Binomial(NoiseDistribution):
         Mass (or density) function
         """
         p = self.gp_link.transf(gp)
-        return p*(1-p)
+        return p*(1.-p)
 
     def _dvariance_dgp(self,gp):
         return self.gp_link.dtransf_df(gp)*(1. - 2.*self.gp_link.transf(gp))

GPy/likelihoods/noise_models/gamma_noise.py

@@ -0,0 +1,71 @@
+# Copyright (c) 2012, 2013 Ricardo Andrade
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+import numpy as np
+from scipy import stats,special
+import scipy as sp
+from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
+import gp_transformations
+from noise_distributions import NoiseDistribution
+
+class Gamma(NoiseDistribution):
+    """
+    Gamma likelihood
+    Y is expected to take values in {0,1,2,...}
+    -----
+    $$
+    L(x) = \exp(\lambda) * \lambda**Y_i / Y_i!
+    $$
+    """
+    def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False,beta=1.):
+        self.beta = beta
+        super(Gamma, self).__init__(gp_link,analytical_mean,analytical_variance)
+
+    def _preprocess_values(self,Y):
+        return Y
+
+    def _mass(self,gp,obs):
+        """
+        Mass (or density) function
+        """
+        #return stats.gamma.pdf(obs,a = self.gp_link.transf(gp)/self.variance,scale=self.variance)
+        alpha = self.gp_link.transf(gp)*self.beta
+        return obs**(alpha - 1.) * np.exp(-self.beta*obs) * self.beta**alpha / special.gamma(alpha)
+
+    def _nlog_mass(self,gp,obs):
+        """
+        Negative logarithm of the un-normalized distribution: factors that are not a function of gp are omitted
+        """
+        alpha = self.gp_link.transf(gp)*self.beta
+        return (1. - alpha)*np.log(obs) + self.beta*obs - alpha * np.log(self.beta) + np.log(special.gamma(alpha))
+
+    def _dnlog_mass_dgp(self,gp,obs):
+        return -self.gp_link.dtransf_df(gp)*self.beta*np.log(obs) + special.psi(self.gp_link.transf(gp)*self.beta) * self.gp_link.dtransf_df(gp)*self.beta
+
+    def _d2nlog_mass_dgp2(self,gp,obs):
+        return -self.gp_link.d2transf_df2(gp)*self.beta*np.log(obs) + special.polygamma(1,self.gp_link.transf(gp)*self.beta)*(self.gp_link.dtransf_df(gp)*self.beta)**2 + special.psi(self.gp_link.transf(gp)*self.beta)*self.gp_link.d2transf_df2(gp)*self.beta
+
+    def _mean(self,gp):
+        """
+        Mass (or density) function
+        """
+        return self.gp_link.transf(gp)
+
+    def _dmean_dgp(self,gp):
+        return self.gp_link.dtransf_df(gp)
+
+    def _d2mean_dgp2(self,gp):
+        return self.gp_link.d2transf_df2(gp)
+
+    def _variance(self,gp):
+        """
+        Mass (or density) function
+        """
+        return self.gp_link.transf(gp)/self.beta
+
+    def _dvariance_dgp(self,gp):
+        return self.gp_link.dtransf_df(gp)/self.beta
+
+    def _d2variance_dgp2(self,gp):
+        return self.gp_link.d2transf_df2(gp)/self.beta

GPy/likelihoods/noise_models/poisson_noise.py

@@ -66,9 +66,6 @@ class Poisson(NoiseDistribution):
         """
         return self.gp_link.transf(gp)
 
-    #def _variance(self,gp):
-    #    return self.gp_link.transf(gp)
-
     def _dmean_dgp(self,gp):
         return self.gp_link.dtransf_df(gp)
 
@@ -81,9 +78,6 @@ class Poisson(NoiseDistribution):
         """
         return self.gp_link.transf(gp)
 
-    #def _variance(self,gp):
-    #    return self.gp_link.transf(gp)
-
     def _dvariance_dgp(self,gp):
         return self.gp_link.dtransf_df(gp)