From 2fdb60287f768db6e08ae3c515ad711cf5f61376 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Fri, 25 Oct 2013 15:08:53 +0100 Subject: [PATCH] Added derivatives for poisson and a couple of examples, need to fix for EP. --- GPy/examples/regression.py | 44 ++++++ GPy/likelihoods/noise_models/poisson_noise.py | 132 +++++++++++++++--- GPy/testing/likelihoods_tests.py | 11 ++ 3 files changed, 169 insertions(+), 18 deletions(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index ca4f506d..2978ebdc 100644 --- a/GPy/examples/regression.py +++ b/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 diff --git a/GPy/likelihoods/noise_models/poisson_noise.py b/GPy/likelihoods/noise_models/poisson_noise.py index 80d7951b..fba00417 100644 --- a/GPy/likelihoods/noise_models/poisson_noise.py +++ b/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 - """ - return stats.poisson.pmf(obs,self.gp_link.transf(gp)) + Likelihood function given link(f) - def _nlog_mass(self,gp,obs): - """ - 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)) + .. math:: + p(y_{i}|\\lambda(f_{i})) = \\frac{\\lambda(f_{i})^{y_{i}}}{y_{i}!}e^{-\\lambda(f_{i})} - def _dnlog_mass_dgp(self,gp,obs): - return self.gp_link.dtransf_df(gp) * (1. - obs/self.gp_link.transf(gp)) + :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.prod(stats.poisson.pmf(y,link_f)) - def _d2nlog_mass_dgp2(self,gp,obs): - 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 logpdf_link(self, link_f, y, extra_data=None): + """ + Log Likelihood Function given link(f) + + .. 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) diff --git a/GPy/testing/likelihoods_tests.py b/GPy/testing/likelihoods_tests.py index c3ea6a43..155842fd 100644 --- a/GPy/testing/likelihoods_tests.py +++ b/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... } }