diff --git a/GPy/core/GP.py b/GPy/core/GP.py
index 01648189..04ea7af1 100644
--- a/GPy/core/GP.py
+++ b/GPy/core/GP.py
@@ -126,7 +126,7 @@ class GP(GPBase):
         Arguments
         ---------
         :param Xnew: The points at which to make a prediction
-        :type Xnew: np.ndarray, Nnew x self.Q
+        :type Xnew: np.ndarray, Nnew x self.input_dim
         :param which_parts:  specifies which outputs kernel(s) to use in prediction
         :type which_parts: ('all', list of bools)
         :param full_cov: whether to return the folll covariance matrix, or just the diagonal
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index 26bed691..aa71b550 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -13,7 +13,7 @@ class GPBase(model.model):
     def __init__(self, X, likelihood, kernel, normalize_X=False):
         self.X = X
         assert len(self.X.shape) == 2
-        self.N, self.Q = self.X.shape
+        self.N, self.input_dim = self.X.shape
         assert isinstance(kernel, kern.kern)
         self.kern = kernel
         self.likelihood = likelihood
@@ -25,8 +25,8 @@ class GPBase(model.model):
             self._Xstd = X.std(0)[None, :]
             self.X = (X.copy() - self._Xmean) / self._Xstd
         else:
-            self._Xmean = np.zeros((1,self.Q))
-            self._Xstd = np.ones((1,self.Q))
+            self._Xmean = np.zeros((1,self.input_dim))
+            self._Xstd = np.ones((1,self.input_dim))
 
         model.model.__init__(self)
 
diff --git a/GPy/core/sparse_GP.py b/GPy/core/sparse_GP.py
index fa96132f..c4fe6763 100644
--- a/GPy/core/sparse_GP.py
+++ b/GPy/core/sparse_GP.py
@@ -13,15 +13,15 @@ class sparse_GP(GPBase):
     Variational sparse GP model
 
     :param X: inputs
-    :type X: np.ndarray (N x Q)
+    :type X: np.ndarray (N 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 Q) | None
+    :type X_variance: np.ndarray (N x input_dim) | None
     :param Z: inducing inputs (optional, see note)
-    :type Z: np.ndarray (M x Q) | None
+    :type Z: np.ndarray (M x input_dim) | None
     :param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
     :type M: int
     :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
@@ -152,7 +152,7 @@ class sparse_GP(GPBase):
         return A + B + C + D + self.likelihood.Z
 
     def _set_params(self, p):
-        self.Z = p[:self.M * self.Q].reshape(self.M, self.Q)
+        self.Z = p[:self.M * self.input_dim].reshape(self.M, self.input_dim)
         self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam])
         self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:])
         self._compute_kernel_matrices()
@@ -252,9 +252,9 @@ class sparse_GP(GPBase):
         Arguments
         ---------
         :param Xnew: The points at which to make a prediction
-        :type Xnew: np.ndarray, Nnew x self.Q
+        :type Xnew: np.ndarray, Nnew x self.input_dim
         :param X_variance_new: The uncertainty in the prediction points
-        :type X_variance_new: np.ndarray, Nnew x self.Q
+        :type X_variance_new: np.ndarray, Nnew x self.input_dim
         :param which_parts:  specifies which outputs kernel(s) to use in prediction
         :type which_parts: ('all', list of bools)
         :param full_cov: whether to return the folll covariance matrix, or just the diagonal
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 22feb28d..441d05a5 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -14,20 +14,20 @@ default_seed = np.random.seed(123344)
 def BGPLVM(seed=default_seed):
     N = 10
     M = 3
-    Q = 2
+    input_dim = 2
     D = 4
     # generate GPLVM-like data
-    X = np.random.rand(N, Q)
-    k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+    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, D).T
 
-    k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) + GPy.kern.white(Q)
-    # k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
-    # 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)
+    k = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.linear(input_dim, ARD=True) + GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.white(input_dim)
+    # k = GPy.kern.rbf(input_dim) + GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim)
+    # k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+    # k = GPy.kern.rbf(input_dim, ARD = False)  + GPy.kern.white(input_dim, 0.00001)
 
-    m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
+    m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel=k, M=M)
     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, M=15, input_dim=4, sigma=.2, plot=False):
     from GPy.util.datasets import swiss_roll
     from GPy.core.transformations import logexp_clipped
 
@@ -79,10 +79,10 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
         from sklearn.manifold.isomap import Isomap
         iso = Isomap().fit(Y)
         X = iso.embedding_
-        if Q > 2:
-            X = np.hstack((X, np.random.randn(N, Q - 2)))
+        if input_dim > 2:
+            X = np.hstack((X, np.random.randn(N, input_dim - 2)))
     except ImportError:
-        X = np.random.randn(N, Q)
+        X = np.random.randn(N, input_dim)
 
     if plot:
         from mpl_toolkits import mplot3d
@@ -98,14 +98,14 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
 
 
     var = .5
