From 7361d311c143e2057299297cda33fd0c18f488e0 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Sun, 30 Jun 2013 23:36:37 +0100 Subject: [PATCH] further corrections --- GPy/likelihoods/binomial_likelihood.py | 62 ++++++++++++++++++++----- GPy/likelihoods/likelihood_functions.py | 11 +++-- GPy/likelihoods/poisson_likelihood.py | 18 ++++--- GPy/models/gp_classification.py | 3 +- 4 files changed, 68 insertions(+), 26 deletions(-) diff --git a/GPy/likelihoods/binomial_likelihood.py b/GPy/likelihoods/binomial_likelihood.py index 420a9607..15b8067e 100644 --- a/GPy/likelihoods/binomial_likelihood.py +++ b/GPy/likelihoods/binomial_likelihood.py @@ -23,17 +23,15 @@ class Binomial(LikelihoodFunction): def __init__(self,link=None): self.discrete = True self.support_limits = (0,1) - self._analytical = link_functions.Probit + if not link: - link = self._analytical + link = link_functions.Probit + if isinstance(link,link_functions.Probit): + self.analytical_moments = True + else: + self.analytical_moments = False super(Binomial, self).__init__(link) - def _mass(self,gp,obs): - pass - - def _nlog_mass(self,gp,obs): - pass - def _preprocess_values(self,Y): """ Check if the values of the observations correspond to the values @@ -66,12 +64,51 @@ class Binomial(LikelihoodFunction): def _predictive_mean_analytical(self,mu,sigma): return stats.norm.cdf(mu/np.sqrt(1+sigma**2)) - def predictive_values(self,mu,var): + def _mass(self,gp,obs): + #NOTE obs must be in {0,1} + p = self.link.inv_transf(gp) + return p**obs * (1.-p)**(1.-obs) + + def _nlog_mass(self,gp,obs): + p = self.link.inv_transf(gp) + return obs*np.log(p) + (1.-obs)*np.log(1-p) + + def _dnlog_mass_dgp(self,gp,obs): + p = self.link.inv_transf(gp) + dp = self.link.dinv_transf_df(gp) + return obs/p * dp - (1.-obs)/(1.-p) * dp + + def _d2nlog_mass_dgp2(self,gp,obs): + p = self.link.inv_transf(gp) + return (obs/p + (1.-obs)/(1.-p))*self.lind.d2inv_transf_df(gp) + ((1.-obs)/(1.-p)**2-obs/p**2)*self.link.dinv_transf_df(gp) + + def _mean(self,gp): """ - Compute mean, variance and conficence interval (percentiles 5 and 95) of the prediction - :param mu: mean of the latent variable - :param var: variance of the latent variable + Mass (or density) function """ + return self.link.inv_transf(gp) + + def _dmean_dgp(self,gp): + return self.link.dinv_transf_df(gp) + + def _d2mean_dgp2(self,gp): + return self.link.d2inv_transf_df2(gp) + + def _variance(self,gp): + """ + Mass (or density) function + """ + p = self.link.inv_transf(gp) + return p*(1-p) + + def _dvariance_dgp(self,gp): + return self.link.dinv_transf_df(gp)*(1. - 2.*self.link.inv_transf(gp)) + + def _d2variance_dgp2(self,gp): + return self.link.d2inv_transf_df2(gp)*(1. - 2.*self.link.inv_transf(gp)) - 2*self.link.dinv_transf_df(gp)**2 + + """ + def predictive_values(self,mu,var): #TODO remove mu = mu.flatten() var = var.flatten() #mean = stats.norm.cdf(mu/np.sqrt(1+var)) @@ -83,3 +120,4 @@ class Binomial(LikelihoodFunction): p_025 = self._predictive_mean_analytical(norm_025,np.sqrt(var)) p_975 = self._predictive_mean_analytical(norm_975,np.sqrt(var)) return mean[:,None], np.nan*var, p_025[:,None], p_975[:,None] # TODO: var + """ diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py index cb0be86a..9b132d72 100644 --- a/GPy/likelihoods/likelihood_functions.py +++ b/GPy/likelihoods/likelihood_functions.py @@ -19,12 +19,13 @@ class LikelihoodFunction(object): ..Note:: Y values allowed depend on the LikelihoodFunction used """ def __init__(self,link): - if link == self._analytical: + #assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME + self.link = link + + if self.analytical_moments: self.moments_match = self._moments_match_analytical self.predictive_mean = self._predictive_mean_analytical else: - assert isinstance(link,link_functions.LinkFunction) - self.link = link self.moments_match = self._moments_match_numerical self.predictive_mean = self._predictive_mean_numerical @@ -258,6 +259,7 @@ class LikelihoodFunction(object): maximum = sp.optimize.fmin_ncg(self._nlog_exp_conditional_variance_scaled,x0=self._variance(mu),fprime=self._dnlog_exp_conditional_variance_dgp,fhess=self._d2nlog_exp_conditional_variance_dgp2,args=(mu,sigma)) exp_var = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))/(np.sqrt(self._d2nlog_exp_conditional_variance_dgp2(maximum,mu,sigma))*sigma) + """ pb.figure() x = np.array([mu + step*sigma for step in np.linspace(-7,7,100)]) f = np.array([np.exp(-self._nlog_exp_conditional_variance_scaled(xi,mu,sigma))/np.sqrt(2*np.pi*sigma**2) for xi in x]) @@ -267,6 +269,7 @@ class LikelihoodFunction(object): k = np.exp(-self._nlog_exp_conditional_variance_scaled(maximum,mu,sigma))*np.sqrt(sigma2)/np.sqrt(sigma**2) pb.plot(x,f2*exp_var,'r--') pb.vlines(maximum,0,f.max()) + """ #V( E(Y_star|f_star) ) = E( E(Y_star|f_star)**2 ) - E( E(Y_star|f_star)**2 ) exp_exp2 = self._predictive_mean_sq(mu,sigma) @@ -323,6 +326,8 @@ class LikelihoodFunction(object): def predictive_values(self,mu,var,sample=True,sample_size=5000): """ Compute mean, variance and conficence interval (percentiles 5 and 95) of the prediction + :param mu: mean of the latent variable + :param var: variance of the latent variable """ if isinstance(mu,float) or isinstance(mu,int): mu = [mu] diff --git a/GPy/likelihoods/poisson_likelihood.py b/GPy/likelihoods/poisson_likelihood.py index 9466d4de..fba89331 100644 --- a/GPy/likelihoods/poisson_likelihood.py +++ b/GPy/likelihoods/poisson_likelihood.py @@ -23,9 +23,8 @@ class Poisson(LikelihoodFunction): def __init__(self,link=None): self.discrete = True self.support_limits = (0,np.inf) - self._analytical = None - if not link: - link = link_functions.Log() + + self.analytical_moments = False super(Poisson, self).__init__(link) def _mass(self,gp,obs): @@ -34,15 +33,14 @@ class Poisson(LikelihoodFunction): """ return stats.poisson.pmf(obs,self.link.inv_transf(gp)) - def _percentile(self,x,gp,*args): #TODO *args - return stats.poisson.ppf(x,self.link.inv_transf(gp)) - 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.link.inv_transf(gp) - obs * np.log(self.link.inv_transf(gp)) + np.log(special.gamma(obs+1)) + #def _preprocess_values(self,Y): #TODO + def _dnlog_mass_dgp(self,gp,obs): return self.link.dinv_transf_df(gp) * (1. - obs/self.link.inv_transf(gp)) @@ -66,8 +64,8 @@ class Poisson(LikelihoodFunction): """ return self.link.inv_transf(gp) - def _variance(self,gp): - return self.link.inv_transf(gp) + #def _variance(self,gp): + # return self.link.inv_transf(gp) def _dmean_dgp(self,gp): return self.link.dinv_transf_df(gp) @@ -81,8 +79,8 @@ class Poisson(LikelihoodFunction): """ return self.link.inv_transf(gp) - def _variance(self,gp): - return self.link.inv_transf(gp) + #def _variance(self,gp): + # return self.link.inv_transf(gp) def _dvariance_dgp(self,gp): return self.link.dinv_transf_df(gp) diff --git a/GPy/models/gp_classification.py b/GPy/models/gp_classification.py index c6012988..2e0d9c4a 100644 --- a/GPy/models/gp_classification.py +++ b/GPy/models/gp_classification.py @@ -31,7 +31,8 @@ class GPClassification(GP): kernel = kern.rbf(X.shape[1]) if likelihood is None: - distribution = likelihoods.likelihood_functions.Binomial() + #distribution = GPy.likelihoods.binomial_likelihood.Binomial(link=link) + distribution = likelihoods.binomial_likelihood.Binomial() likelihood = likelihoods.EP(Y, distribution) elif Y is not None: if not all(Y.flatten() == likelihood.data.flatten()):