diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 60726b1d..f7d9cda4 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -13,7 +13,7 @@ default_seed = np.random.seed(123344)
 
 def BGPLVM(seed=default_seed):
     N = 10
-    M = 3
+    num_inducing = 3
     Q = 2
     D = 4
     # generate GPLVM-like data
@@ -27,7 +27,7 @@ def BGPLVM(seed=default_seed):
     # k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
     # k = GPy.kern.rbf(Q, ARD = False)  + GPy.kern.white(Q, 0.00001)
 
-    m = GPy.models.BayesianGPLVM(Y, Q, kernel=k, M=M)
+    m = GPy.models.BayesianGPLVM(Y, Q, kernel=k, num_inducing=num_inducing)
     m.constrain_positive('(rbf|bias|noise|white|S)')
     # m.constrain_fixed('S', 1)
 
@@ -63,7 +63,7 @@ def GPLVM_oil_100(optimize=True):
     m.plot_latent(labels=m.data_labels)
     return m
 
-def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
+def swiss_roll(optimize=True, N=1000, num_inducing=15, Q=4, sigma=.2, plot=False):
     from GPy.util.datasets import swiss_roll
     from GPy.core.transformations import logexp_clipped
 
@@ -101,11 +101,11 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
     S = (var * np.ones_like(X) + np.clip(np.random.randn(N, Q) * var ** 2,
                                          - (1 - var),
                                          (1 - var))) + .001
-    Z = np.random.permutation(X)[:M]
+    Z = np.random.permutation(X)[:num_inducing]
 
     kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
 
-    m = BayesianGPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
+    m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
     m.data_colors = c
     m.data_t = t
 
@@ -118,7 +118,7 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
         m.optimize('scg', messages=1)
     return m
 
-def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k):
+def BGPLVM_oil(optimize=True, N=100, Q=5, num_inducing=25, max_f_eval=4e3, plot=False, **k):
     np.random.seed(0)
     data = GPy.util.datasets.oil()
     from GPy.core.transformations import logexp_clipped
@@ -129,7 +129,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
     Yn = Y - Y.mean(0)
     Yn /= Yn.std(0)
 
-    m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, M=M, **k)
+    m = GPy.models.BayesianGPLVM(Yn, Q, kernel=kernel, num_inducing=num_inducing, **k)
     m.data_labels = data['Y'][:N].argmax(axis=1)
 
     # m.constrain('variance|leng', logexp_clipped())
@@ -168,7 +168,7 @@ def oil_100():
 
 
 
-def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
+def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
     x = np.linspace(0, 4 * np.pi, N)[:, None]
     s1 = np.vectorize(lambda x: np.sin(x))
     s2 = np.vectorize(lambda x: np.cos(x))
@@ -228,13 +228,13 @@ def bgplvm_simulation_matlab_compare():
     Y = sim_data['Y']
     S = sim_data['S']
     mu = sim_data['mu']
-    M, [_, Q] = 3, mu.shape
+    num_inducing, [_, Q] = 3, mu.shape
 
     from GPy.models import mrd
     from GPy import kern
     reload(mrd); reload(kern)
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
-    m = BayesianGPLVM(Y, Q, init="PCA", M=M, kernel=k,
+    m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k,
 #                        X=mu,
 #                        X_variance=S,
                        _debug=False)
@@ -248,8 +248,8 @@ def bgplvm_simulation(optimize='scg',
                       plot=True,
                       max_f_eval=2e4):
 #     from GPy.core.transformations import logexp_clipped
-    D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
-    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
+    D1, D2, D3, N, num_inducing, Q = 15, 8, 8, 100, 3, 5
+    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot)
 
     from GPy.models import mrd
     from GPy import kern
@@ -259,7 +259,7 @@ def bgplvm_simulation(optimize='scg',
     Y = Ylist[0]
 
     k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
-    m = BayesianGPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
+    m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k, _debug=True)
     # m.constrain('variance|noise', logexp_clipped())
     m.ensure_default_constraints()
     m['noise'] = Y.var() / 100.
@@ -276,8 +276,8 @@ def bgplvm_simulation(optimize='scg',
     return m
 
 def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
-    D1, D2, D3, N, M, Q = 150, 200, 400, 500, 3, 7
-    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
+    D1, D2, D3, N, num_inducing, Q = 150, 200, 400, 500, 3, 7
+    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
 
     from GPy.models import mrd
     from GPy import kern
@@ -285,7 +285,7 @@ def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
     reload(mrd); reload(kern)
 
     k = kern.linear(Q, [.05] * Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
-    m = mrd.MRD(Ylist, input_dim=Q, M=M, kernels=k, initx="", initz='permute', **kw)
+    m = mrd.MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernels=k, initx="", initz='permute', **kw)
 
     for i, Y in enumerate(Ylist):
         m['{}_noise'.format(i + 1)] = Y.var() / 100.
@@ -313,7 +313,7 @@ def brendan_faces():
     Yn /= Yn.std()
 
     m = GPy.models.GPLVM(Yn, Q)
-    # m = GPy.models.BayesianGPLVM(Yn, Q, M=100)
+    # m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=100)
 
     # optimize
     m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
@@ -377,16 +377,16 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True):
 #     X /= X.std(axis=0)
 #
 #     Q = 10