-    S = (var * np.ones_like(X) + np.clip(np.random.randn(N, Q) * var ** 2,
+    S = (var * np.ones_like(X) + np.clip(np.random.randn(N, input_dim) * var ** 2,
                                          - (1 - var),
                                          (1 - var))) + .001
     Z = np.random.permutation(X)[:M]
 
-    kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
+    kernel = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.bias(input_dim, np.exp(-2)) + GPy.kern.white(input_dim, np.exp(-2))
 
-    m = Bayesian_GPLVM(Y, Q, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
+    m = Bayesian_GPLVM(Y, input_dim, X=X, X_variance=S, M=M, Z=Z, kernel=kernel)
     m.data_colors = c
     m.data_t = t
 
@@ -118,18 +118,18 @@ 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, input_dim=5, M=25, max_f_eval=4e3, plot=False, **k):
     np.random.seed(0)
     data = GPy.util.datasets.oil()
     from GPy.core.transformations import logexp_clipped
 
     # create simple GP model
-    kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
+    kernel = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.bias(input_dim, np.exp(-2)) + GPy.kern.white(input_dim, np.exp(-2))
     Y = data['X'][:N]
     Yn = Y - Y.mean(0)
     Yn /= Yn.std(0)
 
-    m = GPy.models.Bayesian_GPLVM(Yn, Q, kernel=kernel, M=M, **k)
+    m = GPy.models.Bayesian_GPLVM(Yn, input_dim, kernel=kernel, M=M, **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, M, input_dim, 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
+    M, [_, input_dim] = 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 = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
+    k = kern.linear(input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
+    m = Bayesian_GPLVM(Y, input_dim, init="PCA", M=M, 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, M, input_dim = 15, 8, 8, 100, 3, 5
+    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, input_dim, plot)
 
     from GPy.models import mrd
     from GPy import kern
@@ -258,8 +258,8 @@ 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 = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
+    k = kern.linear(input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2)) # + kern.bias(input_dim)
+    m = Bayesian_GPLVM(Y, input_dim, init="PCA", M=M, kernel=k, _debug=True)
     # m.constrain('variance|noise', logexp_clipped())
     m.ensure_default_constraints()
     m['noise'] = Y.var() / 100.
@@ -279,16 +279,16 @@ 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, M, input_dim = 150, 200, 400, 500, 3, 7
+    slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, input_dim, plot_sim)
 
     from GPy.models import mrd
     from GPy import kern
 
     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, Q=Q, M=M, kernels=k, initx="", initz='permute', **kw)
+    k = kern.linear(input_dim, [.05] * input_dim, ARD=True) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
+    m = mrd.MRD(Ylist, input_dim=input_dim, M=M, kernels=k, initx="", initz='permute', **kw)
 
     for i, Y in enumerate(Ylist):
         m['{}_noise'.format(i + 1)] = Y.var() / 100.
@@ -309,14 +309,14 @@ def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
 def brendan_faces():
     from GPy import kern
     data = GPy.util.datasets.brendan_faces()
-    Q = 2
+    input_dim = 2
     Y = data['Y'][0:-1:10, :]
     # Y = data['Y']
     Yn = Y - Y.mean()
     Yn /= Yn.std()
 
-    m = GPy.models.GPLVM(Yn, Q)
-    # m = GPy.models.Bayesian_GPLVM(Yn, Q, M=100)
+    m = GPy.models.GPLVM(Yn, input_dim)
+    # m = GPy.models.Bayesian_GPLVM(Yn, input_dim, M=100)
 
     # optimize
     m.constrain('rbf|noise|white', GPy.core.transformations.logexp_clipped())
@@ -379,11 +379,11 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True):
 #     X -= X.mean(axis=0)
 #     X /= X.std(axis=0)
 #
-#     Q = 10
+#     input_dim = 10
 #     M = 30
 #
-#     kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q) + GPy.kern.white(Q)
-#     m = GPy.models.Bayesian_GPLVM(X, Q, kernel=kernel, M=M)
+#     kernel = GPy.kern.rbf(input_dim, ARD=True) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
+#     m = GPy.models.Bayesian_GPLVM(X, input_dim, kernel=kernel, M=M)
 #     # m.scale_factor = 100.0
 #     m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)')
 #     from sklearn import cluster
diff --git a/GPy/inference/SGD.py b/GPy/inference/SGD.py
index 7543a819..c2a77e40 100644
--- a/GPy/inference/SGD.py
+++ b/GPy/inference/SGD.py
@@ -95,11 +95,11 @@ class opt_SGD(Optimizer):
         i = 0
 
         for s in param_shapes:
-            N, Q = s
-            X = x[i:i+N*Q].reshape(N, Q)
+            N, input_dim = s
+            X = x[i:i+N*input_dim].reshape(N, input_dim)
             X = X[samples]
             subset = np.append(subset, X.flatten())
-            i += N*Q
+            i += N*input_dim
 
         subset = np.append(subset, x[i:])
 
@@ -167,17 +167,17 @@ class opt_SGD(Optimizer):
         # self.model.constrained_positive_indices = p
         self.model.constrained_indices = c
 
-    def get_param_shapes(self, N = None, Q = None):
+    def get_param_shapes(self, N = None, input_dim = None):
         model_name = self.model.__class__.__name__
         if model_name == 'GPLVM':
-            return [(N, Q)]
+            return [(N, input_dim)]
         if model_name == 'Bayesian_GPLVM':
-            return [(N, Q), (N, Q)]
+            return [(N, input_dim), (N, input_dim)]
         else:
             raise NotImplementedError
 
     def step_with_missing_data(self, f_fp, X, step, shapes):
-        N, Q = X.shape
+        N, input_dim = X.shape
 
         if not sp.sparse.issparse(self.model.likelihood.Y):
             Y = self.model.likelihood.Y
@@ -269,7 +269,7 @@ class opt_SGD(Optimizer):
         self.model.likelihood._offset = 0.0
         self.model.likelihood._scale = 1.0
 
-        N, Q = self.model.X.shape
+        N, input_dim = self.model.X.shape
         D = self.model.likelihood.Y.shape[1]
         num_params = self.model._get_params()
         self.trace = []
@@ -302,7 +302,7 @@ class opt_SGD(Optimizer):
                 self.model.likelihood._set_params(sigma)
 
                 if missing_data:
-                    shapes = self.get_param_shapes(N, Q)
+                    shapes = self.get_param_shapes(N, input_dim)
                     f, step, Nj = self.step_with_missing_data(f_fp, X, step, shapes)
                 else:
                     self.model.likelihood.YYT = np.dot(self.model.likelihood.Y, self.model.likelihood.Y.T)
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 28e33b4b..0d3b4ccf 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -300,15 +300,15 @@ class kern(parameterised):
         return target
 
     def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
-        """return shapes are N,M,Q"""
+        """return shapes are N,M,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 Q)
-        :param mu, S: np.ndarrays of means and variances (each N x Q)
+        :param Z: np.ndarray of inducing inputs (M x input_dim)
+        :param mu, S: np.ndarrays of means and variances (each N x input_dim)
         :returns psi2: np.ndarray (N,M,M)
         """
         target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py
