diff --git a/GPy/likelihoods/noise_models/gaussian_noise.py b/GPy/likelihoods/noise_models/gaussian_noise.py index 8bce30b7..5811f916 100644 --- a/GPy/likelihoods/noise_models/gaussian_noise.py +++ b/GPy/likelihoods/noise_models/gaussian_noise.py @@ -68,7 +68,7 @@ class Gaussian(NoiseDistribution): def _predictive_variance_analytical(self,mu,sigma,predictive_mean=None): return 1./(1./self.variance + 1./sigma**2) - def _mass(self, link_f, y): + def pdf_link(self, link_f, y, extra_data=None): #FIXME: Careful now passing link_f in not gp (f)! #return std_norm_pdf( (self.gp_link.transf(gp)-obs)/np.sqrt(self.variance) ) #Assumes no covariance, exp, sum, log for numerical stability @@ -76,21 +76,26 @@ class Gaussian(NoiseDistribution): #return np.exp(np.sum(np.log(stats.norm.pdf(y, link_f, np.sqrt(self.variance))))) return np.exp(np.sum(np.log(stats.norm.pdf(y, link_f, np.sqrt(self.variance))))) + def _mass(self, link_f, y, extra_data=None): + NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\ + Please negate your function and use pdf in noise_model.py, if implementing a likelihood\ + rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\ + its derivatives") def _nlog_mass(self, link_f, y, extra_data=None): - NotImplementedError("Deprecated, now doing chain in likelihood.py for link function evaluation\ + NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\ Please negate your function and use logpdf in noise_model.py, if implementing a likelihood\ rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\ its derivatives") def _dnlog_mass_dgp(self, link_f, y, extra_data=None): - NotImplementedError("Deprecated, now doing chain in likelihood.py for link function evaluation\ - Please negate your function and use logpdf in noise_model.py, if implementing a likelihood\ + NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\ + Please negate your function and use dlogpdf_df in noise_model.py, if implementing a likelihood\ rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\ its derivatives") def _d2nlog_mass_dgp2(self, link_f, y, extra_data=None): - NotImplementedError("Deprecated, now doing chain in likelihood.py for link function evaluation\ - Please negate your function and use logpdf in noise_model.py, if implementing a likelihood\ + NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\ + Please negate your function and use d2logpdf_df2 in noise_model.py, if implementing a likelihood\ rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\ its derivatives") diff --git a/GPy/likelihoods/noise_models/gp_transformations.py b/GPy/likelihoods/noise_models/gp_transformations.py index c6e316e8..b9db75ce 100644 --- a/GPy/likelihoods/noise_models/gp_transformations.py +++ b/GPy/likelihoods/noise_models/gp_transformations.py @@ -80,6 +80,10 @@ class Probit(GPTransformation): def d2transf_df2(self,f): return -f * std_norm_pdf(f) + def d3transf_df3(self,f): + f2 = f**2 + return -(1/(np.sqrt(2*np.pi)))*np.exp(-0.5*(f2))*(f2-1) + class Log(GPTransformation): """ .. math:: @@ -96,6 +100,9 @@ class Log(GPTransformation): def d2transf_df2(self,f): return np.exp(f) + def d3transf_df3(self,f): + return np.exp(f) + class Log_ex_1(GPTransformation): """ .. math:: diff --git a/GPy/likelihoods/noise_models/noise_distributions.py b/GPy/likelihoods/noise_models/noise_distributions.py index 6b36f42b..0516a735 100644 --- a/GPy/likelihoods/noise_models/noise_distributions.py +++ b/GPy/likelihoods/noise_models/noise_distributions.py @@ -399,16 +399,82 @@ class NoiseDistribution(object): """ 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),disp=False) - def logpdf(self, f, y, extra_data=None): + def pdf_link(self, link_f, y, extra_data=None): + raise NotImplementedError + + def logpdf_link(self, link_f, y, extra_data=None): + raise NotImplementedError + + def dlogpdf_dlink(self, link_f, y, extra_data=None): + raise NotImplementedError + + def d2logpdf_dlink2(self, link_f, y, extra_data=None): + raise NotImplementedError + + def d3logpdf_dlink3(self, link_f, y, extra_data=None): + raise NotImplementedError + + def dlogpdf_link_dtheta(self, link_f, y, extra_data=None): + raise NotImplementedError + + def dlogpdf_dlink_dtheta(self, link_f, y, extra_data=None): + raise NotImplementedError + + def d2logpdf_dlink2_dtheta(self, link_f, y, extra_data=None): + raise NotImplementedError + + + def pdf(self, f, y, extra_data=None): """ - Evaluates the link function link(f) then computes the log likelihood using it + Evaluates the link function link(f) then computes the likelihood (pdf) using it + + .. math: + p(y|\\lambda(f)) + + :param f: latent variables f + :type f: Nx1 array + :param y: data + :type y: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: likelihood evaluated for this point + :rtype: float """ link_f = self.gp_link.transf(f) - return self.logpdf_link(f, y, extra_data=extra_data) + return self.pdf_link(link_f, y, extra_data=extra_data) + + def logpdf(self, f, y, extra_data=None): + """ + Evaluates the link function link(f) then computes the log likelihood (log pdf) using it + + .. math: + \\log p(y|\\lambda(f)) + + :param f: latent variables f + :type f: Nx1 array + :param y: data + :type y: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: log likelihood evaluated for this point + :rtype: float + """ + link_f = self.gp_link.transf(f) + return self.logpdf_link(link_f, y, extra_data=extra_data) def dlogpdf_df(self, f, y, extra_data=None): """ - TODO: Doc strings + Evaluates the link function link(f) then computes the derivative of log likelihood using it + Uses the Faa di Bruno's formula for the chain rule + + .. math:: + \\frac{d\\log p(y|\\lambda(f))}{df} = \\frac{d\\log p(y|\\lambda(f))}{d\\lambda(f)}\\frac{d\\lambda(f)}{df} + + :param f: latent variables f + :type f: Nx1 array + :param y: data + :type y: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: derivative of log likelihood evaluated for this point + :rtype: float """ link_f = self.gp_link.transf(f) dlogpdf_dlink = self.dlogpdf_dlink(link_f, y, extra_data=extra_data) @@ -417,7 +483,19 @@ class NoiseDistribution(object): def d2logpdf_df2(self, f, y, extra_data=None): """ - TODO: Doc strings + Evaluates the link function link(f) then computes the second derivative of log likelihood using it + Uses the Faa di Bruno's formula for the chain rule + + .. math:: + \\frac{d^{2}\\log p(y|\\lambda(f))}{df^{2}} = \\frac{d^{2}\\log p(y|\\lambda(f))}{d^{2}\\lambda(f)}\\left(\\frac{d\\lambda(f)}{df}\\right)^{2} + \\frac{d\\log p(y|\\lambda(f))}{d\\lambda(f)}\\frac{d^{2}\\lambda(f)}{df^{2}} + + :param f: latent variables f + :type f: Nx1 array + :param y: data + :type y: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: second derivative of log likelihood evaluated for this point + :rtype: float """ link_f = self.gp_link.transf(f) d2logpdf_dlink2 = self.d2logpdf_dlink2(link_f, y, extra_data=extra_data) @@ -428,7 +506,19 @@ class NoiseDistribution(object): def d3logpdf_df3(self, f, y, extra_data=None): """ - TODO: Doc strings + Evaluates the link function link(f) then computes the third derivative of log likelihood using it + Uses the Faa di Bruno's formula for the chain rule + + .. math:: + \\frac{d^{3}\\log p(y|\\lambda(f))}{df^{3}} = \\frac{d^{3}\\log p(y|\\lambda(f)}{d\\lambda(f)^{3}}\\left(\\frac{d\\lambda(f)}{df}\\right)^{3} + 3\\frac{d^{2}\\log p(y|\\lambda(f)}{d\\lambda(f)^{2}}\\frac{d\\lambda(f)}{df}\\frac{d^{2}\\lambda(f)}{df^{2}} + \\frac{d\\log p(y|\\lambda(f)}{d\\lambda(f)}\\frac{d^{3}\\lambda(f)}{df^{3}} + + :param f: latent variables f + :type f: Nx1 array + :param y: data + :type y: Nx1 array + :param extra_data: extra_data which is not used in student t distribution - not used + :returns: third derivative of log likelihood evaluated for this point + :rtype: float """ link_f = self.gp_link.transf(f) d3logpdf_dlink3 = self.d3logpdf_dlink3(link_f, y, extra_data=extra_data) @@ -440,23 +530,33 @@ class NoiseDistribution(object): return chain_3(d3logpdf_dlink3, dlink_df, d2logpdf_dlink2, d2link_df2, dlogpdf_dlink, d3link_df3) def dlogpdf_dtheta(self, f, y, extra_data=None): + """ + TODO: Doc strings + """ link_f = self.gp_link.transf(f) return self.dlogpdf_link_dtheta(link_f, y, extra_data=extra_data) def dlogpdf_df_dtheta(self, f, y, extra_data=None): + """ + TODO: Doc strings + """ link_f = self.gp_link.transf(f) dlink_df = self.gp_link.dtransf_df(f) dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data) return chain_1(dlogpdf_dlink_dtheta, dlink_df) def d2logpdf_df2_dtheta(self, f, y, extra_data=None): + """ + TODO: Doc strings + """ link_f = self.gp_link.transf(f) dlink_df = self.gp_link.dtransf_df(f) - d2link_df2 = self.gp_link.d2transf_df2(f) #FIXME: I THINK ITS THIS + d2link_df2 = self.gp_link.d2transf_df2(f) d2logpdf_dlink2_dtheta = self.d2logpdf_dlink2_dtheta(link_f, y, extra_data=extra_data) dlogpdf_dlink_dtheta = self.dlogpdf_dlink_dtheta(link_f, y, extra_data=extra_data) - return chain_2(d2logpdf_dlink2_dtheta, dlink_df, dlogpdf_dlink_dtheta, d2link_df2) + #FIXME: Why isn't this chain_1? #return chain_1(d2logpdf_dlink2_dtheta, d2link_df2) + return chain_2(d2logpdf_dlink2_dtheta, dlink_df, dlogpdf_dlink_dtheta, d2link_df2) def _laplace_gradients(self, f, y, extra_data=None): #link_f = self.gp_link.transf(f) @@ -508,14 +608,10 @@ class NoiseDistribution(object): q3 = np.vstack(q3) return pred_mean, pred_var, q1, q3 - def samples(self, gp): """ Returns a set of samples of observations based on a given value of the latent variable. :param gp: latent variable """ - pass - - - + raise NotImplementedError diff --git a/GPy/testing/laplace_tests.py b/GPy/testing/laplace_tests.py index dbdd34f3..1f20d9ae 100644 --- a/GPy/testing/laplace_tests.py +++ b/GPy/testing/laplace_tests.py @@ -4,6 +4,7 @@ import GPy from GPy.models import GradientChecker import functools import inspect +from GPy.likelihoods.noise_models import gp_transformations def dparam_partial(inst_func, *args): """ @@ -77,7 +78,7 @@ class LaplaceTests(unittest.TestCase): self.var = np.random.rand(1) self.stu_t = GPy.likelihoods.student_t(deg_free=5, sigma2=self.var) - self.gauss = GPy.likelihoods.gaussian(variance=self.var, D=self.D, N=self.N) + self.gauss = GPy.likelihoods.gaussian(gp_transformations.Log(), variance=self.var, D=self.D, N=self.N) #Make a bigger step as lower bound can be quite curved self.step = 1e-3 @@ -92,7 +93,7 @@ class LaplaceTests(unittest.TestCase): def test_mass_logpdf(self): print "\n{}".format(inspect.stack()[0][3]) np.testing.assert_almost_equal( - np.log(self.gauss._mass(self.f.copy(), self.Y.copy())), + np.log(self.gauss.pdf(self.f.copy(), self.Y.copy())), self.gauss.logpdf(self.f.copy(), self.Y.copy())) @@ -149,7 +150,7 @@ class LaplaceTests(unittest.TestCase): """ dGauss_dlink's """ def test_gaussian_dlogpdf_dlink(self): print "\n{}".format(inspect.stack()[0][3]) - logpdf = functools.partial(self.gauss.logpdf, y=self.Y) + logpdf = functools.partial(self.gauss.logpdf_link, y=self.Y) dlogpdf_dlink = functools.partial(self.gauss.dlogpdf_dlink, y=self.Y) grad = GradientChecker(logpdf, dlogpdf_dlink, self.f.copy(), 'g') grad.randomize() diff --git a/doc/GPy.testing.rst b/doc/GPy.testing.rst index ef25ba60..078a41a2 100644 --- a/doc/GPy.testing.rst +++ b/doc/GPy.testing.rst @@ -76,6 +76,14 @@ GPy.testing.mrd_tests module :undoc-members: :show-inheritance: +GPy.testing.noise_distributions module +-------------------------------------- + +.. automodule:: GPy.testing.noise_distributions + :members: + :undoc-members: + :show-inheritance: + GPy.testing.prior_tests module ------------------------------ diff --git a/doc/GPy.util.rst b/doc/GPy.util.rst index 5aca7cf9..f2aaed7f 100644 --- a/doc/GPy.util.rst +++ b/doc/GPy.util.rst @@ -27,6 +27,14 @@ GPy.util.classification module :undoc-members: :show-inheritance: +GPy.util.config module +---------------------- + +.. automodule:: GPy.util.config + :members: + :undoc-members: + :show-inheritance: + GPy.util.datasets module ------------------------ @@ -91,6 +99,14 @@ GPy.util.multioutput module :undoc-members: :show-inheritance: +GPy.util.netpbmfile module +-------------------------- + +.. automodule:: GPy.util.netpbmfile + :members: + :undoc-members: + :show-inheritance: + GPy.util.plot module --------------------