From 06ffb884abd7b1eb135f0b474f90c981ccc4c9c4 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Tue, 9 Jul 2013 17:54:56 +0100 Subject: [PATCH] files re-organized --- GPy/likelihoods/__init__.py | 2 +- GPy/likelihoods/constructors.py | 42 --------- GPy/likelihoods/ep.py | 28 +++--- GPy/likelihoods/noise_model_constructors.py | 42 +++++++++ GPy/likelihoods/noise_models/__init__.py | 8 +- ...nomial_likelihood.py => binomial_noise.py} | 34 +++---- .../noise_models/gaussian_noise.py | 89 +++++++++++++++++++ ...ink_functions.py => gp_transformations.py} | 56 ++++++------ ...od_functions.py => noise_distributions.py} | 14 +-- ...poisson_likelihood.py => poisson_noise.py} | 43 +++++---- 10 files changed, 223 insertions(+), 135 deletions(-) delete mode 100644 GPy/likelihoods/constructors.py create mode 100644 GPy/likelihoods/noise_model_constructors.py rename GPy/likelihoods/noise_models/{binomial_likelihood.py => binomial_noise.py} (76%) create mode 100644 GPy/likelihoods/noise_models/gaussian_noise.py rename GPy/likelihoods/noise_models/{link_functions.py => gp_transformations.py} (59%) rename GPy/likelihoods/noise_models/{likelihood_functions.py => noise_distributions.py} (97%) rename GPy/likelihoods/noise_models/{poisson_likelihood.py => poisson_noise.py} (58%) diff --git a/GPy/likelihoods/__init__.py b/GPy/likelihoods/__init__.py index 4932bd40..3e6a28d3 100644 --- a/GPy/likelihoods/__init__.py +++ b/GPy/likelihoods/__init__.py @@ -1,5 +1,5 @@ from ep import EP from gaussian import Gaussian +from noise_model_constructors import * # TODO: from Laplace import Laplace -from constructors import * diff --git a/GPy/likelihoods/constructors.py b/GPy/likelihoods/constructors.py deleted file mode 100644 index 0b995894..00000000 --- a/GPy/likelihoods/constructors.py +++ /dev/null @@ -1,42 +0,0 @@ -# Copyright (c) 2013, GPy authors (see AUTHORS.txt). -# Licensed under the BSD 3-clause license (see LICENSE.txt) - -import numpy as np -from likelihood_functions import LikelihoodFunction -import noise_models -import link_functions - -def binomial(link=None): - """ - Construct a binomial likelihood - - :param link: a GPy link function - """ - #self.discrete = True - #self.support_limits = (0,1) - - if link is None: - link = link_functions.Probit() - else: - assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.' - - if isinstance(link,link_functions.Probit): - analytical_moments = True - else: - analytical_moments = False - return noise_models.binomial_likelihood.Binomial(link,analytical_moments) - - -def poisson(link=None): - """ - Construct a Poisson likelihood - - :param link: a GPy link function - """ - if link is None: - link = link_functions.Log_ex_1() - else: - assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.' - #assert isinstance(link,link_functions.LinkFunction), 'link function is not valid.' - analytical_moments = False - return noise_models.poisson_likelihood.Poisson(link,analytical_moments) diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py index 7e90755e..717bfcb7 100644 --- a/GPy/likelihoods/ep.py +++ b/GPy/likelihoods/ep.py @@ -4,23 +4,23 @@ from ..util.linalg import pdinv,mdot,jitchol,chol_inv,DSYR,tdot,dtrtrs from likelihood import likelihood class EP(likelihood): - def __init__(self,data,LikelihoodFunction,epsilon=1e-3,power_ep=[1.,1.]): + def __init__(self,data,noise_model,epsilon=1e-3,power_ep=[1.,1.]): """ Expectation Propagation Arguments --------- epsilon : Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float) - LikelihoodFunction : a likelihood function (see likelihood_functions.py) + noise_model : a likelihood function (see likelihood_functions.py) """ - self.LikelihoodFunction = LikelihoodFunction + self.noise_model = noise_model self.epsilon = epsilon self.eta, self.delta = power_ep self.data = data self.N, self.output_dim = self.data.shape self.is_heteroscedastic = True self.Nparams = 0 - self._transf_data = self.LikelihoodFunction._preprocess_values(data) + self._transf_data = self.noise_model._preprocess_values(data) #Initial values - Likelihood approximation parameters: #p(y|f) = t(f|tau_tilde,v_tilde) @@ -28,9 +28,9 @@ class EP(likelihood): self.tau_tilde = np.zeros(self.N) self.v_tilde = np.zeros(self.N) - #_gp = self.LikelihoodFunction.link.transf(self.data) - #_mean = self.LikelihoodFunction._mean(_gp) - #_variance = self.LikelihoodFunction._variance(_gp) + #_gp = self.noise_model.gp_link.transf(self.data) + #_mean = self.noise_model._mean(_gp) + #_variance = self.noise_model._variance(_gp) #self.tau_tilde = 1./_variance #self.tau_tilde[_variance== 0] = 1. #self.v_tilde = _mean*self.tau_tilde @@ -62,7 +62,7 @@ class EP(likelihood): def predictive_values(self,mu,var,full_cov): if full_cov: raise NotImplementedError, "Cannot make correlated predictions with an EP likelihood" - return self.LikelihoodFunction.predictive_values(mu,var) + return self.noise_model.predictive_values(mu,var) def _get_params(self): return np.zeros(0) @@ -128,7 +128,7 @@ class EP(likelihood): self.tau_[i] = 1./Sigma[i,i] - self.eta*self.tau_tilde[i] self.v_[i] = mu[i]/Sigma[i,i] - self.eta*self.v_tilde[i] #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.LikelihoodFunction.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) + self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) #DELETE @@ -136,11 +136,11 @@ class EP(likelihood): import pylab as pb from scipy import stats import scipy as sp - import link_functions + import gp_transformations from constructors import * - link = link_functions.Log_ex_1() - distribution = poisson(link=link) + gp_link = gp_transformations.Log_ex_1() + distribution = poisson(gp_link=gp_link) gp = np.linspace(-3,50,100) #distribution = binomial() #gp = np.linspace(-3,3,100) @@ -311,7 +311,7 @@ class EP(likelihood): self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i] self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i] #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.LikelihoodFunction.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) + self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) #Site parameters update Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i]) Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i]) @@ -406,7 +406,7 @@ class EP(likelihood): self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i] self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i] #Marginal moments - self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.LikelihoodFunction.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) + self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i]) #Site parameters update Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i]) Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i]) diff --git a/GPy/likelihoods/noise_model_constructors.py b/GPy/likelihoods/noise_model_constructors.py new file mode 100644 index 00000000..4267fc32 --- /dev/null +++ b/GPy/likelihoods/noise_model_constructors.py @@ -0,0 +1,42 @@ +# Copyright (c) 2013, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +import numpy as np +import noise_models +#from likelihood_functions import LikelihoodFunction +#import gp_transformations + +def binomial(gp_link=None): + """ + Construct a binomial likelihood + + :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.' + + if isinstance(gp_link,noise_models.gp_transformations.Probit): + analytical_moments = True + else: + analytical_moments = False + return noise_models.binomial_noise.Binomial(gp_link,analytical_moments) + + +def poisson(gp_link=None): + """ + Construct a Poisson likelihood + + :param gp_link: a GPy gp_link function + """ + 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) diff --git a/GPy/likelihoods/noise_models/__init__.py b/GPy/likelihoods/noise_models/__init__.py index 40282402..c5fc66b0 100644 --- a/GPy/likelihoods/noise_models/__init__.py +++ b/GPy/likelihoods/noise_models/__init__.py @@ -1,4 +1,4 @@ -import likelihood_functions -import binomial_likelihood -import poisson_likelihood -import link_functions +import noise_distributions +import binomial_noise +import poisson_noise +import gp_transformations diff --git a/GPy/likelihoods/noise_models/binomial_likelihood.py b/GPy/likelihoods/noise_models/binomial_noise.py similarity index 76% rename from GPy/likelihoods/noise_models/binomial_likelihood.py rename to GPy/likelihoods/noise_models/binomial_noise.py index d23dd2f7..77fde5eb 100644 --- a/GPy/likelihoods/noise_models/binomial_likelihood.py +++ b/GPy/likelihoods/noise_models/binomial_noise.py @@ -5,10 +5,10 @@ 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 link_functions -from likelihood_functions import NoiseModel +import gp_transformations +from noise_distributions import NoiseDistribution -class Binomial(NoiseModel): +class Binomial(NoiseDistribution): """ Probit likelihood Y is expected to take values in {-1,1} @@ -17,8 +17,8 @@ class Binomial(NoiseModel): L(x) = \\Phi (Y_i*f_i) $$ """ - def __init__(self,link=None,analytical_moments=False): - super(Binomial, self).__init__(link,analytical_moments) + def __init__(self,gp_link=None,analytical_moments=False): + super(Binomial, self).__init__(gp_link,analytical_moments) def _preprocess_values(self,Y): """ @@ -54,46 +54,46 @@ class Binomial(NoiseModel): def _mass(self,gp,obs): #NOTE obs must be in {0,1} - p = self.link.inv_transf(gp) + p = self.gp_link.transf(gp) return p**obs * (1.-p)**(1.-obs) def _nlog_mass(self,gp,obs): - p = self.link.inv_transf(gp) + p = self.gp_link.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) + p = self.gp_link.transf(gp) + dp = self.gp_link.dtransf_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.link.d2inv_transf_df2(gp) + ((1.-obs)/(1.-p)**2-obs/p**2)*self.link.dinv_transf_df(gp) + p = self.gp_link.transf(gp) + return (obs/p + (1.-obs)/(1.-p))*self.gp_link.d2transf_df2(gp) + ((1.-obs)/(1.-p)**2-obs/p**2)*self.gp_link.dtransf_df(gp) def _mean(self,gp): """ Mass (or density) function """ - return self.link.inv_transf(gp) + return self.gp_link.transf(gp) def _dmean_dgp(self,gp): - return self.link.dinv_transf_df(gp) + return self.gp_link.dtransf_df(gp) def _d2mean_dgp2(self,gp): - return self.link.d2inv_transf_df2(gp) + return self.gp_link.d2transf_df2(gp) def _variance(self,gp): """ Mass (or density) function """ - p = self.link.inv_transf(gp) + p = self.gp_link.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)) + return self.gp_link.dtransf_df(gp)*(1. - 2.*self.gp_link.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 + return self.gp_link.d2transf_df2(gp)*(1. - 2.*self.gp_link.transf(gp)) - 2*self.gp_link.dtransf_df(gp)**2 """ def predictive_values(self,mu,var): #TODO remove diff --git a/GPy/likelihoods/noise_models/gaussian_noise.py b/GPy/likelihoods/noise_models/gaussian_noise.py new file mode 100644 index 00000000..a77becb2 --- /dev/null +++ b/GPy/likelihoods/noise_models/gaussian_noise.py @@ -0,0 +1,89 @@ +# 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 Gaussian(NoiseDistribution): + """ + Gaussian likelihood + + :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 + self.variance = variance + super(Gaussian, self).__init__(gp_link,analytical_moments) + + def _preprocess_values(self,Y): + """ + Check if the values of the observations correspond to the values + assumed by the likelihood function. + """ + return Y + + def _moments_match_analytical(self,data_i,tau_i,v_i): + """ + Moments match of the marginal approximation in EP algorithm + + :param i: number of observation (int) + :param tau_i: precision of the cavity distribution (float) + :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 + 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) + + def _mass(self,gp,obs): + p = (self.gp_link.transf(gp)-self.mean)/np.sqrt(self.variance) + return std_norm_pdf(p) + + 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 + + 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 + + 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 + + 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 + """ + p = self.gp_link.transf(gp) + return p*(1-p) + + def _dvariance_dgp(self,gp): + return self.gp_link.dtransf_df(gp)*(1. - 2.*self.gp_link.transf(gp)) + + 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 diff --git a/GPy/likelihoods/noise_models/link_functions.py b/GPy/likelihoods/noise_models/gp_transformations.py similarity index 59% rename from GPy/likelihoods/noise_models/link_functions.py rename to GPy/likelihoods/noise_models/gp_transformations.py index b0cdcd49..b81e88e1 100644 --- a/GPy/likelihoods/noise_models/link_functions.py +++ b/GPy/likelihoods/noise_models/gp_transformations.py @@ -8,7 +8,7 @@ import scipy as sp import pylab as pb from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf,inv_std_norm_cdf -class LinkFunction(object): +class GPTransformation(object): """ Link function class for doing non-Gaussian likelihoods approximation @@ -18,82 +18,82 @@ class LinkFunction(object): def __init__(self): pass -class Identity(LinkFunction): +class Identity(GPTransformation): """ $$ g(f) = f $$ """ - def transf(self,mu): - return mu + #def transf(self,mu): + # return mu - def inv_transf(self,f): + def transf(self,f): return f - def dinv_transf_df(self,f): + def dtransf_df(self,f): return 1. - def d2inv_transf_df2(self,f): + def d2transf_df2(self,f): return 0 -class Probit(LinkFunction): +class Probit(GPTransformation): """ $$ g(f) = \\Phi^{-1} (mu) $$ """ - def transf(self,mu): - return inv_std_norm_cdf(mu) + #def transf(self,mu): + # return inv_std_norm_cdf(mu) - def inv_transf(self,f): + def transf(self,f): return std_norm_cdf(f) - def dinv_transf_df(self,f): + def dtransf_df(self,f): return std_norm_pdf(f) - def d2inv_transf_df2(self,f): + def d2transf_df2(self,f): return -f * std_norm_pdf(f) -class Log(LinkFunction): +class Log(GPTransformation): """ $$ g(f) = \log(\mu) $$ """ - def transf(self,mu): - return np.log(mu) + #def transf(self,mu): + # return np.log(mu) - def inv_transf(self,f): + def transf(self,f): return np.exp(f) - def dinv_transf_df(self,f): + def dtransf_df(self,f): return np.exp(f) - def d2inv_transf_df2(self,f): + def d2transf_df2(self,f): return np.exp(f) -class Log_ex_1(LinkFunction): +class Log_ex_1(GPTransformation): """ $$ g(f) = \log(\exp(\mu) - 1) $$ """ - def transf(self,mu): - """ - function: output space -> latent space - """ - return np.log(np.exp(mu) - 1) + #def transf(self,mu): + # """ + # function: output space -> latent space + # """ + # return np.log(np.exp(mu) - 1) - def inv_transf(self,f): + def transf(self,f): """ function: latent space -> output space """ return np.log(1.+np.exp(f)) - def dinv_transf_df(self,f): + def dtransf_df(self,f): return np.exp(f)/(1.+np.exp(f)) - def d2inv_transf_df2(self,f): + def d2transf_df2(self,f): aux = np.exp(f)/(1.+np.exp(f)) return aux*(1.-aux) diff --git a/GPy/likelihoods/noise_models/likelihood_functions.py b/GPy/likelihoods/noise_models/noise_distributions.py similarity index 97% rename from GPy/likelihoods/noise_models/likelihood_functions.py rename to GPy/likelihoods/noise_models/noise_distributions.py index dd597521..0dc9e03c 100644 --- a/GPy/likelihoods/noise_models/likelihood_functions.py +++ b/GPy/likelihoods/noise_models/noise_distributions.py @@ -8,19 +8,19 @@ import scipy as sp import pylab as pb from GPy.util.plot import gpplot from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf -import link_functions +import gp_transformations -class NoiseModel(object): +class NoiseDistribution(object): """ Likelihood class for doing Expectation propagation :param Y: observed output (Nx1 numpy.darray) ..Note:: Y values allowed depend on the LikelihoodFunction used """ - def __init__(self,link,analytical_moments=False): - #assert isinstance(link,link_functions.LinkFunction), "link is not a valid LinkFunction."#FIXME - self.link = link + def __init__(self,gp_link,analytical_moments=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.moments_match = self._moments_match_analytical @@ -289,7 +289,7 @@ class NoiseModel(object): :predictive_mean: output's predictive mean, if None _predictive_mean function will be called. """ qf = stats.norm.ppf(p,mu,sigma) - return self.link.inv_transf(qf) + return self.gp_link.transf(qf) def _nlog_joint_predictive_scaled(self,x,mu,sigma): """ @@ -334,7 +334,7 @@ class NoiseModel(object): :param mu: latent variable's predictive mean :param sigma: latent variable's predictive standard deviation """ - return sp.optimize.fmin_ncg(self._nlog_joint_predictive_scaled,x0=(mu,self.link.inv_transf(mu)),fprime=self._gradient_nlog_joint_predictive,fhess=self._hessian_nlog_joint_predictive,args=(mu,sigma)) + return sp.optimize.fmin_ncg(self._nlog_joint_predictive_scaled,x0=(mu,self.gp_link.transf(mu)),fprime=self._gradient_nlog_joint_predictive,fhess=self._hessian_nlog_joint_predictive,args=(mu,sigma)) def predictive_values(self,mu,var,sample=True,sample_size=5000): """ diff --git a/GPy/likelihoods/noise_models/poisson_likelihood.py b/GPy/likelihoods/noise_models/poisson_noise.py similarity index 58% rename from GPy/likelihoods/noise_models/poisson_likelihood.py rename to GPy/likelihoods/noise_models/poisson_noise.py index cdae29e3..e90b3ce8 100644 --- a/GPy/likelihoods/noise_models/poisson_likelihood.py +++ b/GPy/likelihoods/noise_models/poisson_noise.py @@ -5,12 +5,11 @@ import numpy as np from scipy import stats,special import scipy as sp -#import pylab as pb from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf -import link_functions -from likelihood_functions import NoiseModel +import gp_transformations +from noise_distributions import NoiseDistribution -class Poisson(NoiseModel): +class Poisson(NoiseDistribution): """ Poisson likelihood Y is expected to take values in {0,1,2,...} @@ -19,12 +18,12 @@ class Poisson(NoiseModel): L(x) = \exp(\lambda) * \lambda**Y_i / Y_i! $$ """ - def __init__(self,link=None,analytical_moments=False): + def __init__(self,gp_link=None,analytical_moments=False): #self.discrete = True #self.support_limits = (0,np.inf) #self.analytical_moments = False - super(Poisson, self).__init__(link,analytical_moments) + super(Poisson, self).__init__(gp_link,analytical_moments) def _preprocess_values(self,Y): #TODO #self.scale = .5*Y.max() @@ -36,57 +35,57 @@ class Poisson(NoiseModel): Mass (or density) function """ #obs = obs*self.scale + self.shift - return stats.poisson.pmf(obs,self.link.inv_transf(gp)) + return stats.poisson.pmf(obs,self.gp_link.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)) + 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): - return self.link.dinv_transf_df(gp) * (1. - obs/self.link.inv_transf(gp)) + return self.gp_link.dtransf_df(gp) * (1. - obs/self.gp_link.transf(gp)) def _d2nlog_mass_dgp2(self,gp,obs): - d2_df = self.link.d2inv_transf_df2(gp) - inv_transf = self.link.inv_transf(gp) - return obs * ((self.link.dinv_transf_df(gp)/inv_transf)**2 - d2_df/inv_transf) + d2_df + 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 _dnlog_mass_dobs(self,obs,gp): #TODO not needed - return special.psi(obs+1) - np.log(self.link.inv_transf(gp)) + return special.psi(obs+1) - np.log(self.gp_link.transf(gp)) def _d2nlog_mass_dobs2(self,obs,gp=None): #TODO not needed return special.polygamma(1,obs) def _d2nlog_mass_dcross(self,obs,gp): #TODO not needed - return -self.link.dinv_transf_df(gp)/self.link.inv_transf(gp) + return -self.gp_link.dtransf_df(gp)/self.gp_link.transf(gp) def _mean(self,gp): """ Mass (or density) function """ - return self.link.inv_transf(gp) + return self.gp_link.transf(gp) #def _variance(self,gp): - # return self.link.inv_transf(gp) + # return self.gp_link.transf(gp) def _dmean_dgp(self,gp): - return self.link.dinv_transf_df(gp) + return self.gp_link.dtransf_df(gp) def _d2mean_dgp2(self,gp): - return self.link.d2inv_transf_df2(gp) + return self.gp_link.d2transf_df2(gp) def _variance(self,gp): """ Mass (or density) function """ - return self.link.inv_transf(gp) + return self.gp_link.transf(gp) #def _variance(self,gp): - # return self.link.inv_transf(gp) + # return self.gp_link.transf(gp) def _dvariance_dgp(self,gp): - return self.link.dinv_transf_df(gp) + return self.gp_link.dtransf_df(gp) def _d2variance_dgp2(self,gp): - return self.link.d2inv_transf_df2(gp) + return self.gp_link.d2transf_df2(gp)