diff --git a/python/examples/laplace_approximations.py b/python/examples/laplace_approximations.py
new file mode 100644
index 00000000..2f059831
--- /dev/null
+++ b/python/examples/laplace_approximations.py
@@ -0,0 +1,37 @@
+import GPy
+import numpy as np
+import scipy as sp
+import scipy.stats
+import matplotlib.pyplot as plt
+
+
+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]
+    Y = np.sin(X)
+
+    #Add some extreme value noise to some of the datapoints
+    percent_corrupted = 0.05
+    corrupted_datums = int(np.round(Y.shape[0] * percent_corrupted))
+    indices = np.arange(Y.shape[0])
+    np.random.shuffle(indices)
+    corrupted_indices = indices[:corrupted_datums]
+    print corrupted_indices
+    noise = np.random.uniform(-10,10,(len(corrupted_indices), 1))
+    Y[corrupted_indices] += noise
+
+    #A GP should completely break down due to the points as they get a lot of weight
+    # create simple GP model
+    m = GPy.models.GP_regression(X,Y)
+
+    # optimize
+    m.ensure_default_constraints()
+    m.optimize()
+    # plot
+    m.plot()
+    print m
+
+    #with a student t distribution, since it has heavy tails it should work well
diff --git a/python/likelihoods/Laplace.py b/python/likelihoods/Laplace.py
new file mode 100644
index 00000000..a0dbc65c
--- /dev/null
+++ b/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
+
diff --git a/python/likelihoods/likelihood_function.py b/python/likelihoods/likelihood_function.py
new file mode 100644
index 00000000..fd19675b
--- /dev/null
+++ b/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
+
diff --git a/python/models/coxGP.py b/python/models/coxGP.py
new file mode 100644
index 00000000..f61a8f46
--- /dev/null
+++ b/python/models/coxGP.py
@@ -0,0 +1,19 @@
+# Copyright (c) 2013, Alan Saul
+
+from GPy.models import GP
+from .. import likelihoods
+from GPy import kern
+
+
+class cox_GP_regression(GP):
+    """
+    Cox Gaussian Process model for regression
+    """
+
+    def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False, Xslices=None):
+        if kernel is None:
+            kernel = kern.rbf(X.shape[1])
+
+        likelihood = likelihoods.cox_piecewise(Y, normalize=normalize_Y)
+
+        GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, Xslices=Xslices)
diff --git a/python/testing/cox_tests.py b/python/testing/cox_tests.py
new file mode 100644
index 00000000..526f5c92
--- /dev/null
+++ b/python/testing/cox_tests.py
@@ -0,0 +1,14 @@
+# Copyright (c) 2013, Alan Saul
+
+import unittest
+import numpy as np
+import GPy
+
+class coxGPTests(unittest.TestCase):
+    def test_laplace_approx(self):
+        pass
+
+if __name__ == "__main__":
+    print "Running unit tests, please be (very) patient..."
+    unittest.main()
+