diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index ce2b59de..faa77e80 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -128,7 +128,7 @@ class GP(GPBase):
             debug_this # @UndefinedVariable
         return mu, var
 
-    def predict(self, Xnew, which_parts='all', full_cov=False):
+    def predict(self, Xnew, which_parts='all', full_cov=False, likelihood_args=dict()):
         """
         Predict the function(s) at the new point(s) Xnew.
         Arguments
@@ -153,6 +153,6 @@ class GP(GPBase):
         mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts)
 
         # now push through likelihood
-        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov)
+        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, **likelihood_args)
 
         return mean, var, _025pm, _975pm
diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py
index 1ce5db1a..f8ef57df 100644
--- a/GPy/kern/parts/rbf.py
+++ b/GPy/kern/parts/rbf.py
@@ -165,9 +165,8 @@ class RBF(Kernpart):
 
     def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, target):
         self._psi_computations(Z, mu, S)
-        denom_deriv = S[:, None, :] / (self.lengthscale ** 3 + self.lengthscale * S[:, None, :])
-        d_length = self._psi1[:, :, None] * (self.lengthscale * np.square(self._psi1_dist / (self.lengthscale2 + S[:, None, :])) + denom_deriv)
         target[0] += np.sum(dL_dpsi1 * self._psi1 / self.variance)
+        d_length = self._psi1[:,:,None] * ((self._psi1_dist_sq - 1.)/(self.lengthscale*self._psi1_denom) +1./self.lengthscale)
         dpsi1_dlength = d_length * dL_dpsi1[:, :, None]
         if not self.ARD:
             target[1] += dpsi1_dlength.sum()
diff --git a/GPy/kern/parts/rbf_inv.py b/GPy/kern/parts/rbf_inv.py
index 05605056..3512f555 100644
--- a/GPy/kern/parts/rbf_inv.py
+++ b/GPy/kern/parts/rbf_inv.py
@@ -159,21 +159,20 @@ class RBFInv(RBF):
 
     def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, target):
         self._psi_computations(Z, mu, S)
-        # #d_length = self._psi1[:, :, None] * (-0.5 * ((np.square((self._psi1_dist)/(self.inv_lengthscale * S[:,None,:] + 1))) + ((S[:, None, :])/(self.inv_lengthscale * S[:, None, :] + 1))))
-        tmp = 1 + S[:, None, :] * self.inv_lengthscale2
-        # inv_len3 = np.power(self.inv_lengthscale,3)
-        d_length = -(self._psi1[:, :, None] * ((np.square(self._psi1_dist) * self.inv_lengthscale) / (tmp ** 2) + (S[:, None, :] * self.inv_lengthscale) / (tmp)))
+        tmp = 1 + S[:,None,:]*self.inv_lengthscale2
+        d_inv_length_old = -self._psi1[:,:,None] * ((self._psi1_dist_sq - 1.)/(self.lengthscale*self._psi1_denom) + self.inv_lengthscale)/self.inv_lengthscale2
+        d_inv_length = -self._psi1[:,:,None] * ((self._psi1_dist_sq - 1.)/self._psi1_denom + self.lengthscale)
         target[0] += np.sum(dL_dpsi1 * self._psi1 / self.variance)
-        dpsi1_dlength = d_length * dL_dpsi1[:, :, None]
+        dpsi1_dlength = d_inv_length * dL_dpsi1[:, :, None]
         if not self.ARD:
-            target[1] += dpsi1_dlength.sum() # *(-self.lengthscale2)
+            target[1] += dpsi1_dlength.sum()#*(-self.lengthscale2)
         else:
-            target[1:] += dpsi1_dlength.sum(0).sum(0) # *(-self.lengthscale2)
-        # target[1:] = target[1:]*(-self.lengthscale2)
+            target[1:] += dpsi1_dlength.sum(0).sum(0)#*(-self.lengthscale2)
+        #target[1:] = target[1:]*(-self.lengthscale2)
 
     def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, target):
         self._psi_computations(Z, mu, S)
-        dpsi1_dZ = -self._psi1[:, :, None] * ((self.inv_lengthscale2 * self._psi1_dist) / self._psi1_denom)
+        dpsi1_dZ = -self._psi1[:, :, None] * ((self.inv_lengthscale2*self._psi1_dist)/self._psi1_denom)
         target += np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0)
 
     def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S, target_mu, target_S):
@@ -186,15 +185,15 @@ class RBFInv(RBF):
         """Shape N,num_inducing,num_inducing,Ntheta"""
         self._psi_computations(Z, mu, S)
         d_var = 2.*self._psi2 / self.variance
-        # d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
+        #d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom)
         d_length = -2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] * self.inv_lengthscale2) / (self.inv_lengthscale * self._psi2_denom)
         target[0] += np.sum(dL_dpsi2 * d_var)
         dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None]
         if not self.ARD:
-            target[1] += dpsi2_dlength.sum() # *(-self.lengthscale2)
+            target[1] += dpsi2_dlength.sum()#*(-self.lengthscale2)
         else:
-            target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0) # *(-self.lengthscale2)
-        # target[1:] = target[1:]*(-self.lengthscale2)
+            target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)#*(-self.lengthscale2)
+        #target[1:] = target[1:]*(-self.lengthscale2)
 
     def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target):
         self._psi_computations(Z, mu, S)