index af3e60ea..8744eda0 100644
--- a/GPy/kern/linear.py
+++ b/GPy/kern/linear.py
@@ -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,Q """
+        """Think N,M,M,input_dim """
         self._psi_computations(Z, mu, S)
         AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]
         AZZA = AZZA + AZZA.swapaxes(1, 2)
@@ -192,11 +192,11 @@ class linear(kernpart):
               else
                 factor = 2.0*dL_dpsi2(n,m,mm);
 
-              for(q=0;q<Q;q++){
+              for(q=0;q<input_dim;q++){
 
                 //take the dot product of mu[n,:] and AZZA[:,m,mm,q] TODO: blas!
                 tmp = 0.0;
-                for(qq=0;qq<Q;qq++){
+                for(qq=0;qq<input_dim;qq++){
                   tmp += mu(n,qq)*AZZA(qq,m,mm,q);
                 }
 
@@ -215,9 +215,9 @@ class linear(kernpart):
                          'extra_compile_args': ['-fopenmp -O3'],  #-march=native'],
                          'extra_link_args'   : ['-lgomp']}
 
-        N,M,Q = mu.shape[0],Z.shape[0],mu.shape[1]
+        N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','Q','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
+                     arg_names=['N','M','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
                      type_converters=weave.converters.blitz,**weave_options)
 
 
@@ -232,7 +232,7 @@ class linear(kernpart):
         int n,m,mm,q;
         #pragma omp parallel for private(n,mm,q)
         for(m=0;m<M;m++){
-          for(q=0;q<Q;q++){
+          for(q=0;q<input_dim;q++){
             for(mm=0;mm<M;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,Q = mu.shape[0],Z.shape[0],mu.shape[1]
+        N,M,input_dim = mu.shape[0],Z.shape[0],mu.shape[1]
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','Q','AZA','target','dL_dpsi2'],
+                     arg_names=['N','M','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,Q
+            # self.ZZ = Z[:,None,:]*Z[None,:,:] # M,M,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 Q]!
+            self.ZAinner = np.dot(self.ZA, self.inner).swapaxes(0, 1)  # NOTE: self.ZAinner \in [M x N x input_dim]!
             self._psi2 = np.dot(self.ZAinner, self.ZA.T)
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index 8b644241..a89b7d45 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -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, Q
-        term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, Q
+        term1 = self._psi2_Zdist/self.lengthscale2 # M, M, input_dim
+        term2 = self._psi2_mudist/self._psi2_denom/self.lengthscale2 # N, M, M, 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,Q """
+        """Think N,M,M,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,Q
-            self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,Q
-            self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,Q
+            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._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,9 +258,9 @@ 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,Q
+            self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,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,Q
+            #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,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
@@ -269,11 +269,11 @@ class rbf(kernpart):
             self._Z, self._mu, self._S = Z, mu,S
 
     def weave_psi2(self,mu,Zhat):
-        N,Q = mu.shape
+        N,input_dim = mu.shape
         M = Zhat.shape[0]
 
-        mudist = np.empty((N,M,M,Q))
-        mudist_sq = np.empty((N,M,M,Q))
+        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))
 
@@ -284,7 +284,7 @@ class rbf(kernpart):
         if self.ARD:
             lengthscale2 = self.lengthscale2
         else:
-            lengthscale2 = np.ones(Q)*self.lengthscale2
+            lengthscale2 = np.ones(input_dim)*self.lengthscale2
         code = """
         double tmp;
 
