diff --git a/GPy/core/fitc.py b/GPy/core/fitc.py
index 670c0354..684d232e 100644
--- a/GPy/core/fitc.py
+++ b/GPy/core/fitc.py
@@ -20,7 +20,7 @@ class FITC(SparseGP):
     sparse FITC approximation
 
     :param X: inputs
-    :type X: np.ndarray (N x Q)
+    :type X: np.ndarray (num_data x Q)
     :param likelihood: a likelihood instance, containing the observed data
     :type likelihood: GPy.likelihood.(Gaussian | EP)
     :param kernel : the kernel (covariance function). See link kernels
@@ -112,7 +112,7 @@ class FITC(SparseGP):
         else:
             if self.likelihood.is_heteroscedastic:
                 assert self.likelihood.output_dim == 1
-            tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.N)))
+            tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.num_data)))
             tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
             self.A = tdot(tmp)
 
@@ -168,7 +168,7 @@ class FITC(SparseGP):
         self._dpsi1_dX_jkj = 0
         self._dpsi1_dtheta_jkj = 0
 
-        for i,V_n,alpha_n,gamma_n,gamma_k in zip(range(self.N),self.V_star,alpha,gamma_2,gamma_3):
+        for i,V_n,alpha_n,gamma_n,gamma_k in zip(range(self.num_data),self.V_star,alpha,gamma_2,gamma_3):
             K_pp_K = np.dot(Kmmipsi1[:,i:(i+1)],Kmmipsi1[:,i:(i+1)].T)
 
             #Diag_dpsi1 = Diag_dA_dpsi1: yT*beta_star*y + Diag_dC_dpsi1 +Diag_dD_dpsi1
@@ -215,7 +215,7 @@ class FITC(SparseGP):
 
     def log_likelihood(self):
         """ Compute the (lower bound on the) log marginal likelihood """
-        A = -0.5 * self.N * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.beta_star)) - 0.5 * np.sum(self.V_star * self.likelihood.Y)
+        A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.beta_star)) - 0.5 * np.sum(self.V_star * self.likelihood.Y)
         C = -self.output_dim * (np.sum(np.log(np.diag(self.LB))))
         D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
         return A + C + D
diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index f0dd4fed..246b8cc9 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -46,12 +46,12 @@ class GP(GPBase):
             #alpha = np.dot(self.Ki, self.likelihood.Y)
             alpha,_ = linalg.lapack.flapack.dpotrs(self.L, self.likelihood.Y,lower=1)
 
-            self.dL_dK = 0.5 * (tdot(alpha) - self.input_dim * self.Ki)
+            self.dL_dK = 0.5 * (tdot(alpha) - self.output_dim * self.Ki)
         else:
             #tmp = mdot(self.Ki, self.likelihood.YYT, self.Ki)
             tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(self.likelihood.YYT), lower=1)
             tmp, _ = linalg.lapack.flapack.dpotrs(self.L, np.asfortranarray(tmp.T), lower=1)
-            self.dL_dK = 0.5 * (tmp - self.input_dim * self.Ki)
+            self.dL_dK = 0.5 * (tmp - self.output_dim * self.Ki)
 
     def _get_params(self):
         return np.hstack((self.kern._get_params_transformed(), self.likelihood._get_params()))
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index c236115e..9188fe6f 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -13,12 +13,12 @@ class GPBase(Model):
     def __init__(self, X, likelihood, kernel, normalize_X=False):
         self.X = X
         assert len(self.X.shape) == 2
-        self.N, self.input_dim = self.X.shape
+        self.num_data, self.input_dim = self.X.shape
         assert isinstance(kernel, kern.kern)
         self.kern = kernel
         self.likelihood = likelihood
         assert self.X.shape[0] == self.likelihood.data.shape[0]
-        self.N, self.output_dim = self.likelihood.data.shape
+        self.num_data, self.output_dim = self.likelihood.data.shape
 
         if normalize_X:
             self._Xmean = X.mean(0)[None, :]
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index 8a10ca28..2cfc8ae4 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -13,13 +13,13 @@ class SparseGP(GPBase):
     Variational sparse GP model
 
     :param X: inputs
