Added log predictive density, student t degrees of freedom gradients and plotting functionality

GPy/core/gp.py

@@ -485,3 +485,38 @@ class GP(Model):
         """
         from ..inference.latent_function_inference.inferenceX import infer_newX
         return infer_newX(self, Y_new, optimize=optimize)
+
+    def log_predictive_density(self, x_test, y_test, Y_metadata=None):
+        """
+        Calculation of the log predictive density
+
+        .. math:
+            p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})
+
+        :param x_test: test locations (x_{*})
+        :type x_test: (Nx1) array
+        :param y_test: test observations (y_{*})
+        :type y_test: (Nx1) array
+        :param Y_metadata: metadata associated with the test points
+        """
+        mu_star, var_star = self._raw_predict(x_test)
+        return self.likelihood.log_predictive_density(y_test, mu_star, var_star, Y_metadata=Y_metadata)
+
+    def log_predictive_density_sampling(self, x_test, y_test, Y_metadata=None, num_samples=1000):
+        """
+        Calculation of the log predictive density by sampling
+
+        .. math:
+            p(y_{*}|D) = p(y_{*}|f_{*})p(f_{*}|\mu_{*}\\sigma^{2}_{*})
+
+        :param x_test: test locations (x_{*})
+        :type x_test: (Nx1) array
+        :param y_test: test observations (y_{*})
+        :type y_test: (Nx1) array
+        :param Y_metadata: metadata associated with the test points
+        :param num_samples: number of samples to use in monte carlo integration
+        :type num_samples: int
+        """
+        mu_star, var_star = self._raw_predict(x_test)
+        return self.likelihood.log_predictive_density_sampling(y_test, mu_star, var_star, Y_metadata=Y_metadata, num_samples=num_samples)
+

GPy/likelihoods/likelihood.py

@@ -41,6 +41,14 @@ class Likelihood(Parameterized):
         self.log_concave = False
         self.not_block_really = False
 
+    def request_num_latent_functions(self, Y):
+        """
+        The likelihood should infer how many latent functions are needed for the likelihood
+
+        Default is the number of outputs
+        """
+        return Y.shape[1]
+
     def _gradients(self,partial):
         return np.zeros(0)
 
@@ -118,15 +126,19 @@ class Likelihood(Parameterized):
             """Generate a function which can be integrated
             to give p(Y*|Y) = int p(Y*|f*)p(f*|Y) df*"""
             def f(fi_star):
-                #exponent = np.exp(-(1./(2*v))*np.square(m-f_star))
+                #exponent = np.exp(-(1./(2*vi))*np.square(mi-fi_star))
                 #from GPy.util.misc import safe_exp
                 #exponent = safe_exp(exponent)
-                #return self.pdf(f_star, y, y_m)*exponent
+                #res = safe_exp(self.logpdf(fi_star, yi, yi_m))*exponent
 
                 #More stable in the log space
-                return np.exp(self.logpdf(fi_star, yi, yi_m)
+                res = np.exp(self.logpdf(fi_star, yi, yi_m)
                               - 0.5*np.log(2*np.pi*vi)
-                              - 0.5*np.square(mi-fi_star)/vi)
+                              - 0.5*np.square(fi_star-mi)/vi)
+                if not np.isfinite(res):
+                    import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
+                return res
+
             return f
 
         p_ystar, _ = zip(*[quad(integral_generator(yi, mi, vi, yi_m), -np.inf, np.inf)
@@ -134,6 +146,36 @@ class Likelihood(Parameterized):
         p_ystar = np.array(p_ystar).reshape(-1, 1)
         return np.log(p_ystar)
 
+    def log_predictive_density_sampling(self, y_test, mu_star, var_star, Y_metadata=None, num_samples=1000):
+        """
+        Calculation of the log predictive density via sampling
+
+        .. math:
+            log p(y_{*}|D) = log 1/num_samples prod^{S}_{s=1} p(y_{*}|f_{*s})
+            f_{*s} ~ p(f_{*}|\mu_{*}\\sigma^{2}_{*})
+
+        :param y_test: test observations (y_{*})
+        :type y_test: (Nx1) array
+        :param mu_star: predictive mean of gaussian p(f_{*}|mu_{*}, var_{*})
+        :type mu_star: (Nx1) array
+        :param var_star: predictive variance of gaussian p(f_{*}|mu_{*}, var_{*})
+        :type var_star: (Nx1) array
+        :param num_samples: num samples of p(f_{*}|mu_{*}, var_{*}) to take
+        :type num_samples: int
+        """
+        assert y_test.shape==mu_star.shape
+        assert y_test.shape==var_star.shape
+        assert y_test.shape[1] == 1
+
+        #Take samples of p(f*|y)
+        #fi_samples = np.random.randn(num_samples)*np.sqrt(var_star) + mu_star
+        fi_samples = np.random.normal(mu_star, np.sqrt(var_star), size=(mu_star.shape[0], num_samples))
+
+        from scipy.misc import logsumexp
+        log_p_ystar = -np.log(num_samples) + logsumexp(self.logpdf(fi_samples, y_test, Y_metadata=Y_metadata), axis=1)
+        return log_p_ystar
+
+
     def _moments_match_ep(self,obs,tau,v):
         """
         Calculation of moments using quadrature

GPy/likelihoods/student_t.py

@@ -10,6 +10,7 @@ from scipy.special import gammaln, gamma
 from .likelihood import Likelihood
 from ..core.parameterization import Param
 from ..core.parameterization.transformations import Logexp