@@ -292,7 +292,7 @@ class rbf(kernpart):
         for (int n=0; n<N; n++){
             for (int m=0; m<M; m++){
                for (int mm=0; mm<(m+1); mm++){
-                   for (int q=0; q<Q; q++){
+                   for (int q=0; q<input_dim; q++){
                        //compute mudist
                        tmp = mu(n,q) - Zhat(m,mm,q);
                        mudist(n,m,mm,q) = tmp;
@@ -325,7 +325,7 @@ class rbf(kernpart):
         #include <math.h>
         """
         weave.inline(code, support_code=support_code, libraries=['gomp'],
-                     arg_names=['N','M','Q','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'],
+                     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'],
                      type_converters=weave.converters.blitz,**self.weave_options)
 
         return mudist,mudist_sq, psi2_exponent, psi2
diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py
index b7f7b42b..ce8b89a3 100644
--- a/GPy/models/Bayesian_GPLVM.py
+++ b/GPy/models/Bayesian_GPLVM.py
@@ -22,13 +22,13 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
     :param Y: observed data (np.ndarray) or GPy.likelihood
     :type Y: np.ndarray| GPy.likelihood instance
-    :param Q: latent dimensionality
-    :type Q: int
+    :param input_dim: latent dimensionality
+    :type input_dim: int
     :param init: initialisation method for the latent space
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, likelihood_or_Y, Q, X=None, X_variance=None, init='PCA', M=10,
+    def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', M=10,
                  Z=None, kernel=None, oldpsave=10, _debug=False,
                  **kwargs):
         if type(likelihood_or_Y) is np.ndarray:
@@ -37,7 +37,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
             likelihood = likelihood_or_Y
 
         if X == None:
-            X = self.initialise_latent(init, Q, likelihood.Y)
+            X = self.initialise_latent(init, input_dim, likelihood.Y)
         self.init = init
 
         if X_variance is None:
@@ -48,7 +48,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         assert Z.shape[1] == X.shape[1]
 
         if kernel is None:
-            kernel = kern.rbf(Q) + kern.white(Q)
+            kernel = kern.rbf(input_dim) + kern.white(input_dim)
 
         self.oldpsave = oldpsave
         self._oldps = []
@@ -78,8 +78,8 @@ class Bayesian_GPLVM(sparse_GP, 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.Q)] for n in range(self.N)], [])
-        S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] 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.N)], [])
+        S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
         return (X_names + S_names + sparse_GP._get_param_names(self))
 
     def _get_params(self):
@@ -101,10 +101,10 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
     def _set_params(self, x, save_old=True, save_count=0):
 #         try:
             x = self._clipped(x)
-            N, Q = self.N, self.Q
-            self.X = x[:self.X.size].reshape(N, Q).copy()
-            self.X_variance = x[(N * Q):(2 * N * Q)].reshape(N, Q).copy()
-            sparse_GP._set_params(self, x[(2 * N * Q):])
+            N, input_dim = self.N, 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()
+            sparse_GP._set_params(self, x[(2 * N * input_dim):])
 #             self.oldps = x
 #         except (LinAlgError, FloatingPointError, ZeroDivisionError):
 #             print "\rWARNING: Caught LinAlgError, continueing without setting            "
@@ -131,7 +131,7 @@ class Bayesian_GPLVM(sparse_GP, 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.Q * self.N
+        return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.N
 
     def log_likelihood(self):
         ll = sparse_GP.log_likelihood(self)
@@ -196,16 +196,16 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         """
         assert not self.likelihood.is_heteroscedastic
         N_test = Y.shape[0]
-        Q = self.Z.shape[1]
-        means = np.zeros((N_test, Q))
-        covars = np.zeros((N_test, Q))
+        input_dim = self.Z.shape[1]
+        means = np.zeros((N_test, input_dim))
+        covars = np.zeros((N_test, input_dim))
 
         dpsi0 = -0.5 * self.D * self.likelihood.precision
         dpsi2 = self.dL_dpsi2[0][None, :, :] # TODO: this may change if we ignore het. likelihoods
         V = self.likelihood.precision * Y
         dpsi1 = np.dot(self.Cpsi1V, V.T)
 
-        start = np.zeros(self.Q * 2)
+        start = np.zeros(self.input_dim * 2)
 
         for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
             args = (self.kern, self.Z, dpsi0, dpsi1_n, dpsi2)
@@ -222,12 +222,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         """
         Plot latent space X in 1D:
 
-            -if fig is given, create Q subplots in fig and plot in these
-            -if ax is given plot Q 1D latent space plots of X into each `axis`
+            -if fig is given, create input_dim subplots in fig and plot in these
+            -if ax is given plot input_dim 1D latent space plots of X into each `axis`
             -if neither fig nor ax is given create a figure with fignum and plot in there
 
         colors:
-            colors of different latent space dimensions Q
+            colors of different latent space dimensions input_dim
         """
         import pylab
         if ax is None:
@@ -245,7 +245,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
             elif isinstance(ax, (tuple, list)):
                 ax = ax[i]
             else:
-                raise ValueError("Need one ax per latent dimnesion Q")
+                raise ValueError("Need one ax per latent dimnesion input_dim")
             ax.plot(self.X, c='k', alpha=.3)
             plots.extend(ax.plot(x, self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
             ax.fill_between(x,
@@ -262,7 +262,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
         return fig
 
     def __getstate__(self):
-        return (self.likelihood, self.Q, self.X, self.X_variance,
+        return (self.likelihood, self.input_dim, self.X, self.X_variance,
                 self.init, self.M, self.Z, self.kern,
                 self.oldpsave, self._debug)
 
@@ -271,12 +271,12 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 
     def _debug_filter_params(self, x):
         start, end = 0, self.X.size,
-        X = x[start:end].reshape(self.N, self.Q)
+        X = x[start:end].reshape(self.N, self.input_dim)
         start, end = end, end + self.X_variance.size
-        X_v = x[start:end].reshape(self.N, self.Q)
-        start, end = end, end + (self.M * self.Q)
-        Z = x[start:end].reshape(self.M, self.Q)
-        start, end = end, end + self.Q
+        X_v = x[start:end].reshape(self.N, self.input_dim)
+        start, end = end, end + (self.M * self.input_dim)
+        Z = x[start:end].reshape(self.M, self.input_dim)
+        start, end = end, end + self.input_dim
         theta = x[start:]
         return X, X_v, Z, theta
 
@@ -414,24 +414,24 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
 #         caxkmmdl = divider.append_axes("right", "5%", pad="1%")
 #         cbarkmmdl = pylab.colorbar(imkmmdl, cax=caxkmmdl)
 
-#         Qleg = ax1.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
-#                    loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.15, 1, 1.15),
+#         input_dimleg = ax1.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
+#                    loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.15, 1, 1.15),
 #                    borderaxespad=0, mode="expand")
-        ax2.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
-                   loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
+        ax2.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
+                   loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
                    borderaxespad=0, mode="expand")
-        ax3.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
-                   loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
+        ax3.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
+                   loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
                    borderaxespad=0, mode="expand")
-        ax4.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
-                   loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
+        ax4.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
+                   loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
                    borderaxespad=0, mode="expand")
-        ax5.legend(Xlatentplts, [r"$Q_{}$".format(i + 1) for i in range(self.Q)],
-                   loc=3, ncol=self.Q, bbox_to_anchor=(0, 1.1, 1, 1.1),
+        ax5.legend(Xlatentplts, [r"$input_dim_{}$".format(i + 1) for i in range(self.input_dim)],
+                   loc=3, ncol=self.input_dim, bbox_to_anchor=(0, 1.1, 1, 1.1),
                    borderaxespad=0, mode="expand")
         Lleg = ax1.legend()
         Lleg.draggable()
-#         ax1.add_artist(Qleg)
+#         ax1.add_artist(input_dimleg)
 
         indicatorKL, = ax1.plot(kllls[self.showing, 0], kllls[self.showing, 2], 'o', c=KL.get_color())
         indicatorLL, = ax1.plot(kllls[self.showing, 0], kllls[self.showing, 1] - kllls[self.showing, 2], 'o', c=LL.get_color())
diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py
index bfe422c8..5589304a 100644
--- a/GPy/models/GPLVM.py
+++ b/GPy/models/GPLVM.py
@@ -20,35 +20,35 @@ class GPLVM(GP):
 
     :param Y: observed data
     :type Y: np.ndarray
-    :param Q: latent dimensionality
-    :type Q: int
+    :param input_dim: latent dimensionality
+    :type input_dim: int
     :param init: initialisation method for the latent space
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, Q, init='PCA', X = None, kernel=None, normalize_Y=False):
+    def __init__(self, Y, input_dim, init='PCA', X = None, kernel=None, normalize_Y=False):
         if X is None:
-            X = self.initialise_latent(init, Q, Y)
+            X = self.initialise_latent(init, input_dim, Y)
         if kernel is None:
-            kernel = kern.rbf(Q, ARD=Q>1) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
+            kernel = kern.rbf(input_dim, ARD=input_dim>1) + kern.bias(input_dim, np.exp(-2)) + kern.white(input_dim, np.exp(-2))
         likelihood = Gaussian(Y, normalize=normalize_Y)
         GP.__init__(self, X, likelihood, kernel, normalize_X=False)
         self._set_params(self._get_params())
 
-    def initialise_latent(self, init, Q, Y):
+    def initialise_latent(self, init, input_dim, Y):
         if init == 'PCA':
-            return PCA(Y, Q)[0]
+            return PCA(Y, input_dim)[0]
         else:
-            return np.random.randn(Y.shape[0], Q)
+            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.Q)] 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.N)],[]) + 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.Q].reshape(self.N,self.Q).copy()
+        self.X = x[:self.N*self.input_dim].reshape(self.N,self.input_dim).copy()
         GP._set_params(self, x[self.X.size:])
 
     def _log_likelihood_gradients(self):
diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py
index 145b4b11..6e44baf6 100644
--- a/GPy/models/generalized_FITC.py
+++ b/GPy/models/generalized_FITC.py
@@ -20,15 +20,15 @@ class generalized_FITC(sparse_GP):
     Naish-Guzman, A. and Holden, S. (2008) implemantation of EP with FITC.
 
     :param X: inputs
-    :type X: np.ndarray (N x Q)
+    :type X: np.ndarray (N x input_dim)
     :param likelihood: a likelihood instance, containing the observed data
     :type likelihood: GPy.likelihood.(Gaussian | EP)
     :param kernel : the kernel/covariance function. See link kernels
     :type kernel: a GPy kernel
     :param X_variance: The variance in the measurements of X (Gaussian variance)
-    :type X_variance: np.ndarray (N x Q) | None
+    :type X_variance: np.ndarray (N x input_dim) | None
     :param Z: inducing inputs (optional, see note)
-    :type Z: np.ndarray (M x Q) | None
+    :type Z: np.ndarray (M x input_dim) | None
     :param M : Number of inducing points (optional, default 10. Ignored if Z is not None)
     :type M: int
     :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales)
@@ -44,7 +44,7 @@ class generalized_FITC(sparse_GP):
         self._set_params(self._get_params())
 
     def _set_params(self, p):
-        self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
+        self.Z = p[:self.M*self.input_dim].reshape(self.M, self.input_dim)
         self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
         self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:])
         self._compute_kernel_matrices()
diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py
index e1bfa947..641438b0 100644
--- a/GPy/models/mrd.py
+++ b/GPy/models/mrd.py
@@ -23,8 +23,8 @@ class MRD(model):
     :type likelihood_list: [GPy.likelihood] | [Y1..Yy]
     :param names: names for different gplvm models
     :type names: [str]
-    :param Q: latent dimensionality (will raise
-    :type Q: int
+    :param input_dim: latent dimensionality (will raise
+    :type input_dim: int
     :param initx: initialisation method for the latent space
     :type initx: 'PCA'|'random'
     :param initz: initialisation method for inducing inputs
@@ -45,7 +45,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, Q, M=10, names=None,
+    def __init__(self, likelihood_or_Y_list, input_dim, M=10, names=None,
                  kernels=None, initx='PCA',
                  initz='permute', _debug=False, **kw):
         if names is None:
@@ -61,14 +61,14 @@ class MRD(model):
             assert all([isinstance(k, kern) for k in kernels]), "invalid kernel object detected!"
         assert not ('kernel' in kw), "pass kernels through `kernels` argument"
 
-        self.Q = Q
+        self.input_dim = input_dim
         self.M = M
         self._debug = _debug
 
         self._init = True
         X = self._init_X(initx, likelihood_or_Y_list)
         Z = self._init_Z(initz, X)
-        self.bgplvms = [Bayesian_GPLVM(l, Q=Q, kernel=k, X=X, Z=Z, M=self.M, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
+        self.bgplvms = [Bayesian_GPLVM(l, input_dim=input_dim, kernel=k, X=X, Z=Z, M=self.M, **kw) for l, k in zip(likelihood_or_Y_list, kernels)]
         del self._init
 
         self.gref = self.bgplvms[0]
@@ -76,8 +76,8 @@ class MRD(model):
         self.nparams = nparams.cumsum()
 
         self.N = self.gref.N
-        self.NQ = self.N * self.Q
-        self.MQ = self.M * self.Q
+        self.NQ = self.N * self.input_dim
+        self.MQ = self.M * self.input_dim
 
         model.__init__(self) # @UndefinedVariable
         self._set_params(self._get_params())
@@ -143,8 +143,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.Q)] for n in range(self.N)], [])
-        # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.Q)] 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.N)], [])
+        # S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.N)], [])
         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),
@@ -170,9 +170,9 @@ class MRD(model):
         return params
 
 #     def _set_var_params(self, g, X, X_var, Z):
-#         g.X = X.reshape(self.N, self.Q)
-#         g.X_variance = X_var.reshape(self.N, self.Q)
-#         g.Z = Z.reshape(self.M, self.Q)
+#         g.X = X.reshape(self.N, self.input_dim)
+#         g.X_variance = X_var.reshape(self.N, self.input_dim)
+#         g.Z = Z.reshape(self.M, self.input_dim)
 #
 #     def _set_kern_params(self, g, p):
 #         g.kern._set_params(p[:g.kern.Nparam])
@@ -235,13 +235,13 @@ class MRD(model):
                 Ylist.append(likelihood_or_Y.Y)
         del likelihood_list
         if init in "PCA_concat":
-            X = PCA(numpy.hstack(Ylist), self.Q)[0]
+            X = PCA(numpy.hstack(Ylist), self.input_dim)[0]
         elif init in "PCA_single":
-            X = numpy.zeros((Ylist[0].shape[0], self.Q))
-            for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.Q), len(Ylist)), Ylist):
+            X = numpy.zeros((Ylist[0].shape[0], self.input_dim))
+            for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.input_dim), len(Ylist)), Ylist):
                 X[:, qs] = PCA(Y, len(qs))[0]
         else: # init == 'random':
-            X = numpy.random.randn(Ylist[0].shape[0], self.Q)
+            X = numpy.random.randn(Ylist[0].shape[0], self.input_dim)
         self.X = X
         return X
 
@@ -252,7 +252,7 @@ class MRD(model):
         if init in "permute":
             Z = numpy.random.permutation(X.copy())[:self.M]
         elif init in "random":