-    :type X: np.ndarray (N x input_dim)
+    :type X: np.ndarray (num_data x input_dim)
     :param likelihood: a likelihood instance, containing the observed data
     :type likelihood: GPy.likelihood.(Gaussian | EP | Laplace)
     :param kernel : the kernel (covariance function). See link kernels
     :type kernel: a GPy.kern.kern instance
     :param X_variance: The uncertainty in the measurements of X (Gaussian variance)
-    :type X_variance: np.ndarray (N x input_dim) | None
+    :type X_variance: np.ndarray (num_data x input_dim) | None
     :param Z: inducing inputs (optional, see note)
     :type Z: np.ndarray (num_inducing x input_dim) | None
     :param num_inducing : Number of inducing points (optional, default 10. Ignored if Z is not None)
@@ -69,7 +69,7 @@ class SparseGP(GPBase):
         # The rather complex computations of self.A
         if self.has_uncertain_inputs:
             if self.likelihood.is_heteroscedastic:
-                psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0)
+                psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.num_data, 1, 1))).sum(0)
             else:
                 psi2_beta = self.psi2.sum(0) * self.likelihood.precision
             evals, evecs = linalg.eigh(psi2_beta)
@@ -77,7 +77,7 @@ class SparseGP(GPBase):
             tmp = evecs * np.sqrt(clipped_evals)
         else:
             if self.likelihood.is_heteroscedastic:
-                tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)))
+                tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.num_data)))
             else:
                 tmp = self.psi1 * (np.sqrt(self.likelihood.precision))
         tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
@@ -99,28 +99,28 @@ class SparseGP(GPBase):
 
         # Compute dL_dKmm
         tmp = tdot(self._LBi_Lmi_psi1V)
-        self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.input_dim * np.eye(self.num_inducing) + tmp)
+        self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.output_dim * np.eye(self.num_inducing) + tmp)
         tmp = -0.5 * self.DBi_plus_BiPBi
-        tmp += -0.5 * self.B * self.input_dim
-        tmp += self.input_dim * np.eye(self.num_inducing)
+        tmp += -0.5 * self.B * self.output_dim
+        tmp += self.output_dim * np.eye(self.num_inducing)
         self.dL_dKmm = backsub_both_sides(self.Lm, tmp)
 
         # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertain inputs case
-        self.dL_dpsi0 = -0.5 * self.input_dim * (self.likelihood.precision * np.ones([self.N, 1])).flatten()
+        self.dL_dpsi0 = -0.5 * self.output_dim * (self.likelihood.precision * np.ones([self.num_data, 1])).flatten()
         self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T)
-        dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.input_dim * np.eye(self.num_inducing) - self.DBi_plus_BiPBi)
+        dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.output_dim * np.eye(self.num_inducing) - self.DBi_plus_BiPBi)
 
         if self.likelihood.is_heteroscedastic:
             if self.has_uncertain_inputs:
                 self.dL_dpsi2 = self.likelihood.precision.flatten()[:, None, None] * dL_dpsi2_beta[None, :, :]
             else:
-                self.dL_dpsi1 += 2.*np.dot(dL_dpsi2_beta, self.psi1 * self.likelihood.precision.reshape(1, self.N))
+                self.dL_dpsi1 += 2.*np.dot(dL_dpsi2_beta, self.psi1 * self.likelihood.precision.reshape(1, self.num_data))
                 self.dL_dpsi2 = None
         else:
             dL_dpsi2 = self.likelihood.precision * dL_dpsi2_beta
             if self.has_uncertain_inputs:
                 # repeat for each of the N psi_2 matrices
-                self.dL_dpsi2 = np.repeat(dL_dpsi2[None, :, :], self.N, axis=0)
+                self.dL_dpsi2 = np.repeat(dL_dpsi2[None, :, :], self.num_data, axis=0)
             else:
                 # subsume back into psi1 (==Kmn)
                 self.dL_dpsi1 += 2.*np.dot(dL_dpsi2, self.psi1)
@@ -135,17 +135,17 @@ class SparseGP(GPBase):
             raise NotImplementedError, "heteroscedatic derivates not implemented"
         else:
             # likelihood is not heterscedatic