+from scipy.special import psi as digamma
 
 class StudentT(Likelihood):
     """
@@ -28,10 +29,10 @@ class StudentT(Likelihood):
         super(StudentT, self).__init__(gp_link, name='Student_T')
         # sigma2 is not a noise parameter, it is a squared scale.
         self.sigma2 = Param('t_scale2', float(sigma2), Logexp())
-        self.v = Param('deg_free', float(deg_free))
+        self.v = Param('deg_free', float(deg_free), Logexp())
         self.link_parameter(self.sigma2)
         self.link_parameter(self.v)
-        self.v.constrain_fixed()
+        #self.v.constrain_fixed()
 
         self.log_concave = False
 
@@ -224,20 +225,47 @@ class StudentT(Likelihood):
                            )
         return d2logpdf_dlink2_dvar
 
+    def dlogpdf_link_dv(self, inv_link_f, y, Y_metadata=None):
+        e = y - inv_link_f
+        e2 = np.square(e)
+        df = float(self.v[:])
+        s2 = float(self.sigma2[:])
+        dlogpdf_dv =  0.5*digamma(0.5*(df+1)) - 0.5*digamma(0.5*df) - 1.0/(2*df)
+        dlogpdf_dv += (1.0/(2*df))*(df+1)*e/(e2 + s2*df)
+        dlogpdf_dv -= np.log(1 + e2/(s2*df))
+        return dlogpdf_dv
+
+    def dlogpdf_dlink_dv(self, inv_link_f, y, Y_metadata=None):
+        e = y - inv_link_f
+        e2 = np.square(e)
+        df = float(self.v[:])
+        s2 = float(self.sigma2[:])
+        dlogpdf_df_dv = e*(e2 - self.sigma2)/(e2 + s2*df)**2
+        return dlogpdf_df_dv
+
+    def d2logpdf_dlink2_dv(self, inv_link_f, y, Y_metadata=None):
+        e = y - inv_link_f
+        e2 = np.square(e)
+        df = float(self.v[:])
+        s2 = float(self.sigma2[:])
+        #derivative of hess = ((self.v + 1)*(e**2 - self.v*self.sigma2)) / ((self.sigma2*self.v + e**2)**2)
+        e2_s2v = e**2 + s2*df
+        d2logpdf_df2_dv = (e2 - s2*df - s2*(df + 1))/e2_s2v**2 - 2*s2*(df+1)*(e2 - s2*df)/e2_s2v
+        return d2logpdf_df2_dv
+
     def dlogpdf_link_dtheta(self, f, y, Y_metadata=None):
         dlogpdf_dvar = self.dlogpdf_link_dvar(f, y, Y_metadata=Y_metadata)
-        dlogpdf_dv = np.zeros_like(dlogpdf_dvar) #FIXME: Not done yet
+        dlogpdf_dv = self.dlogpdf_link_dv(f, y, Y_metadata=Y_metadata)
         return np.array((dlogpdf_dvar, dlogpdf_dv))
 
     def dlogpdf_dlink_dtheta(self, f, y, Y_metadata=None):
         dlogpdf_dlink_dvar = self.dlogpdf_dlink_dvar(f, y, Y_metadata=Y_metadata)
-        dlogpdf_dlink_dv = np.zeros_like(dlogpdf_dlink_dvar) #FIXME: Not done yet
+        dlogpdf_dlink_dv = self.dlogpdf_dlink_dv(f, y, Y_metadata=Y_metadata)
         return np.array((dlogpdf_dlink_dvar, dlogpdf_dlink_dv))
 
     def d2logpdf_dlink2_dtheta(self, f, y, Y_metadata=None):
         d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, Y_metadata=Y_metadata)
-        d2logpdf_dlink2_dv = np.zeros_like(d2logpdf_dlink2_dvar) #FIXME: Not done yet
+        d2logpdf_dlink2_dv = self.d2logpdf_dlink2_dv(f, y, Y_metadata=Y_metadata)
-
         return np.array((d2logpdf_dlink2_dvar, d2logpdf_dlink2_dv))
 
     def predictive_mean(self, mu, sigma, Y_metadata=None):

GPy/plotting/matplot_dep/models_plots.py

@@ -216,7 +216,7 @@ def plot_fit_f(model, *args, **kwargs):
     kwargs['plot_raw'] = True
     plot_fit(model,*args, **kwargs)
 
-def fixed_inputs(model, non_fixed_inputs, fix_routine='median'):
+def fixed_inputs(model, non_fixed_inputs, fix_routine='median', as_list=True):
     """
     Convenience function for returning back fixed_inputs where the other inputs
     are fixed using fix_routine
@@ -226,6 +226,8 @@ def fixed_inputs(model, non_fixed_inputs, fix_routine='median'):
     :type non_fixed_inputs: list
     :param fix_routine: fixing routine to use, 'mean', 'median', 'zero'
     :type fix_routine: string
+    :param as_list: if true, will return a list of tuples with (dimension, fixed_val) otherwise it will create the corresponding X matrix
+    :type as_list: boolean
     """
     f_inputs = []
     if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs():
@@ -238,6 +240,11 @@ def fixed_inputs(model, non_fixed_inputs, fix_routine='median'):
                 f_inputs.append( (i, np.mean(X[:,i])) )
             if fix_routine == 'median':
                 f_inputs.append( (i, np.median(X[:,i])) )
-            elif fix_routine == 'zero':
+            else: # set to zero zero
                 f_inputs.append( (i, 0) )
-    return f_inputs
+            if not as_list:
+                X[:,i] = f_inputs[-1][1]
+    if as_list:
+        return f_inputs
+    else:
+        return X