From ddf64629ae9977c13569c6e608d5f7d45862929c Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 16 Jul 2013 18:38:13 +0100 Subject: [PATCH] gamma noise added --- GPy/likelihoods/noise_model_constructors.py | 20 +++++- GPy/likelihoods/noise_models/__init__.py | 1 + .../noise_models/binomial_noise.py | 2 +- GPy/likelihoods/noise_models/gamma_noise.py | 71 +++++++++++++++++++ GPy/likelihoods/noise_models/poisson_noise.py | 6 -- 5 files changed, 91 insertions(+), 9 deletions(-) create mode 100644 GPy/likelihoods/noise_models/gamma_noise.py diff --git a/GPy/likelihoods/noise_model_constructors.py b/GPy/likelihoods/noise_model_constructors.py index f64fb235..cc205c6d 100644 --- a/GPy/likelihoods/noise_model_constructors.py +++ b/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_variance = True + analytical_mean = False + 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) + + diff --git a/GPy/likelihoods/noise_models/__init__.py b/GPy/likelihoods/noise_models/__init__.py index b16f8fbd..65a94e1e 100644 --- a/GPy/likelihoods/noise_models/__init__.py +++ b/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 diff --git a/GPy/likelihoods/noise_models/binomial_noise.py b/GPy/likelihoods/noise_models/binomial_noise.py index 200625d9..e47d9211 100644 --- a/GPy/likelihoods/noise_models/binomial_noise.py +++ b/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)) diff --git a/GPy/likelihoods/noise_models/gamma_noise.py b/GPy/likelihoods/noise_models/gamma_noise.py new file mode 100644 index 00000000..6bf0dd7b --- /dev/null +++ b/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 diff --git a/GPy/likelihoods/noise_models/poisson_noise.py b/GPy/likelihoods/noise_models/poisson_noise.py index f7dd93b9..e4ce90d3 100644 --- a/GPy/likelihoods/noise_models/poisson_noise.py +++ b/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)