-            self.partial_for_likelihood = -0.5 * self.N * self.input_dim * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
-            self.partial_for_likelihood += 0.5 * self.input_dim * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
+            self.partial_for_likelihood = -0.5 * self.num_data * self.output_dim * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
+            self.partial_for_likelihood += 0.5 * self.output_dim * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
             self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V)))
 
     def log_likelihood(self):
         """ Compute the (lower bound on the) log marginal likelihood """
         if self.likelihood.is_heteroscedastic:
-            A = -0.5 * self.N * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V * self.likelihood.Y)
+            A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V * self.likelihood.Y)
             B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
         else:
-            A = -0.5 * self.N * self.output_dim * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
+            A = -0.5 * self.num_data * self.output_dim * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
             B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
         C = -self.output_dim * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.num_inducing * np.log(sf2))
         D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py
index f4a44380..8043c635 100644
--- a/GPy/models/bayesian_gplvm.py
+++ b/GPy/models/bayesian_gplvm.py
@@ -73,8 +73,8 @@ class BayesianGPLVM(SparseGP, GPLVM):
         self._oldps.insert(0, p.copy())
 
     def _get_param_names(self):
-        X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
-        S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
+        X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
+        S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
         return (X_names + S_names + SparseGP._get_param_names(self))
 
     def _get_params(self):
@@ -96,7 +96,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
     def _set_params(self, x, save_old=True, save_count=0):
 #         try:
             x = self._clipped(x)
-            N, input_dim = self.N, self.input_dim
+            N, input_dim = self.num_data, self.input_dim
             self.X = x[:self.X.size].reshape(N, input_dim).copy()
             self.X_variance = x[(N * input_dim):(2 * N * input_dim)].reshape(N, input_dim).copy()
             SparseGP._set_params(self, x[(2 * N * input_dim):])
@@ -126,7 +126,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
     def KL_divergence(self):
         var_mean = np.square(self.X).sum()
         var_S = np.sum(self.X_variance - np.log(self.X_variance))
-        return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.N
+        return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.num_data
 
     def log_likelihood(self):
         ll = SparseGP.log_likelihood(self)
@@ -146,11 +146,11 @@ class BayesianGPLVM(SparseGP, GPLVM):
                 self._savedpsiKmm.append([self.f_call, [self.Kmm, self.dL_dKmm]])
 #                 sf2 = self.scale_factor ** 2
                 if self.likelihood.is_heteroscedastic:
-                    A = -0.5 * self.N * self.input_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.V * self.likelihood.Y)
+                    A = -0.5 * self.num_data * self.input_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.V * self.likelihood.Y)
 #                     B = -0.5 * self.input_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
                     B = -0.5 * self.input_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
                 else:
-                    A = -0.5 * self.N * self.input_dim * (np.log(2.*np.pi) + np.log(self.likelihood._variance)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
+                    A = -0.5 * self.num_data * self.input_dim * (np.log(2.*np.pi) + np.log(self.likelihood._variance)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
 #                     B = -0.5 * self.input_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2)
                     B = -0.5 * self.input_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
                 C = -self.input_dim * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.num_inducing * np.log(sf2))
@@ -266,9 +266,9 @@ class BayesianGPLVM(SparseGP, GPLVM):
 
     def _debug_filter_params(self, x):
         start, end = 0, self.X.size,
-        X = x[start:end].reshape(self.N, self.input_dim)
+        X = x[start:end].reshape(self.num_data, self.input_dim)
         start, end = end, end + self.X_variance.size
-        X_v = x[start:end].reshape(self.N, self.input_dim)
+        X_v = x[start:end].reshape(self.num_data, self.input_dim)
         start, end = end, end + (self.num_inducing * self.input_dim)
         Z = x[start:end].reshape(self.num_inducing, self.input_dim)
         start, end = end, end + self.input_dim
diff --git a/GPy/models/fitc.py b/GPy/models/fitc.py
index 785ea46c..5df1a7b5 100644
--- a/GPy/models/fitc.py
+++ b/GPy/models/fitc.py
@@ -52,7 +52,7 @@ class FITC(SparseGP):
         else:
             if self.likelihood.is_heteroscedastic:
                 assert self.likelihood.input_dim == 1
-            tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.N)))
+            tmp = self.psi1 * (np.sqrt(self.beta_star.flatten().reshape(1, self.num_data)))
             tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
             self.A = tdot(tmp)
 
@@ -108,7 +108,7 @@ class FITC(SparseGP):
         self._dpsi1_dX_jkj = 0
         self._dpsi1_dtheta_jkj = 0
 
