From 68e37e8684703fe26f71ef693547b525c84c9421 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 10 Jul 2013 19:33:43 +0100 Subject: [PATCH] The next step is to optimize the noise models' parameters --- GPy/likelihoods/ep.py | 12 +++-- GPy/likelihoods/noise_model_constructors.py | 40 ++++++++++------ GPy/likelihoods/noise_models/__init__.py | 1 + .../noise_models/binomial_noise.py | 4 +- .../noise_models/gaussian_noise.py | 48 ++++++++++--------- .../noise_models/noise_distributions.py | 28 +++++++++-- GPy/likelihoods/noise_models/poisson_noise.py | 6 +-- 7 files changed, 90 insertions(+), 49 deletions(-) diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py index 717bfcb7..50df599a 100644 --- a/GPy/likelihoods/ep.py +++ b/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... diff --git a/GPy/likelihoods/noise_model_constructors.py b/GPy/likelihoods/noise_model_constructors.py index 4267fc32..f64fb235 100644 --- a/GPy/likelihoods/noise_model_constructors.py +++ b/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) diff --git a/GPy/likelihoods/noise_models/__init__.py b/GPy/likelihoods/noise_models/__init__.py index c5fc66b0..b16f8fbd 100644 --- a/GPy/likelihoods/noise_models/__init__.py +++ b/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 diff --git a/GPy/likelihoods/noise_models/binomial_noise.py b/GPy/likelihoods/noise_models/binomial_noise.py index 77fde5eb..200625d9 100644 --- a/GPy/likelihoods/noise_models/binomial_noise.py +++ b/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): """ diff --git a/GPy/likelihoods/noise_models/gaussian_noise.py b/GPy/likelihoods/noise_models/gaussian_noise.py index a77becb2..389363a3 100644 --- a/GPy/likelihoods/noise_models/gaussian_noise.py +++ b/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 diff --git a/GPy/likelihoods/noise_models/noise_distributions.py b/GPy/likelihoods/noise_models/noise_distributions.py index 0dc9e03c..45f92950 100644 --- a/GPy/likelihoods/noise_models/noise_distributions.py +++ b/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) ) diff --git a/GPy/likelihoods/noise_models/poisson_noise.py b/GPy/likelihoods/noise_models/poisson_noise.py index e90b3ce8..f7dd93b9 100644 --- a/GPy/likelihoods/noise_models/poisson_noise.py +++ b/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()