-            Z = numpy.random.randn(self.M, self.Q) * X.var()
+            Z = numpy.random.randn(self.M, self.input_dim) * X.var()
         self.Z = Z
         return Z
 
@@ -265,7 +265,7 @@ class MRD(model):
             elif isinstance(axes, (tuple, list)):
                 axes = axes[i]
             else:
-                raise ValueError("Need one axes per latent dimension Q")
+                raise ValueError("Need one axes per latent dimension input_dim")
             plotf(i, g, axes)
         pylab.draw()
         if axes is None:
diff --git a/GPy/models/sparse_GPLVM.py b/GPy/models/sparse_GPLVM.py
index 8760ed2e..ce71cd9b 100644
--- a/GPy/models/sparse_GPLVM.py
+++ b/GPy/models/sparse_GPLVM.py
@@ -17,25 +17,25 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
 
     :param Y: observed data
     :type Y: np.ndarray
-    :param Q: latent dimensionality
-    :type Q: int
+    :param input_dim: latent dimensionality
+    :type input_dim: int
     :param init: initialisation method for the latent space
     :type init: 'PCA'|'random'
 
     """
-    def __init__(self, Y, Q, kernel=None, init='PCA', M=10):
-        X = self.initialise_latent(init, Q, Y)
+    def __init__(self, Y, input_dim, kernel=None, init='PCA', M=10):
+        X = self.initialise_latent(init, input_dim, Y)
         sparse_GP_regression.__init__(self, X, Y, kernel=kernel,M=M)
 
     def _get_param_names(self):
-        return (sum([['X_%i_%i'%(n,q) for q in range(self.Q)] 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.N)],[])
                 + sparse_GP_regression._get_param_names(self))
 
     def _get_params(self):
         return np.hstack((self.X.flatten(), sparse_GP_regression._get_params(self)))
 
     def _set_params(self,x):
-        self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
+        self.X = x[:self.X.size].reshape(self.N,self.input_dim).copy()
         sparse_GP_regression._set_params(self, x[self.X.size:])
 
     def log_likelihood(self):
diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py
index 5396e175..ae72983a 100644
--- a/GPy/testing/bgplvm_tests.py
+++ b/GPy/testing/bgplvm_tests.py
@@ -7,67 +7,67 @@ import GPy
 
 class BGPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
         Y -= Y.mean(axis=0)
-        k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k,  M=M)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_linear_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
         Y -= Y.mean(axis=0)
-        k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k,  M=M)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
         Y -= Y.mean(axis=0)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k,  M=M)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) +  GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
         Y -= Y.mean(axis=0)
-        k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        k = GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k,  M=M)
         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, Q, D = 30, 5, 4, 30
-        X = np.random.rand(N, Q)
-        k = GPy.kern.linear(Q) +  GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
         Y -= Y.mean(axis=0)
-        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.Bayesian_GPLVM(Y, Q, kernel = k,  M=M)
+        k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel = k,  M=M)
         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 51828768..b1721a3c 100644
--- a/GPy/testing/gplvm_tests.py
+++ b/GPy/testing/gplvm_tests.py
@@ -7,37 +7,37 @@ import GPy
 
 class GPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.GPLVM(Y, Q, kernel = k)
+        k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.GPLVM(Y, input_dim, kernel = k)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_linear_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.GPLVM(Y, Q, kernel = k)
+        k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.GPLVM(Y, input_dim, kernel = k)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.GPLVM(Y, Q, kernel = k)
+        k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.GPLVM(Y, input_dim, kernel = k)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
diff --git a/GPy/testing/mrd_tests.py b/GPy/testing/mrd_tests.py
index 903cddfc..25adbca6 100644
--- a/GPy/testing/mrd_tests.py
+++ b/GPy/testing/mrd_tests.py
@@ -14,16 +14,16 @@ class MRDTests(unittest.TestCase):
 
     def test_gradients(self):
         num_m = 3
-        N, M, Q, D = 20, 8, 6, 20
-        X = np.random.rand(N, Q)
+        N, M, input_dim, D = 20, 8, 6, 20
+        X = np.random.rand(N, input_dim)
 
-        k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q)
+        k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
         K = k.K(X)
 
         Ylist = [np.random.multivariate_normal(np.zeros(N), K, D).T for _ in range(num_m)]
         likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
 
-        m = GPy.models.MRD(likelihood_list, Q=Q, kernels=k, M=M)
+        m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, M=M)
         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 7c41098f..23f841b5 100644
--- a/GPy/testing/psi_stat_gradient_tests.py
+++ b/GPy/testing/psi_stat_gradient_tests.py
@@ -16,25 +16,25 @@ class PsiStatModel(model):
         self.X = X
         self.X_variance = X_variance
         self.Z = Z
-        self.N, self.Q = X.shape
-        self.M, Q = Z.shape
-        assert self.Q == Q, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
+        self.N, self.input_dim = X.shape
+        self.M, input_dim = Z.shape
+        assert self.input_dim == input_dim, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
         self.kern = kernel
         super(PsiStatModel, self).__init__()
         self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance)
     def _get_param_names(self):
-        Xnames = ["{}_{}_{}".format(what, i, j) for what, i, j in itertools.product(['X', 'X_variance'], range(self.N), range(self.Q))]
-        Znames = ["Z_{}_{}".format(i, j) for i, j in itertools.product(range(self.M), range(self.Q))]
+        Xnames = ["{}_{}_{}".format(what, i, j) for what, i, j in itertools.product(['X', 'X_variance'], range(self.N), range(self.input_dim))]
+        Znames = ["Z_{}_{}".format(i, j) for i, j in itertools.product(range(self.M), range(self.input_dim))]
         return Xnames + Znames + self.kern._get_param_names()
     def _get_params(self):
         return numpy.hstack([self.X.flatten(), self.X_variance.flatten(), self.Z.flatten(), self.kern._get_params()])
     def _set_params(self, x, save_old=True, save_count=0):
         start, end = 0, self.X.size
-        self.X = x[start:end].reshape(self.N, self.Q)
+        self.X = x[start:end].reshape(self.N, self.input_dim)
         start, end = end, end + self.X_variance.size
-        self.X_variance = x[start: end].reshape(self.N, self.Q)
+        self.X_variance = x[start: end].reshape(self.N, self.input_dim)
         start, end = end, end + self.Z.size
-        self.Z = x[start: end].reshape(self.M, self.Q)
+        self.Z = x[start: end].reshape(self.M, self.input_dim)
         self.kern._set_params(x[end:])
     def log_likelihood(self):
         return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum()
@@ -43,24 +43,24 @@ class PsiStatModel(model):
         try:
             psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance)
         except AttributeError:
-            psiZ = numpy.zeros(self.M * self.Q)
+            psiZ = numpy.zeros(self.M * self.input_dim)
         thetagrad = self.kern.__getattribute__("d" + self.which + "_dtheta")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance).flatten()
         return numpy.hstack((psimu.flatten(), psiS.flatten(), psiZ.flatten(), thetagrad))
 
 class DPsiStatTest(unittest.TestCase):
-    Q = 5
+    input_dim = 5
     N = 50
     M = 10
     D = 20
-    X = numpy.random.randn(N, Q)
+    X = numpy.random.randn(N, input_dim)
     X_var = .5 * numpy.ones_like(X) + .4 * numpy.clip(numpy.random.randn(*X.shape), 0, 1)
     Z = numpy.random.permutation(X)[:M]
-    Y = X.dot(numpy.random.randn(Q, D))
-#     kernels = [GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)), GPy.kern.rbf(Q, ARD=True), GPy.kern.bias(Q)]
+    Y = X.dot(numpy.random.randn(input_dim, D))
+#     kernels = [GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)), GPy.kern.rbf(input_dim, ARD=True), GPy.kern.bias(input_dim)]
 
-    kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
-               GPy.kern.linear(Q) + GPy.kern.bias(Q),
-               GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
+    kernels = [GPy.kern.linear(input_dim), GPy.kern.rbf(input_dim), GPy.kern.bias(input_dim),
+               GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim),
+               GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim)]
 
     def testPsi0(self):
         for k in self.kernels:
@@ -108,31 +108,31 @@ if __name__ == "__main__":
     import sys
     interactive = 'i' in sys.argv
     if interactive:
-#         N, M, Q, D = 30, 5, 4, 30
-#         X = numpy.random.rand(N, Q)
-#         k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
+#         N, M, input_dim, D = 30, 5, 4, 30
+#         X = numpy.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 = numpy.random.multivariate_normal(numpy.zeros(N), K, D).T
 #         Y -= Y.mean(axis=0)
-#         k = GPy.kern.linear(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-#         m = GPy.models.Bayesian_GPLVM(Y, Q, kernel=k, M=M)
+#         k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+#         m = GPy.models.Bayesian_GPLVM(Y, input_dim, kernel=k, M=M)
 #         m.ensure_default_constraints()
 #         m.randomize()
 # #         self.assertTrue(m.checkgrad())
         numpy.random.seed(0)
-        Q = 5
+        input_dim = 5
         N = 50
         M = 10
         D = 15
-        X = numpy.random.randn(N, Q)
+        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]
-        Y = X.dot(numpy.random.randn(Q, D))
-#         kernel = GPy.kern.bias(Q)
+        Y = X.dot(numpy.random.randn(input_dim, D))
+#         kernel = GPy.kern.bias(input_dim)
 #
-#         kernels = [GPy.kern.linear(Q), GPy.kern.rbf(Q), GPy.kern.bias(Q),
-#                GPy.kern.linear(Q) + GPy.kern.bias(Q),
-#                GPy.kern.rbf(Q) + GPy.kern.bias(Q)]
+#         kernels = [GPy.kern.linear(input_dim), GPy.kern.rbf(input_dim), GPy.kern.bias(input_dim),
+#                GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim),
+#                GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim)]
 
 #         for k in kernels:
 #             m = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
@@ -140,18 +140,18 @@ if __name__ == "__main__":
 #             assert m.checkgrad(), "{} x psi1".format("+".join(map(lambda x: x.name, k.parts)))
 #
 #         m0 = PsiStatModel('psi0', X=X, X_variance=X_var, Z=Z,
-#                          M=M, kernel=GPy.kern.linear(Q))
+#                          M=M, kernel=GPy.kern.linear(input_dim))
 #         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
 #                          M=M, kernel=kernel)
 #         m1 = PsiStatModel('psi1', X=X, X_variance=X_var, Z=Z,
 #                          M=M, kernel=kernel)
 #         m2 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
-#                          M=M, kernel=GPy.kern.rbf(Q))
+#                          M=M, kernel=GPy.kern.rbf(input_dim))
         m3 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
-                         M=M, kernel=GPy.kern.linear(Q, ARD=True, variances=numpy.random.rand(Q)))
+                         M=M, kernel=GPy.kern.linear(input_dim, ARD=True, variances=numpy.random.rand(input_dim)))
         m3.ensure_default_constraints()
-        # + GPy.kern.bias(Q))
+        # + GPy.kern.bias(input_dim))
 #         m4 = PsiStatModel('psi2', X=X, X_variance=X_var, Z=Z,
-#                          M=M, kernel=GPy.kern.rbf(Q) + GPy.kern.bias(Q))
+#                          M=M, kernel=GPy.kern.rbf(input_dim) + GPy.kern.bias(input_dim))
     else:
         unittest.main()
diff --git a/GPy/testing/sparse_gplvm_tests.py b/GPy/testing/sparse_gplvm_tests.py
index 35fa4fcf..a790ce54 100644
--- a/GPy/testing/sparse_gplvm_tests.py
+++ b/GPy/testing/sparse_gplvm_tests.py
@@ -7,38 +7,38 @@ import GPy
 
 class sparse_GPLVMTests(unittest.TestCase):
     def test_bias_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
+        k = GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
         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, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.linear(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
+        k = GPy.kern.linear(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
 
     def test_rbf_kern(self):
-        N, M, Q, D = 10, 3, 2, 4
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
+        N, M, 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,D).T
-        k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
-        m = GPy.models.sparse_GPLVM(Y, Q, kernel = k, M=M)
+        k = GPy.kern.rbf(input_dim) + GPy.kern.white(input_dim, 0.00001)
+        m = GPy.models.sparse_GPLVM(Y, input_dim, kernel = k, M=M)
         m.ensure_default_constraints()
         m.randomize()
         self.assertTrue(m.checkgrad())
diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py
index b2c8196b..8224c7a8 100644
--- a/GPy/testing/unit_tests.py
+++ b/GPy/testing/unit_tests.py
@@ -144,23 +144,23 @@ class GradientTests(unittest.TestCase):
 
     def test_GPLVM_rbf_bias_white_kern_2D(self):
         """ Testing GPLVM with rbf + bias and white kernel """
-        N, Q, D = 50, 1, 2
-        X = np.random.rand(N, Q)
-        k = GPy.kern.rbf(Q, 0.5, 0.9*np.ones((1,))) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
+        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)
         Y = np.random.multivariate_normal(np.zeros(N),K,D).T
-        m = GPy.models.GPLVM(Y, Q, kernel = k)
+        m = GPy.models.GPLVM(Y, input_dim, kernel = k)
         m.ensure_default_constraints()
         self.assertTrue(m.checkgrad())
 
     def test_GPLVM_rbf_linear_white_kern_2D(self):
         """ Testing GPLVM with rbf + bias and white kernel """
-        N, Q, D = 50, 1, 2
-        X = np.random.rand(N, Q)
-        k = GPy.kern.linear(Q) + GPy.kern.bias(Q, 0.1) + GPy.kern.white(Q, 0.05)
+        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)
         Y = np.random.multivariate_normal(np.zeros(N),K,D).T
-        m = GPy.models.GPLVM(Y, Q, init = 'PCA', kernel = k)
+        m = GPy.models.GPLVM(Y, input_dim, init = 'PCA', kernel = k)
         m.ensure_default_constraints()
         self.assertTrue(m.checkgrad())
 
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index da23a0a8..6e7d26d8 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -162,19 +162,19 @@ def multiple_pdinv(A):
     return np.dstack(invs),np.array(halflogdets)
 
 
-def PCA(Y, Q):
+def PCA(Y, input_dim):
     """
     Principal component analysis: maximum likelihood solution by SVD
 
     Arguments
     ---------
     :param Y: NxD np.array of data
