Added derivatives for poisson and a couple of examples,

need to fix for EP.

GPy/examples/regression.py

@@ -270,6 +270,50 @@ def toy_rbf_1d_50(max_iters=100):
     print(m)
     return m
 
+def toy_poisson_rbf_1d(optimizer='bfgs', max_nb_eval_optim=100):
+    """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
+    X = np.linspace(0,10)[:, None]
+    F = np.round(X*3-4)
+    F = np.where(F > 0, F, 0)
+    eps = np.random.randint(0,4, F.shape[0])[:, None]
+    Y = F + eps
+
+    noise_model = GPy.likelihoods.poisson()
+    likelihood = GPy.likelihoods.EP(Y,noise_model)
+
+    # create simple GP Model
+    m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
+
+    # optimize
+    m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
+    # plot
+    m.plot()
+    print(m)
+    return m
+
+def toy_poisson_rbf_1d_laplace(optimizer='bfgs', max_nb_eval_optim=100):
+    """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
+    X = np.linspace(0,10)[:, None]
+    F = np.round(X*3-4)
+    F = np.where(F > 0, F, 0)
+    eps = np.random.randint(0,4, F.shape[0])[:, None]
+    Y = F + eps
+
+    noise_model = GPy.likelihoods.poisson()
+    likelihood = GPy.likelihoods.Laplace(Y,noise_model)
+
+    # create simple GP Model
+    m = GPy.models.GPRegression(X, Y, likelihood=likelihood)
+
+    # optimize
+    m.optimize(optimizer, max_f_eval=max_nb_eval_optim)
+    # plot
+    m.plot()
+    print(m)
+    return m
+
+
+
 def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4):
     # Create an artificial dataset where the values in the targets (Y)
     # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to

GPy/likelihoods/noise_models/poisson_noise.py

@@ -1,7 +1,7 @@
+from __future__ import division
 # 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
@@ -14,9 +14,10 @@ class Poisson(NoiseDistribution):
     Poisson likelihood
 
     .. math::
-        L(x) = \\exp(\\lambda) * \\frac{\\lambda^Y_i}{Y_i!}
+        p(y_{i}|\\lambda(f_{i})) = \\frac{\\lambda(f_{i})^{y_{i}}}{y_{i}!}e^{-\\lambda(f_{i})}
 
-    ..Note: Y is expected to take values in {0,1,2,...}
+    .. Note::
+        Y is expected to take values in {0,1,2,...}
     """
     def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
         super(Poisson, self).__init__(gp_link,analytical_mean,analytical_variance)
@@ -24,25 +25,108 @@ class Poisson(NoiseDistribution):
     def _preprocess_values(self,Y): #TODO
         return Y
 
-    def _mass(self,gp,obs):
+    def pdf_link(self, link_f, y, extra_data=None):
         """
-        Mass (or density) function
+        Likelihood function given link(f)
-        """
-        return stats.poisson.pmf(obs,self.gp_link.transf(gp))
 
-    def _nlog_mass(self,gp,obs):
+        .. math::
-        """
+            p(y_{i}|\\lambda(f_{i})) = \\frac{\\lambda(f_{i})^{y_{i}}}{y_{i}!}e^{-\\lambda(f_{i})}
-        Negative logarithm of the un-normalized distribution: factors that are not a function of gp are omitted
-        """
-        return self.gp_link.transf(gp) - obs * np.log(self.gp_link.transf(gp)) + np.log(special.gamma(obs+1))
 
-    def _dnlog_mass_dgp(self,gp,obs):
+        :param link_f: latent variables link(f)
-        return self.gp_link.dtransf_df(gp) * (1. - obs/self.gp_link.transf(gp))
+        :type link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param extra_data: extra_data which is not used in poisson distribution
+        :returns: likelihood evaluated for this point
+        :rtype: float
+        """
+        assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
+        return np.prod(stats.poisson.pmf(y,link_f))
 
