Initial commit, setting up the laplace approximation for a student t

python/likelihoods/Laplace.py

@@ -0,0 +1,54 @@
+import nump as np
+import GPy
+from GPy.util.linalg import jitchol
+
+class Laplace(GPy.likelihoods.likelihood):
+    """Laplace approximation to a posterior"""
+
+    def __init__(self,data,likelihood_function):
+        """
+        Laplace Approximation
+
+        First find the moments \hat{f} and the hessian at this point (using Newton-Raphson)
+        then find the z^{prime} which allows this to be a normalised gaussian instead of a
+        non-normalized gaussian
+
+        Finally we must compute the GP variables (i.e. generate some Y^{squiggle} and z^{squiggle}
+        which makes a gaussian the same as the laplace approximation
+
+        Arguments
+        ---------
+
+        :data: @todo
+        :likelihood_function: @todo
+
+        """
+        GPy.likelihoods.likelihood.__init__(self)
+
+        self.data = data
+        self.likelihood_function = likelihood_function
+
+        #Inital values
+        self.N, self.D = self.data.shape
+
+    def _compute_GP_variables(self):
+        """
+        Generates data Y which would give the normal distribution identical to the laplace approximation
+
+        GPy expects a likelihood to be gaussian, so need to caluclate the points Y^{squiggle} and Z^{squiggle}
+        that makes the posterior match that found by a laplace approximation to a non-gaussian likelihood
+        """
+        raise NotImplementedError
+
+    def fit_full(self, K):
+        """
+        The laplace approximation algorithm
+        For nomenclature see Rasmussen & Williams 2006
+        :K: Covariance matrix
+        """
+        self.f = np.zeros(self.N)
+
+        #Find \hat(f) using a newton raphson optimizer for example
+
+        #At this point get the hessian matrix
+

python/likelihoods/likelihood_function.py

@@ -0,0 +1,51 @@
+import GPy
+from scipy.special import gamma, gammaln
+
+class student_t(GPy.likelihoods.likelihood_function):
+    """Student t likelihood distribution
+    For nomanclature see Bayesian Data Analysis 2003 p576
+
+    Laplace:
+    Needs functions to calculate
+    ln p(yi|fi)
+    dln p(yi|fi)_dfi
+    d2ln p(yi|fi)_d2fi
+    """
+    def __init__(self, deg_free, sigma=1):
+        self.v = deg_free
+        self.sigma = 1
+
+    def link_function(self, y_i, f_i):
+        """link_function $\ln p(y_i|f_i)$
+
+        :y_i: datum number i
+        :f_i: latent variable f_i
+        :returns: float(likelihood evaluated for this point)
+
+        """
+        e = y_i - f_i
+        return gammaln((v+1)*0.5) - gammaln(v*0.5) - np.ln(v*np.pi*sigma)*0.5 - (v+1)*0.5*np.ln(1 + ((e/sigma)**2)/v)
+
+    def link_grad(self, y_i, f_i):
+        """gradient of the link function at y_i, given f_i w.r.t f_i
+
+        :y_i: datum number i
+        :f_i: latent variable f_i
+        :returns: float(gradient of likelihood evaluated at this point)
+
+        """
+        pass
+
+    def link_hess(self, y_i, f_i, f_j):
+        """hessian at this point (the hessian will be 0 unless i == j)
+        i.e. second derivative w.r.t f_i and f_j
+
+        :y_i: @todo
+        :f_i: @todo
+        :f_j: @todo
+        :returns: @todo
+
+        """
+        if f_i =
+        pass
+