-#     M = 30
+#     num_inducing = 30
 #
 #     kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
-#     m = GPy.models.BayesianGPLVM(X, Q, kernel=kernel, M=M)
+#     m = GPy.models.BayesianGPLVM(X, Q, kernel=kernel, num_inducing=num_inducing)
 #     # m.scale_factor = 100.0
 #     m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)')
 #     from sklearn import cluster
-#     km = cluster.KMeans(M, verbose=10)
+#     km = cluster.KMeans(num_inducing, verbose=10)
 #     Z = km.fit(m.X).cluster_centers_
-#     # Z = GPy.util.misc.kmm_init(m.X, M)
+#     # Z = GPy.util.misc.kmm_init(m.X, num_inducing)
 #     m.set('iip', Z)
 #     m.set('bias', 1e-4)
 #     # optimize
diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index be3a71bd..64a2d12c 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -151,8 +151,8 @@ def coregionalisation_sparse(optim_iters=100):
     Y2 = -np.sin(X2) + np.random.randn(*X2.shape)*0.05
     Y = np.vstack((Y1,Y2))
 
-    M = 40
-    Z = np.hstack((np.random.rand(M,1)*8,np.random.randint(0,2,M)[:,None]))
+    num_inducing = 40
+    Z = np.hstack((np.random.rand(num_inducing,1)*8,np.random.randint(0,2,num_inducing)[:,None]))
 
     k1 = GPy.kern.rbf(1)
     k2 = GPy.kern.Coregionalise(2,2)
@@ -261,7 +261,7 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf):
 
     return np.array(lls)
 
-def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
+def sparse_GP_regression_1D(N = 400, num_inducing = 5, optim_iters=100):
     """Run a 1D example of a sparse GP regression."""
     # sample inputs and outputs
     X = np.random.uniform(-3.,3.,(N,1))
@@ -271,7 +271,7 @@ def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
     noise = GPy.kern.white(1)
     kernel = rbf + noise
     # create simple GP model
-    m = GPy.models.SparseGPRegression(X, Y, kernel, M=M)
+    m = GPy.models.SparseGPRegression(X, Y, kernel, num_inducing=num_inducing)
 
     m.ensure_default_constraints()
 
@@ -280,7 +280,7 @@ def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100):
     m.plot()
     return m
 
-def sparse_GP_regression_2D(N = 400, M = 50, optim_iters=100):
+def sparse_GP_regression_2D(N = 400, num_inducing = 50, optim_iters=100):
     """Run a 2D example of a sparse GP regression."""
     X = np.random.uniform(-3.,3.,(N,2))
     Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(N,1)*0.05
@@ -291,7 +291,7 @@ def sparse_GP_regression_2D(N = 400, M = 50, optim_iters=100):
     kernel = rbf + noise
 
     # create simple GP model
-    m = GPy.models.SparseGPRegression(X,Y,kernel, M = M)
+    m = GPy.models.SparseGPRegression(X,Y,kernel, num_inducing = num_inducing)
 
     # contrain all parameters to be positive (but not inducing inputs)
     m.ensure_default_constraints()
diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py
index 7fe922fd..329d1015 100644
--- a/GPy/kern/coregionalise.py
+++ b/GPy/kern/coregionalise.py
@@ -69,14 +69,14 @@ class Coregionalise(kernpart):
         else:
             index2 = np.asarray(index2,dtype=np.int)
             code="""
-            for(int i=0;i<M; i++){
+            for(int i=0;i<num_inducing; i++){
               for(int j=0; j<N; j++){
-                  target[i+j*M] += B[Nout*index[j]+index2[i]];
+                  target[i+j*num_inducing] += B[Nout*index[j]+index2[i]];
                 }
               }
             """
-            N,M,B,Nout = index.size,index2.size, self.B, self.Nout
-            weave.inline(code,['target','index','index2','N','M','B','Nout'])
+            N,num_inducing,B,Nout = index.size,index2.size, self.B, self.Nout
+            weave.inline(code,['target','index','index2','N','num_inducing','B','Nout'])
 
 
     def Kdiag(self,index,target):
@@ -91,14 +91,14 @@ class Coregionalise(kernpart):
             index2 = np.asarray(index2,dtype=np.int)
 
         code="""
-        for(int i=0; i<M; i++){
+        for(int i=0; i<num_inducing; i++){
           for(int j=0; j<N; j++){
-            dL_dK_small[index[j] + Nout*index2[i]] += dL_dK[i+j*M];
+            dL_dK_small[index[j] + Nout*index2[i]] += dL_dK[i+j*num_inducing];
           }
         }
         """
