diff --git a/python/examples/laplace_approximations.py b/python/examples/laplace_approximations.py
index 0e1d3305..5642d8a4 100644
--- a/python/examples/laplace_approximations.py
+++ b/python/examples/laplace_approximations.py
@@ -1,7 +1,7 @@
 import GPy
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.stats import t
+from scipy.stats import t, norm
 from coxGP.python.likelihoods.Laplace import Laplace
 from coxGP.python.likelihoods.likelihood_function import student_t
 
@@ -11,12 +11,13 @@ def student_t_approx():
     Example of regressing with a student t likelihood
     """
     #Start a function, any function
-    X = np.sort(np.random.uniform(0, 15, 70))[:, None]
+    X = np.sort(np.random.uniform(0, 15, 100))[:, None]
     Y = np.sin(X)
 
     #Add student t random noise to datapoints
-    deg_free = 1
-    noise = t.rvs(deg_free, loc=1.8, scale=1, size=Y.shape)
+    deg_free = 2.5
+    t_rv = t(deg_free, loc=5, scale=1)
+    noise = t_rv.rvs(size=Y.shape)
     Y += noise
 
     # Kernel object
@@ -39,6 +40,19 @@ def student_t_approx():
     lap = Laplace(Y, likelihood_function)
     cov = kernel.K(X)
     lap.fit_full(cov)
+    #Get one sample (just look at a single Y
+    mode = float(lap.f_hat[0])
+    variance = float((deg_free/(deg_free-2))) #BUG: Not convinced this is giving reasonable variables
+    #variance = float((deg_free/(deg_free-2)) + np.diagonal(lap.hess_hat)[0]) #BUG: Not convinced this is giving reasonable variables
+    normalised_approx = norm(loc=mode, scale=variance)
+    print "Normal with mode %f, and variance %f" % (mode, variance)
+    print lap.height_unnormalised
+
+    test_range = np.arange(0, 10, 0.1)
+    print np.diagonal(lap.hess_hat)
+    plt.plot(test_range, t_rv.pdf(test_range))
+    plt.plot(test_range, normalised_approx.pdf(test_range))
+    plt.show()
 
 
 def noisy_laplace_approx():
diff --git a/python/likelihoods/Laplace.py b/python/likelihoods/Laplace.py
index 6efbfa30..08ae0e6f 100644
--- a/python/likelihoods/Laplace.py
+++ b/python/likelihoods/Laplace.py
@@ -5,13 +5,13 @@ from GPy.util.linalg import jitchol
 from functools import partial
 from GPy.likelihoods.likelihood import likelihood
 from GPy.util.linalg import pdinv,mdot
-
+from scipy.stats import norm
 
 
 class Laplace(likelihood):
     """Laplace approximation to a posterior"""
 
-    def __init__(self,data,likelihood_function):
+    def __init__(self, data, likelihood_function):
         """
         Laplace Approximation
 
@@ -42,7 +42,13 @@ class Laplace(likelihood):
         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
         """
-        z_hat = N(f_hat|f_hat, hess_hat) / self.height_unnormalised
+        #z_hat = N(f_hat|f_hat, hess_hat) / self.height_unnormalised
+        normalised_approx = norm(loc=self.f_hat, scale=self.hess_hat)
+        self.Z = normalised_approx.pdf(self.f_hat)/self.height_unnormalised
+        #self.Y =
+        #self.YYT =
+        #self.covariance_matrix =
+        #self.precision =
 
     def fit_full(self, K):
         """