-        for i, V_n, alpha_n, gamma_n, gamma_k in zip(range(self.N), self.V_star, alpha, gamma_2, gamma_3):
+        for i, V_n, alpha_n, gamma_n, gamma_k in zip(range(self.num_data), self.V_star, alpha, gamma_2, gamma_3):
             K_pp_K = np.dot(Kmmipsi1[:, i:(i + 1)], Kmmipsi1[:, i:(i + 1)].T)
 
             # Diag_dpsi1 = Diag_dA_dpsi1: yT*beta_star*y + Diag_dC_dpsi1 +Diag_dD_dpsi1
@@ -155,7 +155,7 @@ class FITC(SparseGP):
 
     def log_likelihood(self):
         """ Compute the (lower bound on the) log marginal likelihood """
-        A = -0.5 * self.N * self.input_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.beta_star)) - 0.5 * np.sum(self.V_star * self.likelihood.Y)
+        A = -0.5 * self.num_data * self.input_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.beta_star)) - 0.5 * np.sum(self.V_star * self.likelihood.Y)
         C = -self.input_dim * (np.sum(np.log(np.diag(self.LB))))
         D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
         return A + C + D
diff --git a/GPy/models/fitc_classification.py b/GPy/models/fitc_classification.py
index 7de73636..4ff441c6 100644
--- a/GPy/models/fitc_classification.py
+++ b/GPy/models/fitc_classification.py
@@ -16,7 +16,7 @@ class FITCClassification(FITC):
 
     :param X: input observations
     :param Y: observed values
-    :param likelihood: a GPy likelihood, defaults to binomial with probit link_function
+    :param likelihood: a GPy likelihood, defaults to Binomial with probit link_function
     :param kernel: a GPy kernel, defaults to rbf+white
     :param normalize_X:  whether to normalize the input data before computing (predictions will be in original scales)
     :type normalize_X: False|True
diff --git a/GPy/models/generalized_fitc.py b/GPy/models/generalized_fitc.py
index 1960d0e5..70fedcbc 100644
--- a/GPy/models/generalized_fitc.py
+++ b/GPy/models/generalized_fitc.py
@@ -101,9 +101,9 @@ class GeneralizedFITC(SparseGP):
             self.mu = self.w + np.dot(self.P, self.Gamma)
 
             # Remove extra term from dL_dpsi1
-            self.dL_dpsi1 -= mdot(self.Lmi.T, Lmipsi1 * self.likelihood.precision.flatten().reshape(1, self.N))
-            # self.Kmmi, Lm, Lmi, Kmm_logdet = pdinv(self.Kmm)
-            # self.dL_dpsi1 -= mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB
+            self.dL_dpsi1 -= mdot(self.Lmi.T,Lmipsi1 * self.likelihood.precision.flatten().reshape(1,self.num_data))
+            #self.Kmmi, Lm, Lmi, Kmm_logdet = pdinv(self.Kmm)
+            #self.dL_dpsi1 -= mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.num_data)) #dB
 
             #########333333
             # self.Bi, self.LB, self.LBi, self.B_logdet = pdinv(self.B)
@@ -130,11 +130,11 @@ class GeneralizedFITC(SparseGP):
             raise NotImplementedError, "heteroscedastic derivates not implemented"
         else:
             raise NotImplementedError, "homoscedastic derivatives not implemented"