-    :param Q: int, dimension of projection
+    :param input_dim: int, dimension of projection
 
     Returns
     -------
-    :rval X: - NxQ np.array of dimensionality reduced data
-    W - QxD mapping from X to Y
+    :rval X: - Nxinput_dim np.array of dimensionality reduced data
+    W - input_dimxD mapping from X to Y
     """
     if not np.allclose(Y.mean(axis=0), 0.0):
         print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)"
@@ -182,7 +182,7 @@ def PCA(Y, Q):
         #Y -= Y.mean(axis=0)
 
     Z = linalg.svd(Y-Y.mean(axis=0), full_matrices = False)
-    [X, W] = [Z[0][:,0:Q], np.dot(np.diag(Z[1]), Z[2]).T[:,0:Q]]
+    [X, W] = [Z[0][:,0:input_dim], np.dot(np.diag(Z[1]), Z[2]).T[:,0:input_dim]]
     v = X.std(axis=0)
     X /= v;
     W *= v;
diff --git a/GPy/util/plot_latent.py b/GPy/util/plot_latent.py
index 47896e48..c09b5c93 100644
--- a/GPy/util/plot_latent.py
+++ b/GPy/util/plot_latent.py
@@ -14,10 +14,10 @@ def plot_latent(model, labels=None, which_indices=None, resolution=50, ax=None,
     if labels is None:
         labels = np.ones(model.N)
     if which_indices is None:
-        if model.Q==1:
+        if model.input_dim==1:
             input_1 = 0
             input_2 = None
-        if model.Q==2:
+        if model.input_dim==2:
             input_1, input_2 = 0,1
         else:
             try:
@@ -55,7 +55,7 @@ def plot_latent(model, labels=None, which_indices=None, resolution=50, ax=None,
             m = marker
 
         index = np.nonzero(labels==ul)[0]
-        if model.Q==1:
+        if model.input_dim==1:
             x = model.X[index,input_1]
             y = np.zeros(index.size)
         else:
diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py
index b3429850..47e632bc 100644
--- a/GPy/util/visualize.py
+++ b/GPy/util/visualize.py
@@ -74,7 +74,7 @@ class lvm(data_show):
         self.called = False
         self.move_on = False
         self.latent_index = latent_index
-        self.latent_dim = model.Q
+        self.latent_dim = model.input_dim
 
         # The red cross which shows current latent point.
         self.latent_values = vals
@@ -113,7 +113,7 @@ class lvm(data_show):
         # A click in the bar chart axis for selection a dimension.
         if self.sense_axes != None:
             self.sense_axes.cla()
-            self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
+            self.sense_axes.bar(np.arange(self.model.input_dim),1./self.model.input_sensitivity(),color='b')
 
             if self.latent_index[1] == self.latent_index[0]:
                 self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
@@ -128,11 +128,11 @@ class lvm(data_show):
 
 class lvm_subplots(lvm):
     """