-        N, M, Nout = index.size, index2.size, self.Nout
-        weave.inline(code, ['N','M','Nout','dL_dK','dL_dK_small','index','index2'])
+        N, num_inducing, Nout = index.size, index2.size, self.Nout
+        weave.inline(code, ['N','num_inducing','Nout','dL_dK','dL_dK_small','index','index2'])
 
         dkappa = np.diag(dL_dK_small)
         dL_dK_small += dL_dK_small.T
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 05af1997..efd49020 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -223,11 +223,11 @@ class kern(parameterised):
     def dK_dtheta(self, dL_dK, X, X2=None):
         """
         :param dL_dK: An array of dL_dK derivaties, dL_dK
-        :type dL_dK: Np.ndarray (N x M)
+        :type dL_dK: Np.ndarray (N x num_inducing)
         :param X: Observed data inputs
         :type X: np.ndarray (N x input_dim)
         :param X2: Observed dara inputs (optional, defaults to X)
-        :type X2: np.ndarray (M x input_dim)
+        :type X2: np.ndarray (num_inducing x input_dim)
         """
         assert X.shape[1] == self.input_dim
         target = np.zeros(self.Nparam)
@@ -300,16 +300,16 @@ class kern(parameterised):
         return target
 
     def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
-        """return shapes are N,M,input_dim"""
+        """return shapes are N,num_inducing,input_dim"""
         target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
         [p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
         return target_mu, target_S
 
     def psi2(self, Z, mu, S):
         """
-        :param Z: np.ndarray of inducing inputs (M x input_dim)
+        :param Z: np.ndarray of inducing inputs (num_inducing x input_dim)
         :param mu, S: np.ndarrays of means and variances (each N x input_dim)
-        :returns psi2: np.ndarray (N,M,M)
+        :returns psi2: np.ndarray (N,num_inducing,num_inducing)
         """
         target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
         [p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
@@ -328,7 +328,7 @@ class kern(parameterised):
 
             prod = np.multiply(tmp1, tmp2)
             crossterms += prod[:,:,None] + prod[:, None, :]
-            
+
         target += crossterms
         return target
 
diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py
index 425b81bb..40dd2007 100644
--- a/GPy/kern/linear.py
+++ b/GPy/kern/linear.py
@@ -138,7 +138,7 @@ class linear(kernpart):
 
     def psi2(self, Z, mu, S, target):
         """
-        returns N,M,M matrix
+        returns N,num_inducing,num_inducing matrix
         """
         self._psi_computations(Z, mu, S)
 #         psi2_old = self.ZZ * np.square(self.variances) * self.mu2_S[:, None, None, :]
@@ -168,7 +168,7 @@ class linear(kernpart):
             target += tmp.sum()
 
     def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
-        """Think N,M,M,input_dim """
+        """Think N,num_inducing,num_inducing,input_dim """
         self._psi_computations(Z, mu, S)
         AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
         AZZA = AZZA + AZZA.swapaxes(1, 2)
@@ -184,7 +184,7 @@ class linear(kernpart):
         double factor,tmp;
         #pragma omp parallel for private(m,mm,q,qq,factor,tmp)
         for(n=0;n<N;n++){
-          for(m=0;m<M;m++){
+          for(m=0;m<num_inducing;m++){
             for(mm=0;mm<=m;mm++){
               //add in a factor of 2 for the off-diagonal terms (and then count them only once)
               if(m==mm)
@@ -215,9 +215,9 @@ class linear(kernpart):
                          'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
                          'extra_link_args'   : ['-lgomp']}
 
-        N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+        N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
+                     arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
 
 
@@ -231,9 +231,9 @@ class linear(kernpart):
         code="""
         int n,m,mm,q;
         #pragma omp parallel for private(n,mm,q)
-        for(m=0;m<M;m++){
+        for(m=0;m<num_inducing;m++){
           for(q=0;q<input_dim;q++){
-            for(mm=0;mm<M;mm++){
+            for(mm=0;mm<num_inducing;mm++){
               for(n=0;n<N;n++){
                 target(m,q) += dL_dpsi2(n,m,mm)*AZA(n,mm,q);
               }
@@ -249,9 +249,9 @@ class linear(kernpart):
                          'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
                          'extra_link_args'   : ['-lgomp']}
 
-        N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
+        N,num_inducing,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','input_dim','AZA','target','dL_dpsi2'],
+                     arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
 
 
@@ -278,7 +278,7 @@ class linear(kernpart):
         muS_changed = not (np.array_equal(mu, self._mu) and np.array_equal(S, self._S))
         if Zv_changed:
             # Z has changed, compute Z specific stuff
-            # self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,input_dim
+            # self.ZZ = Z[:,None,:]*Z[None,:,:] # num_inducing,num_inducing,input_dim
 #             self.ZZ = np.empty((Z.shape[0], Z.shape[0], Z.shape[1]), order='F')
 #             [tdot(Z[:, i:i + 1], self.ZZ[:, :, i].T) for i in xrange(Z.shape[1])]
             self.ZA = Z * self.variances
@@ -291,5 +291,5 @@ class linear(kernpart):
             self.inner[:, diag_indices[0], diag_indices[1]] += S
             self._mu, self._S = mu.copy(), S.copy()
         if Zv_changed or muS_changed:
-            self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1)  # NOTE: self.ZAinner \in [M x N x input_dim]!
+            self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1)  # NOTE: self.ZAinner \in [num_inducing x N x input_dim]!
             self._psi2 = np.dot(self.ZAinner, self.ZA.T)
diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py
index 95684a02..fd85f45a 100644
--- a/GPy/kern/periodic_Matern32.py
+++ b/GPy/kern/periodic_Matern32.py
@@ -64,7 +64,7 @@ class periodic_Matern32(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
-    
+
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -113,7 +113,7 @@ class periodic_Matern32(kernpart):
 
     @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
+        """derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
         if X2 is None: X2 = X
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py
index 1e55ab62..226e7f3b 100644
--- a/GPy/kern/periodic_Matern52.py
+++ b/GPy/kern/periodic_Matern52.py
@@ -64,7 +64,7 @@ class periodic_Matern52(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
-    
+
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -115,7 +115,7 @@ class periodic_Matern52(kernpart):
 
     @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
+        """derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
         if X2 is None: X2 = X
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py
index 50575ca9..d256cb4a 100644
--- a/GPy/kern/periodic_exponential.py
+++ b/GPy/kern/periodic_exponential.py
@@ -64,7 +64,7 @@ class periodic_exponential(kernpart):
     def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscale,self.period))
-    
+
     def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==3
@@ -111,7 +111,7 @@ class periodic_exponential(kernpart):
 
     @silence_errors
     def dK_dtheta(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)"""
+        """derivative of the covariance matrix with respect to the parameters (shape is Nxnum_inducingxNparam)"""
         if X2 is None: X2 = X
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index 39cf878e..bfa20f1d 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -110,7 +110,7 @@ class rbf(kernpart):
                   target(q+1) += var_len3(q)*tmp;
                 }
                 """
-                N, M, input_dim = X.shape[0], X.shape[0], self.input_dim
+                N, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
             else:
                 code = """
                 int q,i,j;
@@ -118,16 +118,16 @@ class rbf(kernpart):
                 for(q=0; q<input_dim; q++){
                   tmp = 0;
                   for(i=0; i<N; i++){
-                    for(j=0; j<M; j++){
+                    for(j=0; j<num_inducing; j++){
                       tmp += (X(i,q)-X2(j,q))*(X(i,q)-X2(j,q))*dvardLdK(i,j);
                     }
                   }
                   target(q+1) += var_len3(q)*tmp;
                 }
                 """
-                N, M, input_dim = X.shape[0], X2.shape[0], self.input_dim
+                N, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
             #[np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)]
-            weave.inline(code, arg_names=['N','M','input_dim','X','X2','target','dvardLdK','var_len3'],
+            weave.inline(code, arg_names=['N','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3'],
                  type_converters=weave.converters.blitz,**self.weave_options)
         else:
             target[1] += (self.variance/self.lengthscale)*np.sum(self._K_dvar*self._K_dist2*dL_dK)
@@ -191,7 +191,7 @@ class rbf(kernpart):
         target += self._psi2
 
     def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target):
-        """Shape N,M,M,Ntheta"""
+        """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)
@@ -205,13 +205,13 @@ class rbf(kernpart):
 
     def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target):
         self._psi_computations(Z,mu,S)
-        term1 = self._psi2_Zdist/self.lengthscale2 # M, M, input_dim
-        term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, input_dim
+        term1 = self._psi2_Zdist/self.lengthscale2 # num_inducing, num_inducing, input_dim
+        term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, num_inducing, num_inducing, input_dim
         dZ = self._psi2[:,:,:,None] * (term1[None] + term2)
         target += (dL_dpsi2[:,:,:,None]*dZ).sum(0).sum(0)
 
     def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S):
-        """Think N,M,M,input_dim """
+        """Think N,num_inducing,num_inducing,input_dim """
         self._psi_computations(Z,mu,S)
         tmp = self._psi2[:,:,:,None]/self.lengthscale2/self._psi2_denom
         target_mu += -2.*(dL_dpsi2[:,:,:,None]*tmp*self._psi2_mudist).sum(1).sum(1)
@@ -242,9 +242,9 @@ class rbf(kernpart):
         #here are the "statistics" for psi1 and psi2
         if not np.array_equal(Z, self._Z):
             #Z has changed, compute Z specific stuff
-            self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,input_dim
-            self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,input_dim
-            self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,input_dim
+            self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # num_inducing,num_inducing,input_dim
+            self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # num_inducing,num_inducing,input_dim
+            self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # num_inducing,num_inducing,input_dim
             self._Z = Z
 
         if not (np.array_equal(Z, self._Z) and np.array_equal(mu, self._mu) and np.array_equal(S, self._S)):
@@ -258,24 +258,24 @@ class rbf(kernpart):
             self._psi1 = self.variance*np.exp(self._psi1_exponent)
 
             #psi2
-            self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,input_dim
+            self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,num_inducing,num_inducing,input_dim
             self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu,self._psi2_Zhat)
-            #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,input_dim
+            #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,num_inducing,num_inducing,input_dim
             #self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom)
-            #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M
-            self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
+            #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,num_inducing,num_inducing
+            self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,num_inducing,num_inducing
 
             #store matrices for caching
             self._Z, self._mu, self._S = Z, mu,S
 
     def weave_psi2(self,mu,Zhat):
         N,input_dim = mu.shape
-        M = Zhat.shape[0]
+        num_inducing = Zhat.shape[0]
 
-        mudist = np.empty((N,M,M,input_dim))
-        mudist_sq = np.empty((N,M,M,input_dim))
-        psi2_exponent = np.zeros((N,M,M))
-        psi2 = np.empty((N,M,M))
+        mudist = np.empty((N,num_inducing,num_inducing,input_dim))
+        mudist_sq = np.empty((N,num_inducing,num_inducing,input_dim))
+        psi2_exponent = np.zeros((N,num_inducing,num_inducing))
+        psi2 = np.empty((N,num_inducing,num_inducing))
 
         psi2_Zdist_sq = self._psi2_Zdist_sq
         _psi2_denom = self._psi2_denom.squeeze().reshape(N,self.input_dim)
@@ -290,7 +290,7 @@ class rbf(kernpart):
 
         #pragma omp parallel for private(tmp)
         for (int n=0; n<N; n++){
-            for (int m=0; m<M; m++){
+            for (int m=0; m<num_inducing; m++){
                for (int mm=0; mm<(m+1); mm++){
                    for (int q=0; q<input_dim; q++){
                        //compute mudist
@@ -325,7 +325,7 @@ class rbf(kernpart):
         #include <math.h>
         """
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','input_dim','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
+                     arg_names=['N','num_inducing','input_dim','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
                      type_converters=weave.converters.blitz,**self.weave_options)
 
         return mudist,mudist_sq, psi2_exponent, psi2
diff --git a/GPy/kern/sympykern.py b/GPy/kern/sympykern.py
index db3cc976..0ec86acf 100644
--- a/GPy/kern/sympykern.py
+++ b/GPy/kern/sympykern.py
@@ -122,12 +122,12 @@ class spkern(kernpart):
         int i;
         int j;
         int N = target_array->dimensions[0];
-        int M = target_array->dimensions[1];
+        int num_inducing = target_array->dimensions[1];
         int D = X_array->dimensions[1];
         //#pragma omp parallel for private(j)
         for (i=0;i<N;i++){
-            for (j=0;j<M;j++){
-                target[i*M+j] = k(%s);
+            for (j=0;j<num_inducing;j++){
+                target[i*num_inducing+j] = k(%s);
             }
         }
         %s
@@ -149,17 +149,17 @@ class spkern(kernpart):
         """%(diag_arglist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
 
         #here's some code to compute gradients
-        funclist = '\n'.join([' '*16 + 'target[%i] += partial[i*M+j]*dk_d%s(%s);'%(i,theta.name,arglist) for i,theta in  enumerate(self._sp_theta)])
+        funclist = '\n'.join([' '*16 + 'target[%i] += partial[i*num_inducing+j]*dk_d%s(%s);'%(i,theta.name,arglist) for i,theta in  enumerate(self._sp_theta)])
         self._dK_dtheta_code =\
         """
         int i;
         int j;
         int N = partial_array->dimensions[0];
-        int M = partial_array->dimensions[1];
+        int num_inducing = partial_array->dimensions[1];
         int D = X_array->dimensions[1];
         //#pragma omp parallel for private(j)
         for (i=0;i<N;i++){
-            for (j=0;j<M;j++){
+            for (j=0;j<num_inducing;j++){
 %s
             }
         }
@@ -169,7 +169,7 @@ class spkern(kernpart):
         #here's some code to compute gradients for Kdiag TODO: thius is yucky.
         diag_funclist = re.sub('Z','X',funclist,count=0)
         diag_funclist = re.sub('j','i',diag_funclist)
-        diag_funclist = re.sub('partial\[i\*M\+i\]','partial[i]',diag_funclist)
+        diag_funclist = re.sub('partial\[i\*num_inducing\+i\]','partial[i]',diag_funclist)
         self._dKdiag_dtheta_code =\
         """
         int i;
@@ -182,17 +182,17 @@ class spkern(kernpart):
         """%(diag_funclist,"/*"+str(self._sp_k)+"*/") #adding a string representation forces recompile when needed
 
         #Here's some code to do gradients wrt x
-        gradient_funcs = "\n".join(["target[i*D+%i] += partial[i*M+j]*dk_dx%i(%s);"%(q,q,arglist) for q in range(self.D)])
+        gradient_funcs = "\n".join(["target[i*D+%i] += partial[i*num_inducing+j]*dk_dx%i(%s);"%(q,q,arglist) for q in range(self.D)])
         self._dK_dX_code = \
         """
         int i;
         int j;
         int N = partial_array->dimensions[0];
-        int M = partial_array->dimensions[1];
+        int num_inducing = partial_array->dimensions[1];
         int D = X_array->dimensions[1];
         //#pragma omp parallel for private(j)
         for (i=0;i<N; i++){
-            for (j=0; j<M; j++){
+            for (j=0; j<num_inducing; j++){
                 %s
                 //if(isnan(target[i*D+2])){printf("%%f\\n",dk_dx2(X[i*D+0], X[i*D+1], X[i*D+2], Z[j*D+0], Z[j*D+1], Z[j*D+2], param[0], param[1], param[2], param[3], param[4], param[5]));}
                 //if(isnan(target[i*D+2])){printf("%%f,%%f,%%i,%%i\\n", X[i*D+2], Z[j*D+2],i,j);}
@@ -208,7 +208,7 @@ class spkern(kernpart):
         int i;
         int j;
         int N = partial_array->dimensions[0];
-        int M = 0;
+        int num_inducing = 0;
         int D = X_array->dimensions[1];
         for (i=0;i<N; i++){
             j = i;
diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py
index 6409ed58..70f30fd7 100644
--- a/GPy/likelihoods/ep.py
+++ b/GPy/likelihoods/ep.py
@@ -139,7 +139,7 @@ class EP(likelihood):
         The expectation-propagation algorithm with sparse pseudo-input.
         For nomenclature see ... 2013.
         """
-        M = Kmm.shape[0]
+        num_inducing = Kmm.shape[0]
 
         #TODO: this doesn't work with uncertain inputs!
 
@@ -235,7 +235,7 @@ class EP(likelihood):
         The expectation-propagation algorithm with sparse pseudo-input.
         For nomenclature see Naish-Guzman and Holden, 2008.
         """
-        M = Kmm.shape[0]
+        num_inducing = Kmm.shape[0]
 
         """
         Prior approximation parameters:
@@ -258,7 +258,7 @@ class EP(likelihood):
         mu = w + P*Gamma
         """
         self.w = np.zeros(self.N)
-        self.Gamma = np.zeros(M)
+        self.Gamma = np.zeros(num_inducing)
         mu = np.zeros(self.N)
         P = P0.copy()
         R = R0.copy()
@@ -305,10 +305,10 @@ class EP(likelihood):
                 dtd1 = Delta_tau*Diag[i] + 1.
                 dii = Diag[i]
                 Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
-                pi_ = P[i,:].reshape(1,M)
+                pi_ = P[i,:].reshape(1,num_inducing)
                 P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
                 Rp_i = np.dot(R,pi_.T)
-                RTR = np.dot(R.T,np.dot(np.eye(M) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
+                RTR = np.dot(R.T,np.dot(np.eye(num_inducing) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
                 R = jitchol(RTR).T
                 self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1
                 self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T)
@@ -321,7 +321,7 @@ class EP(likelihood):
             Diag = Diag0 * Iplus_Dprod_i
             P = Iplus_Dprod_i[:,None] * P0
             safe_diag = np.where(Diag0 < self.tau_tilde, self.tau_tilde/(1.+Diag0*self.tau_tilde), (1. - Iplus_Dprod_i)/Diag0)
-            L = jitchol(np.eye(M) + np.dot(RPT0,safe_diag[:,None]*RPT0.T))
+            L = jitchol(np.eye(num_inducing) + np.dot(RPT0,safe_diag[:,None]*RPT0.T))
             R,info = linalg.lapack.flapack.dtrtrs(L,R0,lower=1)
             RPT = np.dot(R,P.T)
             Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py
index 7008e9a3..f4a44380 100644
--- a/GPy/models/bayesian_gplvm.py
+++ b/GPy/models/bayesian_gplvm.py
@@ -23,7 +23,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', M=10,
+    def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
                  Z=None, kernel=None, oldpsave=10, _debug=False,
                  **kwargs):
         if type(likelihood_or_Y) is np.ndarray:
@@ -39,7 +39,7 @@ class BayesianGPLVM(SparseGP, GPLVM):
             X_variance = np.clip((np.ones_like(X) * 0.5) + .01 * np.random.randn(*X.shape), 0.001, 1)
 
         if Z is None:
-            Z = np.random.permutation(X.copy())[:M]
+            Z = np.random.permutation(X.copy())[:num_inducing]
         assert Z.shape[1] == X.shape[1]
 
         if kernel is None:
diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py
index 4300f113..63dae3ce 100644
--- a/GPy/models/mrd.py
+++ b/GPy/models/mrd.py
@@ -33,7 +33,7 @@ class MRD(model):
     :param X_variance:
         Initial latent space variance
     :param init: [cooncat|single|random]
-        initialization method to use: 
+        initialization method to use:
             *concat: PCA on concatenated outputs
             *single: PCA on each output
             *random: random
@@ -44,7 +44,7 @@ class MRD(model):
     :param kernels: list of kernels or kernel shared for all BGPLVMS
     :type kernels: [GPy.kern.kern] | GPy.kern.kern | None (default)
     """
-    def __init__(self, likelihood_or_Y_list, input_dim, M=10, names=None,
+    def __init__(self, likelihood_or_Y_list, input_dim, num_inducing=10, names=None,
                  kernels=None, initx='PCA',
                  initz='permute', _debug=False, **kw):
         if names is None:
@@ -61,13 +61,13 @@ class MRD(model):
         assert not ('kernel' in kw), "pass kernels through `kernels` argument"
 
         self.input_dim = input_dim
-        self.num_inducing = M
+        self.num_inducing = num_inducing
         self._debug = _debug
 
         self._init = True
         X = self._init_X(initx, likelihood_or_Y_list)
         Z = self._init_Z(initz, X)
-        self.bgplvms = [BayesianGPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, M=self.num_inducing, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
+        self.bgplvms = [BayesianGPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, num_inducing=self.num_inducing, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
         del self._init
 
         self.gref = self.bgplvms[0]
diff --git a/GPy/models/sparse_gp_classification.py b/GPy/models/sparse_gp_classification.py
index f82de00f..9027ef07 100644
--- a/GPy/models/sparse_gp_classification.py
+++ b/GPy/models/sparse_gp_classification.py
@@ -26,7 +26,7 @@ class SparseGPClassification(SparseGP):
 
     """
 
-    def __init__(self, X, Y=None, likelihood=None, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10):
+    def __init__(self, X, Y=None, likelihood=None, kernel=None, normalize_X=False, normalize_Y=False, Z=None, num_inducing=10):
         if kernel is None:
             kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1],1e-3)
 
@@ -38,7 +38,7 @@ class SparseGPClassification(SparseGP):
                 raise Warning, 'likelihood.data and Y are different.'
 
         if Z is None:
-            i = np.random.permutation(X.shape[0])[:M]
+            i = np.random.permutation(X.shape[0])[:num_inducing]
             Z = X[i].copy()
         else:
             assert Z.shape[1]==X.shape[1]
diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py
index d7907c5f..432d6e18 100644
--- a/GPy/models/sparse_gp_regression.py
+++ b/GPy/models/sparse_gp_regression.py
@@ -26,14 +26,14 @@ class SparseGPRegression(SparseGP):
 
     """
 
-    def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10, X_variance=None):
+    def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, num_inducing=10, X_variance=None):
         # kern defaults to rbf (plus white for stability)
         if kernel is None:
             kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1], 1e-3)
 
         # Z defaults to a subset of the data
         if Z is None:
-            i = np.random.permutation(X.shape[0])[:M]
+            i = np.random.permutation(X.shape[0])[:num_inducing]
             Z = X[i].copy()
         else:
             assert Z.shape[1] == X.shape[1]
diff --git a/GPy/models/sparse_gplvm.py b/GPy/models/sparse_gplvm.py
index bec8c2da..63368875 100644
--- a/GPy/models/sparse_gplvm.py
+++ b/GPy/models/sparse_gplvm.py
@@ -23,9 +23,9 @@ class SparseGPLVM(SparseGPRegression, GPLVM):
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, input_dim, kernel=None, init='PCA', M=10):
+    def __init__(self, Y, input_dim, kernel=None, init='PCA', num_inducing=10):
         X = self.initialise_latent(init, input_dim, Y)
-        SparseGPRegression.__init__(self, X, Y, kernel=kernel, M=M)
+        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)], [])
diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py
index 3ca32660..ff558f6d 100644
--- a/GPy/testing/bgplvm_tests.py
+++ b/GPy/testing/bgplvm_tests.py
@@ -8,67 +8,67 @@ from GPy.models.bayesian_gplvm import BayesianGPLVM
 
 class BGPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         Y -= Y.mean(axis=0)
         k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
+        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_linear_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         Y -= Y.mean(axis=0)
         k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
+        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         Y -= Y.mean(axis=0)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
+        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_bias_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) +  GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         Y -= Y.mean(axis=0)
         k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
+        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     #@unittest.skip('psi2 cross terms are NotImplemented for this combination')
     def test_linear_bias_kern(self):
-        N, M, input_dim, D = 30, 5, 4, 30
+        N, num_inducing, input_dim, D = 30, 5, 4, 30
         X = np.random.rand(N, input_dim)
         k = GPy.kern.linear(input_dim) +  GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         Y -= Y.mean(axis=0)
         k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, M=M)
+        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
diff --git a/GPy/testing/gplvm_tests.py b/GPy/testing/gplvm_tests.py
index 4194698f..8c2ba9fc 100644
--- a/GPy/testing/gplvm_tests.py
+++ b/GPy/testing/gplvm_tests.py
@@ -7,7 +7,7 @@ import GPy
 
 class GPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
@@ -19,7 +19,7 @@ class GPLVMTests(unittest.TestCase):
         self.assertTrue(m.checkgrad())
 
     def test_linear_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
@@ -31,7 +31,7 @@ class GPLVMTests(unittest.TestCase):
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
diff --git a/GPy/testing/mrd_tests.py b/GPy/testing/mrd_tests.py
index df4c2723..b0137709 100644
--- a/GPy/testing/mrd_tests.py
+++ b/GPy/testing/mrd_tests.py
@@ -14,7 +14,7 @@ class MRDTests(unittest.TestCase):
 
     def test_gradients(self):
         num_m = 3
-        N, M, input_dim, D = 20, 8, 6, 20
+        N, num_inducing, input_dim, D = 20, 8, 6, 20
         X = np.random.rand(N, input_dim)
 
         k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
@@ -23,7 +23,7 @@ class MRDTests(unittest.TestCase):
         Ylist = [np.random.multivariate_normal(np.zeros(N), K, input_dim).T for _ in range(num_m)]
         likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
 
-        m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, M=M)
+        m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
 
         self.assertTrue(m.checkgrad())
diff --git a/GPy/testing/psi_stat_gradient_tests.py b/GPy/testing/psi_stat_gradient_tests.py
index 2257d16c..ce16bf9b 100644
--- a/GPy/testing/psi_stat_gradient_tests.py
+++ b/GPy/testing/psi_stat_gradient_tests.py
@@ -11,7 +11,7 @@ import itertools
 from GPy.core import model
 
 class PsiStatModel(model):
-    def __init__(self, which, X, X_variance, Z, M, kernel):
+    def __init__(self, which, X, X_variance, Z, num_inducing, kernel):
         self.which = which
         self.X = X
         self.X_variance = X_variance
@@ -64,8 +64,8 @@ class DPsiStatTest(unittest.TestCase):
 
     def testPsi0(self):
         for k in self.kernels:
-            m = PsiStatModel('psi0', X=self.X, X_variance=self.X_var, Z=self.Z,
-                         M=self.num_inducing, kernel=k)
+            m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
+                         num_inducing=self.num_inducing, kernel=k)
             try:
                 assert m.checkgrad(), "{} x psi0".format("+".join(map(lambda x: x.name, k.parts)))
             except:
@@ -74,33 +74,33 @@ class DPsiStatTest(unittest.TestCase):
 #     def testPsi1(self):
 #         for k in self.kernels:
 #             m = PsiStatModel('psi1', X=self.X, X_variance=self.X_var, Z=self.Z,
-#                      M=self.M, kernel=k)
+#                      num_inducing=self.num_inducing, kernel=k)
 #             assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
 
     def testPsi2_lin(self):
         k = self.kernels[0]
         m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
-                     M=self.num_inducing, kernel=k)
+                     num_inducing=self.num_inducing, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
     def testPsi2_lin_bia(self):
         k = self.kernels[3]
         m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
-                     M=self.num_inducing, kernel=k)
+                     num_inducing=self.num_inducing, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
     def testPsi2_rbf(self):
         k = self.kernels[1]
         m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
-                     M=self.num_inducing, kernel=k)
+                     num_inducing=self.num_inducing, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
     def testPsi2_rbf_bia(self):
         k = self.kernels[-1]
         m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
-                     M=self.num_inducing, kernel=k)
+                     num_inducing=self.num_inducing, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
     def testPsi2_bia(self):
         k = self.kernels[2]
         m = PsiStatModel('psi2', X=self.X, X_variance=self.X_var, Z=self.Z,
-                     M=self.num_inducing, kernel=k)
+                     num_inducing=self.num_inducing, kernel=k)
         assert m.checkgrad(), "{} x psi2".format("+".join(map(lambda x: x.name, k.parts)))
 
 
@@ -122,11 +122,11 @@ if __name__ == "__main__":
         numpy.random.seed(0)
         input_dim = 5
         N = 50
-        M = 10
+        num_inducing = 10
         D = 15
         X = numpy.random.randn(N, input_dim)
         X_var = .5 * numpy.ones_like(X) + .1 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
-        Z = numpy.random.permutation(X)[:M]
+        Z = numpy.random.permutation(X)[:num_inducing]
         Y = X.dot(numpy.random.randn(input_dim, D))
 #         kernel = GPy.kern.bias(input_dim)
 #
@@ -148,7 +148,7 @@ if __name__ == "__main__":
 #         m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
 #                          num_inducing=num_inducing, kernel=GPy.kern.rbf(input_dim))
         m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
-                         M=M, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
+                         num_inducing=num_inducing, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
         m3.ensure_default_constraints()
         # + GPy.kern.bias(input_dim))
 #         m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
diff --git a/GPy/testing/sparse_gplvm_tests.py b/GPy/testing/sparse_gplvm_tests.py
index 38f95730..e7f714b4 100644
--- a/GPy/testing/sparse_gplvm_tests.py
+++ b/GPy/testing/sparse_gplvm_tests.py
@@ -8,38 +8,38 @@ from GPy.models.sparse_gplvm import SparseGPLVM
 
 class sparse_GPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
+        m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     @unittest.skip('linear kernels do not have dKdiag_dX')
     def test_linear_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
+        m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, input_dim, D = 10, 3, 2, 4
+        N, num_inducing, input_dim, D = 10, 3, 2, 4
         X = np.random.rand(N, input_dim)
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
         K = k.K(X)
         Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
         k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
-        m = SparseGPLVM(Y, input_dim, kernel=k, M=M)
+        m = SparseGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())