-    def _d2nlog_mass_dgp2(self,gp,obs):
+    def logpdf_link(self, link_f, y, extra_data=None):
-        d2_df = self.gp_link.d2transf_df2(gp)
+        """
-        transf = self.gp_link.transf(gp)
+        Log Likelihood Function given link(f)
-        return obs * ((self.gp_link.dtransf_df(gp)/transf)**2 - d2_df/transf) + d2_df
+
+        .. math::
+            \\ln p(y_{i}|\lambda(f_{i})) = -\\lambda(f_{i}) + y_{i}\\log \\lambda(f_{i}) - \\log y_{i}!
+
+        :param link_f: latent variables (link(f))
+        :type link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param extra_data: extra_data which is not used in poisson distribution
+        :returns: likelihood evaluated for this point
+        :rtype: float
+
+        """
+        assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
+        return np.sum(-link_f + y*np.log(link_f) - special.gammaln(y+1))
+
+    def dlogpdf_dlink(self, link_f, y, extra_data=None):
+        """
+        Gradient of the log likelihood function at y, given link(f) w.r.t link(f)
+
+        .. math::
+            \\frac{d \\ln p(y_{i}|\lambda(f_{i}))}{d\\lambda(f)} = \\frac{y_{i}}{\\lambda(f_{i})} - 1
+
+        :param link_f: latent variables (f)
+        :type link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param extra_data: extra_data which is not used in poisson distribution
+        :returns: gradient of likelihood evaluated at points
+        :rtype: Nx1 array
+
+        """
+        assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
+        return y/link_f - 1
+
+    def d2logpdf_dlink2(self, link_f, y, extra_data=None):
+        """
+        Hessian at y, given link(f), w.r.t link(f)
+        i.e. second derivative logpdf at y given link(f_i) and link(f_j)  w.r.t link(f_i) and link(f_j)
+        The hessian will be 0 unless i == j
+
+        .. math::
+            \\frac{d^{2} \\ln p(y_{i}|\lambda(f_{i}))}{d^{2}\\lambda(f)} = \\frac{-y_{i}}{\\lambda(f_{i})^{2}}
+
+        :param link_f: latent variables link(f)
+        :type link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param extra_data: extra_data which is not used in poisson distribution
+        :returns: Diagonal of hessian matrix (second derivative of likelihood evaluated at points f)
+        :rtype: Nx1 array
+
+        .. Note::
+            Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
+            (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i))
+        """
+        assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
+        hess = -y/(link_f**2)
+        return hess
+        #d2_df = self.gp_link.d2transf_df2(gp)
+        #transf = self.gp_link.transf(gp)
+        #return obs * ((self.gp_link.dtransf_df(gp)/transf)**2 - d2_df/transf) + d2_df
+
+    def d3logpdf_dlink3(self, link_f, y, extra_data=None):
+        """
+        Third order derivative log-likelihood function at y given link(f) w.r.t link(f)
+
+        .. math::
+            \\frac{d^{3} \\ln p(y_{i}|\lambda(f_{i}))}{d^{3}\\lambda(f)} = \\frac{2y_{i}}{\\lambda(f_{i})^{3}}
+
+        :param link_f: latent variables link(f)
+        :type link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param extra_data: extra_data which is not used in poisson distribution
+        :returns: third derivative of likelihood evaluated at points f
+        :rtype: Nx1 array
+        """
+        assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
+        d3lik_dlink3 = 2*y/(link_f)**3
+        return d3lik_dlink3
 
     def _mean(self,gp):
         """
@@ -55,3 +139,15 @@ class Poisson(NoiseDistribution):
         Mass (or density) function
         """
         return self.gp_link.transf(gp)
+
+    def samples(self, gp):
+        """
+        Returns a set of samples of observations based on a given value of the latent variable.
+
+        :param size: number of samples to compute
+        :param gp: latent variable
+        """
+        orig_shape = gp.shape
+        gp = gp.flatten()
+        Ysim = np.array([np.random.poisson(self.gp_link.transf(gpj),size=1) for gpj in gp])
+        return Ysim.reshape(orig_shape)

GPy/testing/likelihoods_tests.py

@@ -84,6 +84,10 @@ class TestNoiseModels(object):
         self.f = np.random.rand(self.N, 1)
         self.binary_Y = np.asarray(np.random.rand(self.N) > 0.5, dtype=np.int)[:, None]
         self.positive_Y = np.exp(self.Y.copy())
+        self.integer_Y = np.round(self.X[:, 0]*3-3)[:, None] + np.random.randint(0,3, self.X.shape[0])[:, None]
+        self.integer_Y = np.where(self.integer_Y > 0, self.integer_Y, 0)
+        print self.integer_Y
+        print self.Y
 
         self.var = 0.2
 
@@ -223,6 +227,13 @@ class TestNoiseModels(object):
                             "link_f_constraints": [constrain_positive],
                             "Y": self.positive_Y,
                             "laplace": True,
+                        },
+                        "Poisson_default": {
+                            "model": GPy.likelihoods.poisson(),
+                            "link_f_constraints": [constrain_positive],
+                            "Y": self.integer_Y,
+                            "laplace": True,
+                            "ep": False #Should work though...
                         }
                     }