-            # likelihood is not heterscedatic
-            # self.partial_for_likelihood =   - 0.5 * self.N*self.input_dim*self.likelihood.precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self.likelihood.precision**2
-            # self.partial_for_likelihood += 0.5 * self.input_dim * trace_dot(self.Bi,self.A)*self.likelihood.precision
-            # self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1))
-        # TODO partial derivative vector for the likelihood not implemented
+            #likelihood is not heterscedatic
+            #self.partial_for_likelihood =   - 0.5 * self.num_data*self.input_dim*self.likelihood.precision + 0.5 * np.sum(np.square(self.likelihood.Y))*self.likelihood.precision**2
+            #self.partial_for_likelihood += 0.5 * self.input_dim * trace_dot(self.Bi,self.A)*self.likelihood.precision
+            #self.partial_for_likelihood += self.likelihood.precision*(0.5*trace_dot(self.psi2_beta_scaled,self.E*sf2) - np.trace(self.Cpsi1VVpsi1))
+        #TODO partial derivative vector for the likelihood not implemented
 
     def dL_dtheta(self):
         """
@@ -153,14 +153,14 @@ class GeneralizedFITC(SparseGP):
         """ Compute the (lower bound on the) log marginal likelihood """
         sf2 = self.scale_factor ** 2
         if self.likelihood.is_heteroscedastic:
-            A = -0.5 * self.N * self.input_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.V * self.likelihood.Y)
+            A = -0.5*self.num_data*self.input_dim*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y)
         else:
-            A = -0.5 * self.N * self.input_dim * (np.log(2.*np.pi) + np.log(self.likelihood._variance)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
-        C = -self.input_dim * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.num_inducing * np.log(sf2))
-        # C = -0.5*self.input_dim * (self.B_logdet + self.num_inducing*np.log(sf2))
-        D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
-        # self.Cpsi1VVpsi1 = np.dot(self.Cpsi1V,self.psi1V.T)
-        # D_ = 0.5*np.trace(self.Cpsi1VVpsi1)
+            A = -0.5*self.num_data*self.input_dim*(np.log(2.*np.pi) + np.log(self.likelihood._variance)) -0.5*self.likelihood.precision*self.likelihood.trYYT
+        C = -self.input_dim * (np.sum(np.log(np.diag(self.LB))) + 0.5*self.num_inducing*np.log(sf2))
+        #C = -0.5*self.input_dim * (self.B_logdet + self.num_inducing*np.log(sf2))
+        D = 0.5*np.sum(np.square(self._LBi_Lmi_psi1V))
+        #self.Cpsi1VVpsi1 = np.dot(self.Cpsi1V,self.psi1V.T)
+        #D_ = 0.5*np.trace(self.Cpsi1VVpsi1)
         return A + C + D
 
     def _raw_predict(self, Xnew, which_parts, full_cov=False):
@@ -180,7 +180,7 @@ class GeneralizedFITC(SparseGP):
             #   = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
             #   = I - V.T * V
             U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)
-            V, info = linalg.lapack.dtrtrs(U, self.RPT0.T, lower=1)
+            V, info = linalg.flapack.dtrtrs(U, self.RPT0.T, lower=1)
             C = np.eye(self.num_inducing) - np.dot(V.T, V)
             mu_u = np.dot(C, self.RPT0) * (1. / self.Diag0[None, :])
             # self.C = C
diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py
index 4bd37022..e602a59a 100644
--- a/GPy/models/gplvm.py
+++ b/GPy/models/gplvm.py
@@ -42,13 +42,13 @@ class GPLVM(GP):
             return np.random.randn(Y.shape[0], input_dim)
 
     def _get_param_names(self):
-        return sum([['X_%i_%i'%(n,q) for q in range(self.input_dim)] for n in range(self.N)],[]) + GP._get_param_names(self)
+        return sum([['X_%i_%i'%(n,q) for q in range(self.input_dim)] for n in range(self.num_data)],[]) + GP._get_param_names(self)
 
     def _get_params(self):
         return np.hstack((self.X.flatten(), GP._get_params(self)))
 
     def _set_params(self,x):
-        self.X = x[:self.N*self.input_dim].reshape(self.N,self.input_dim).copy()
+        self.X = x[:self.num_data*self.input_dim].reshape(self.num_data,self.input_dim).copy()
         GP._set_params(self, x[self.X.size:])
 
     def _log_likelihood_gradients(self):
diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py
index 2c4d090f..b078fd27 100644
--- a/GPy/models/mrd.py
+++ b/GPy/models/mrd.py
@@ -74,8 +74,8 @@ class MRD(Model):
         nparams = numpy.array([0] + [SparseGP._get_params(g).size - g.Z.size for g in self.bgplvms])
         self.nparams = nparams.cumsum()
 
-        self.N = self.gref.N
-        self.NQ = self.N * self.input_dim
+        self.num_data = self.gref.num_data
+        self.NQ = self.num_data * self.input_dim
         self.MQ = self.num_inducing * self.input_dim
 
         Model.__init__(self) # @UndefinedVariable
@@ -142,8 +142,8 @@ class MRD(Model):
         self._init_Z(initz, self.X)
 
     def _get_param_names(self):
-        # X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
-        # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
+        # X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
+        # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
         n1 = self.gref._get_param_names()
         n1var = n1[:self.NQ * 2 + self.MQ]
         map_names = lambda ns, name: map(lambda x: "{1}_{0}".format(*x),
@@ -169,8 +169,8 @@ class MRD(Model):
         return params
 
 #     def _set_var_params(self, g, X, X_var, Z):
-#         g.X = X.reshape(self.N, self.input_dim)
-#         g.X_variance = X_var.reshape(self.N, self.input_dim)
+#         g.X = X.reshape(self.num_data, self.input_dim)
+#         g.X_variance = X_var.reshape(self.num_data, self.input_dim)
 #         g.Z = Z.reshape(self.num_inducing, self.input_dim)
 #
 #     def _set_kern_params(self, g, p):
diff --git a/GPy/models/sparse_gplvm.py b/GPy/models/sparse_gplvm.py
index 63368875..76fe65f1 100644
--- a/GPy/models/sparse_gplvm.py
+++ b/GPy/models/sparse_gplvm.py
@@ -28,14 +28,14 @@ class SparseGPLVM(SparseGPRegression, GPLVM):
         SparseGPRegression.__init__(self, X, Y, kernel=kernel, num_inducing=num_inducing)
 
     def _get_param_names(self):
-        return (sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
+        return (sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], [])
                 + SparseGPRegression._get_param_names(self))
 
     def _get_params(self):
         return np.hstack((self.X.flatten(), SparseGPRegression._get_params(self)))
 
     def _set_params(self, x):
-        self.X = x[:self.X.size].reshape(self.N, self.input_dim).copy()
+        self.X = x[:self.X.size].reshape(self.num_data, self.input_dim).copy()
         SparseGPRegression._set_params(self, x[self.X.size:])
 
     def log_likelihood(self):
diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py
index 1507381c..7ee9ef40 100644
--- a/GPy/testing/unit_tests.py
+++ b/GPy/testing/unit_tests.py
@@ -5,6 +5,7 @@
 import unittest
 import numpy as np
 import GPy
+from GPy.likelihoods.likelihood_functions import Binomial
 
 class GradientTests(unittest.TestCase):
     def setUp(self):
@@ -143,7 +144,7 @@ class GradientTests(unittest.TestCase):
 
     def test_GPLVM_rbf_bias_white_kern_2D(self):
         """ Testing GPLVM with rbf + bias and white kernel """
-        N, input_dim = 50, 1
+        N, input_dim, D = 50, 1, 2
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim, 0.5, 0.9 * np.ones((1,))) + GPy.kern.bias(input_dim, 0.1) + GPy.kern.white(input_dim, 0.05)
         K = k.K(X)
@@ -154,7 +155,7 @@ class GradientTests(unittest.TestCase):
 
     def test_GPLVM_rbf_linear_white_kern_2D(self):
         """ Testing GPLVM with rbf + bias and white kernel """
-        N, input_dim = 50, 1
+        N, input_dim, D = 50, 1, 2
         X = np.random.rand(N, input_dim)
         k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim, 0.1) + GPy.kern.white(input_dim, 0.05)
         K = k.K(X)
@@ -193,12 +194,12 @@ class GradientTests(unittest.TestCase):
         N = 20
         X = np.hstack([np.random.rand(N / 2) + 1, np.random.rand(N / 2) - 1])[:, None]
         k = GPy.kern.rbf(1) + GPy.kern.white(1)
-        Y = np.hstack([np.ones(N / 2), -np.ones(N / 2)])[:, None]
+        Y = np.hstack([np.ones(N/2),-np.ones(N/2)])[:,None]
 
         distribution = GPy.likelihoods.likelihood_functions.Binomial()
         likelihood = GPy.likelihoods.EP(Y, distribution)
-        # likelihood = GPy.inference.likelihoods.binomial(Y)
-        m = GPy.models.generalized_fitc(X, likelihood, k, inducing=4)
+        #likelihood = GPy.inference.likelihoods.Binomial(Y)
+        m = GPy.models.generalized_FITC(X,likelihood,k,inducing=4)
         m.constrain_positive('(var|len)')
         m.approximate_likelihood()
         self.assertTrue(m.checkgrad())