diff --git a/GPy/likelihoods/binomial.py b/GPy/likelihoods/binomial.py
index 22009968..32fab8a5 100644
--- a/GPy/likelihoods/binomial.py
+++ b/GPy/likelihoods/binomial.py
@@ -66,14 +66,17 @@ class Binomial(Likelihood):
         :rtype: float
         """
         N = Y_metadata['trials']
+        assert N.shape == y.shape
         nchoosey = special.gammaln(N+1) - special.gammaln(y+1) - special.gammaln(N-y+1)
-
         return nchoosey + y*np.log(inv_link_f) + (N-y)*np.log(1.-inv_link_f)
 
     def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None):
         """
         Gradient of the pdf at y, given inverse link of f w.r.t inverse link of f.
 
+        .. math::
+            \\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{y_{i}}{\\lambda(f)} - \\frac{(N-y_{i})}{(1-\\lambda(f))}
+
         :param inv_link_f: latent variables inverse link of f.
         :type inv_link_f: Nx1 array
         :param y: data
@@ -83,7 +86,8 @@ class Binomial(Likelihood):
         :rtype: Nx1 array
         """
         N = Y_metadata['trials']
-        return y/inv_link_f - (N-y)/(1-inv_link_f)
+        assert N.shape == y.shape
+        return y/inv_link_f - (N-y)/(1.-inv_link_f)
 
     def d2logpdf_dlink2(self, inv_link_f, y, Y_metadata=None):
         """
@@ -92,7 +96,7 @@ class Binomial(Likelihood):
 
 
         .. math::
-            \\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{-y_{i}}{\\lambda(f)^{2}} - \\frac{(1-y_{i})}{(1-\\lambda(f))^{2}}
+            \\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{-y_{i}}{\\lambda(f)^{2}} - \\frac{(N-y_{i})}{(1-\\lambda(f))^{2}}
 
         :param inv_link_f: latent variables inverse link of f.
         :type inv_link_f: Nx1 array
@@ -107,7 +111,32 @@ class Binomial(Likelihood):
             (the distribution for y_i depends only on inverse link of f_i not on inverse link of f_(j!=i)
         """
         N = Y_metadata['trials']
-        return -y/np.square(inv_link_f) - (N-y)/np.square(1-inv_link_f)
+        assert N.shape == y.shape
+        return -y/np.square(inv_link_f) - (N-y)/np.square(1.-inv_link_f)
+
+    def d3logpdf_dlink3(self, inv_link_f, y, Y_metadata=None):
+        """
+        Third order derivative log-likelihood function at y given inverse link of f w.r.t inverse link of f
+
+        .. math::
+            \\frac{d^{2}\\ln p(y_{i}|\\lambda(f_{i}))}{d\\lambda(f)^{2}} = \\frac{2y_{i}}{\\lambda(f)^{3}} - \\frac{2(N-y_{i})}{(1-\\lambda(f))^{3}}
+
+        :param inv_link_f: latent variables inverse link of f.
+        :type inv_link_f: Nx1 array
+        :param y: data
+        :type y: Nx1 array
+        :param Y_metadata: Y_metadata not used in binomial
+        :returns: Diagonal of log hessian matrix (second derivative of log likelihood evaluated at points inverse link of f.
+        :rtype: Nx1 array
+
+        .. Note::
+            Will return diagonal of hessian, since every where else it is 0, as the likelihood factorizes over cases
+            (the distribution for y_i depends only on inverse link of f_i not on inverse link of f_(j!=i)
+        """
+        N = Y_metadata['trials']
+        assert N.shape == y.shape
+        inv_link_f2 = np.square(inv_link_f)
+        return 2*y/inv_link_f**3 - 2*(N-y)/(1.-inv_link_f)**3
 
     def samples(self, gp, Y_metadata=None):
         """
diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py
index 2fb255c9..346b4b83 100644
--- a/GPy/testing/likelihood_tests.py
+++ b/GPy/testing/likelihood_tests.py
@@ -117,7 +117,10 @@ class TestNoiseModels(object):
         self.positive_Y = np.exp(self.Y.copy())
         tmp = np.round(self.X[:, 0]*3-3)[:, None] + np.random.randint(0,3, self.X.shape[0])[:, None]
         self.integer_Y = np.where(tmp > 0, tmp, 0)
-
+        self.ns = np.random.poisson(50, size=self.N)[:, None]
+        p = np.abs(np.cos(2*np.pi*self.X + np.random.normal(scale=.2, size=(self.N, self.D)))).mean(1)
+        self.binomial_Y = np.array([np.random.binomial(self.ns[i], p[i]) for i in range(p.shape[0])])[:, None]
+        
         self.var = 0.2
         self.deg_free = 4.0
 
@@ -204,15 +207,6 @@ class TestNoiseModels(object):
                 },
                 "laplace": True
             },
-            #"Student_t_log": {
-            #"model": GPy.likelihoods.StudentT(gp_link=link_functions.Log(), deg_free=5, sigma2=self.var),
-            #"grad_params": {
-            #"names": [".*t_noise"],
-            #"vals": [self.var],
-            #"constraints": [(".*t_noise", self.constrain_positive), (".*deg_free", self.constrain_fixed)]
-            #},
-            #"laplace": True
-            #},
             "Gaussian_default": {
                 "model": GPy.likelihoods.Gaussian(variance=self.var),
                 "grad_params": {
@@ -273,6 +267,13 @@ class TestNoiseModels(object):
                 "laplace": True,
                 "ep": False #Should work though...
             },
+            "Binomial_default": {
+                "model": GPy.likelihoods.Binomial(),
+                "link_f_constraints": [partial(self.constrain_bounded, lower=0, upper=1)],
+                "Y": self.binomial_Y,
+                "Y_metadata": {'trials': self.ns},
+                "laplace": True,
+            },
             #,
             #GAMMA needs some work!"Gamma_default": {
             #"model": GPy.likelihoods.Gamma(),