-    latent_axes is a np array of dimension np.ceil(Q/2),
+    latent_axes is a np array of dimension np.ceil(input_dim/2),
     one for each pair of the latent dimensions.
     """
     def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None):
-        self.nplots = int(np.ceil(model.Q/2.))+1
+        self.nplots = int(np.ceil(model.input_dim/2.))+1
         assert len(latent_axes)==self.nplots
         if vals==None:
             vals = model.X[0, :]
@@ -140,7 +140,7 @@ class lvm_subplots(lvm):
 
         for i, axis in enumerate(latent_axes):
             if i == self.nplots-1:
-                if self.nplots*2!=model.Q:
+                if self.nplots*2!=model.input_dim:
                     latent_index = [i*2, i*2]
                 lvm.__init__(self, self.latent_vals, model, data_visualize, axis, sense_axes, latent_index=latent_index)
             else:
@@ -174,7 +174,7 @@ class lvm_dimselect(lvm):
     def on_click(self, event):
 
         if event.inaxes==self.sense_axes:
-            new_index = max(0,min(int(np.round(event.xdata-0.5)),self.model.Q-1))
+            new_index = max(0,min(int(np.round(event.xdata-0.5)),self.model.input_dim-1))
             if event.button == 1:
                 # Make it red if and y-axis (red=port=left) if it is a left button click
                 self.latent_index[1] = new_index