@@ -51,11 +57,9 @@ class Laplace(likelihood):
         :K: Covariance matrix
         """
         f = np.zeros((self.N, 1))
-        print K.shape
-        print f.shape
-        print self.data.shape
+        #K = np.diag(np.ones(self.N))
         (Ki, _, _, log_Kdet) = pdinv(K)
-        obj_constant = (0.5 * log_Kdet) - ((0.5 * self.N) * np.log(2*np.pi))
+        obj_constant = (0.5 * log_Kdet) - ((0.5 * self.N) * np.log(2 * np.pi))
 
         #Find \hat(f) using a newton raphson optimizer for example
         #TODO: Add newton-raphson as subclass of optimizer class
@@ -77,11 +81,12 @@ class Laplace(likelihood):
             return np.squeeze(res)
 
         self.f_hat = sp.optimize.fmin_ncg(obj, f, fprime=obj_grad, fhess=obj_hess)
+        print self.f_hat
 
         #At this point get the hessian matrix
-        self.hess_hat = obj_hess(f_hat)
+        self.hess_hat = obj_hess(self.f_hat)
 
         #Need to add the constant as we previously were trying to avoid computing it (seems like a small overhead though...)
-        self.height_unnormalised = obj(f_hat) #FIXME: Is it -1?
+        self.height_unnormalised = obj(self.f_hat) #FIXME: Is it -1?
 
-        return _compute_GP_variables()
+        return self._compute_GP_variables()
diff --git a/python/likelihoods/likelihood_function.py b/python/likelihoods/likelihood_function.py
index 78731199..46128de7 100644
--- a/python/likelihoods/likelihood_function.py
+++ b/python/likelihoods/likelihood_function.py
@@ -15,27 +15,27 @@ class student_t(likelihood_function):
     dln p(yi|fi)_dfi
     d2ln p(yi|fi)_d2fifj
     """
-    def __init__(self, deg_free, sigma=1):
+    def __init__(self, deg_free, sigma=2):
         self.v = deg_free
-        self.sigma = 1
+        self.sigma = sigma
 
     def link_function(self, y, f):
         """link_function $\ln p(y|f)$
         $$\ln p(y_{i}|f_{i}) = \ln \Gamma(\frac{v+1}{2}) - \ln \Gamma(\frac{v}{2})\sqrt{v \pi}\sigma - \frac{v+1}{2}\ln (1 + \frac{1}{v}\left(\frac{y_{i} - f_{i}}{\sigma}\right)^2$$
 
-        :y: datum number i
-        :f: latent variable f
+        :y: data
+        :f: latent variables f
         :returns: float(likelihood evaluated for this point)
 
         """
+        assert y.shape[0] == f.shape[0]
         e = y - f
-        #print "Link ", y.shape, f.shape, e.shape
         objective = (gammaln((self.v + 1) * 0.5)
-                - gammaln(self.v * 0.5)
-                + np.log(self.sigma * np.sqrt(self.v * np.pi))
-                - (self.v + 1) * 0.5
-                * np.log(1 + ((e**2 / self.sigma**2) / self.v))
-                )
+                     - gammaln(self.v * 0.5)
+                     + np.log(self.sigma * np.sqrt(self.v * np.pi))
+                     - (self.v + 1) * 0.5
+                     * np.log(1 + ((e**2 / self.sigma**2) / self.v))
+                     )
         return np.sum(objective)
 
     def link_grad(self, y, f):
@@ -44,13 +44,13 @@ class student_t(likelihood_function):
 
         $$\frac{d}{df}p(y_{i}|f_{i}) = \frac{(v + 1)(y - f)}{v \sigma^{2} + (y_{i} - f_{i})^{2}}$$
 
-        :y: datum number i
-        :f: latent variable f
-        :returns: float(gradient of likelihood evaluated at this point)
+        :y: data
+        :f: latent variables f
+        :returns: gradient of likelihood evaluated at points
 
         """
+        assert y.shape[0] == f.shape[0]
         e = y - f
-        #print "Grad ", y.shape, f.shape, e.shape
         grad = ((self.v + 1) * e) / (self.v * (self.sigma**2) + (e**2))
         return grad
 
@@ -63,10 +63,11 @@ class student_t(likelihood_function):
 
         $$\frac{d^{2}p(y_{i}|f_{i})}{df^{2}} = \frac{(v + 1)(y - f)}{v \sigma^{2} + (y_{i} - f_{i})^{2}}$$
 
-        :y: datum number i
-        :f: latent variable f
-        :returns: float(second derivative of likelihood evaluated at this point)
+        :y: data
+        :f: latent variables f
+        :returns: array which is diagonal of covariance matrix (second derivative of likelihood evaluated at points)
         """
+        assert y.shape[0] == f.shape[0]
         e = y - f
-        hess = ((self.v + 1) * e) / ((((self.sigma**2)*self.v) + e**2)**2)
+        hess = ((self.v + 1) * e) / ((((self.sigma**2) * self.v) + e**2)**2)
         return hess