diff --git a/GPy/core/model.py b/GPy/core/model.py
index 46cf6ac9..d3c6e582 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -14,14 +14,14 @@ from ..inference import optimization
 class model(parameterised):
     def __init__(self):
         parameterised.__init__(self)
-        self.priors = [None for i in range(self.get_param().size)]
+        self.priors = [None for i in range(self._get_params().size)]
         self.optimization_runs = []
         self.sampling_runs = []
-        self.set_param(self.get_param())
+        self._set_params(self._get_params())
         self.preferred_optimizer = 'tnc'
-    def get_param(self):
+    def _get_params(self):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
-    def set_param(self,x):
+    def _set_params(self,x):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
     def log_likelihood(self):
         raise NotImplementedError, "this needs to be implemented to utilise the model class"
@@ -67,7 +67,7 @@ class model(parameterised):
             unconst = np.setdiff1d(which, self.constrained_positive_indices)
             if len(unconst):
                 print "Warning: constraining parameters to be positive:"
-                print '\n'.join([n for i,n in enumerate(self.get_param_names()) if i in unconst])
+                print '\n'.join([n for i,n in enumerate(self._get_param_names()) if i in unconst])
                 print '\n'
                 self.constrain_positive(unconst)
         elif isinstance(what,priors.Gaussian):
@@ -86,7 +86,7 @@ class model(parameterised):
         """
         matches = self.grep_param_names(name)
         if len(matches):
-            return self.get_param()[matches]
+            return self._get_params()[matches]
         else:
             raise AttributeError, "no parameter matches %s"%name
 
@@ -96,9 +96,9 @@ class model(parameterised):
         """
         matches = self.grep_param_names(name)
         if len(matches):
-            x = self.get_param()
+            x = self._get_params()
             x[matches] = val
-            self.set_param(x)
+            self._set_params(x)
         else:
             raise AttributeError, "no parameter matches %s"%name
 
@@ -106,11 +106,11 @@ class model(parameterised):
 
     def log_prior(self):
         """evaluate the prior"""
-        return np.sum([p.lnpdf(x) for p, x in zip(self.priors,self.get_param()) if p is not None])
+        return np.sum([p.lnpdf(x) for p, x in zip(self.priors,self._get_params()) if p is not None])
 
     def log_prior_gradients(self):
         """evaluate the gradients of the priors"""
-        x = self.get_param()
+        x = self._get_params()
         ret = np.zeros(x.size)
         [np.put(ret,i,p.lnpdf_grad(xx)) for i,(p,xx) in enumerate(zip(self.priors,x)) if not p is None]
         return ret
@@ -121,7 +121,7 @@ class model(parameterised):
         Adjust the gradient for constraints and ties, return.
         """
         g = self.log_likelihood_gradients() + self.log_prior_gradients()
-        x = self.get_param()
+        x = self._get_params()
         g[self.constrained_positive_indices] = g[self.constrained_positive_indices]*x[self.constrained_positive_indices]
         g[self.constrained_negative_indices] = g[self.constrained_negative_indices]*x[self.constrained_negative_indices]
         [np.put(g,i,g[i]*(x[i]-l)*(h-x[i])/(h-l)) for i,l,h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
@@ -142,9 +142,9 @@ class model(parameterised):
         x = np.random.randn(x.size)
         self.expand_param(x)
         #now draw from prior where possible
-        x = self.get_param()
+        x = self._get_params()
         [np.put(x,i,p.rvs(1)) for i,p in enumerate(self.priors) if not p is None]
-        self.set_param(x)
+        self._set_params(x)
         self.expand_param(self.extract_param())#makes sure all of the tied parameters get the same init (since there's only one prior object...)
 
 
@@ -194,7 +194,7 @@ class model(parameterised):
         for s in positive_strings:
             for i in self.grep_param_names(s):
                 if not (i in self.all_constrained_indices()):
-                    name = self.get_param_names()[i]
+                    name = self._get_param_names()[i]
                     self.constrain_positive(name)
                     if warn:
                         print "Warning! constraining %s postive"%name
@@ -248,13 +248,13 @@ class model(parameterised):
 
         else:
             print "numerically calculating hessian. please be patient!"
-            x = self.get_param()
+            x = self._get_params()
             def f(x):
-                self.set_param(x)
+                self._set_params(x)
                 return self.log_likelihood()
             h = ndt.Hessian(f)
             A = -h(x)
-            self.set_param(x)
+            self._set_params(x)
         # check for almost zero components on the diagonal which screw up the cholesky
         aa = np.nonzero((np.diag(A)<1e-6) & (np.diag(A)>0.))[0]
         A[aa,aa] = 0.
@@ -268,7 +268,7 @@ class model(parameterised):
             hld = np.sum(np.log(np.diag(jitchol(A)[0])))
         except:
             return np.nan
-        return 0.5*self.get_param().size*np.log(2*np.pi) + self.log_likelihood() - hld
+        return 0.5*self._get_params().size*np.log(2*np.pi) + self.log_likelihood() - hld
 
     def __str__(self):
         s = parameterised.__str__(self).split('\n')
diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py
index da0b6056..738bde5b 100644
--- a/GPy/core/parameterised.py
+++ b/GPy/core/parameterised.py
@@ -87,7 +87,7 @@ class parameterised(object):
 
         Returns
         -------
-        the indices of self.get_param_names which match the regular expression.
+        the indices of self._get_param_names which match the regular expression.
 
         Notes
         -----
@@ -96,9 +96,9 @@ class parameterised(object):
 
         if type(expr) is str:
             expr = re.compile(expr)
-            return np.nonzero([expr.search(name) for name in self.get_param_names()])[0]
+            return np.nonzero([expr.search(name) for name in self._get_param_names()])[0]
         elif type(expr) is re._pattern_type:
-            return np.nonzero([expr.search(name) for name in self.get_param_names()])[0]
+            return np.nonzero([expr.search(name) for name in self._get_param_names()])[0]
         else:
             return expr
 
@@ -115,11 +115,11 @@ class parameterised(object):
         assert not np.any(matches[:,None]==self.all_constrained_indices()), "Some indices are already constrained"
         self.constrained_positive_indices = np.hstack((self.constrained_positive_indices, matches))
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if (xx<0) & (i in matches):
                 x[i] = -xx
-        self.set_param(x)
+        self._set_params(x)
 
 
     def unconstrain(self,which):
@@ -163,11 +163,11 @@ class parameterised(object):
         assert not np.any(matches[:,None]==self.all_constrained_indices()), "Some indices are already constrained"
         self.constrained_negative_indices = np.hstack((self.constrained_negative_indices, matches))
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if (xx>0.) and (i in matches):
                 x[i] = -xx
-        self.set_param(x)
+        self._set_params(x)
 
 
 
@@ -187,11 +187,11 @@ class parameterised(object):
         self.constrained_bounded_uppers.append(upper)
         self.constrained_bounded_lowers.append(lower)
         #check to ensure constraint is in place
-        x = self.get_param()
+        x = self._get_params()
         for i,xx in enumerate(x):
             if ((xx<=lower)|(xx>=upper)) & (i in matches):
                 x[i] = sigmoid(xx)*(upper-lower) + lower
-        self.set_param(x)
+        self._set_params(x)
 
 
     def constrain_fixed(self, which, value = None):
@@ -213,14 +213,14 @@ class parameterised(object):
         if value != None:
             self.constrained_fixed_values.append(value)
         else:
-            self.constrained_fixed_values.append(self.get_param()[self.constrained_fixed_indices[-1]])
+            self.constrained_fixed_values.append(self._get_params()[self.constrained_fixed_indices[-1]])
 
         #self.constrained_fixed_values.append(value)
         self.expand_param(self.extract_param())
 
     def extract_param(self):
-        """use self.get_param to get the 'true' parameters of the model, which are then tied, constrained and fixed"""
-        x = self.get_param()
+        """use self._get_params to get the 'true' parameters of the model, which are then tied, constrained and fixed"""
+        x = self._get_params()
         x[self.constrained_positive_indices] = np.log(x[self.constrained_positive_indices])
         x[self.constrained_negative_indices] = np.log(-x[self.constrained_negative_indices])
         [np.put(x,i,np.log(np.clip(x[i]-l,1e-10,np.inf)/np.clip(h-x[i],1e-10,np.inf))) for i,l,h in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
@@ -233,7 +233,7 @@ class parameterised(object):
 
 
     def expand_param(self,x):
-        """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self.set_param"""
+        """ takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params"""
 
         #work out how many places are fixed, and where they are. tricky logic!
         Nfix_places = 0.
@@ -257,14 +257,14 @@ class parameterised(object):
         xx[self.constrained_positive_indices] = np.exp(xx[self.constrained_positive_indices])
         xx[self.constrained_negative_indices] = -np.exp(xx[self.constrained_negative_indices])
         [np.put(xx,i,low+sigmoid(xx[i])*(high-low)) for i,low,high in zip(self.constrained_bounded_indices, self.constrained_bounded_lowers, self.constrained_bounded_uppers)]
-        self.set_param(xx)
+        self._set_params(xx)
 
     def extract_param_names(self):
         """
         Returns the parameter names as propagated after constraining,
         tying or fixing, i.e. a list of the same length as extract_param()
         """
-        n =  self.get_param_names()
+        n =  self._get_param_names()
 
         #remove/concatenate the tied parameter names
         if len(self.tied_indices):
@@ -294,13 +294,13 @@ class parameterised(object):
         """
         Return a string describing the parameter names and their ties and constraints
         """
-        names = self.get_param_names()
+        names = self._get_param_names()
         N = len(names)
 
         if not N:
             return "This object has no free parameters."
         header = ['Name','Value','Constraints','Ties']
-        values = self.get_param() #map(str,self.get_param())
+        values = self._get_params() #map(str,self._get_params())
         #sort out the constraints
         constraints = ['']*len(names)
         for i in self.constrained_positive_indices:
diff --git a/GPy/examples/GPLVM_demo.py b/GPy/examples/GPLVM_demo.py
deleted file mode 100644
index 651af6be..00000000
--- a/GPy/examples/GPLVM_demo.py
+++ /dev/null
@@ -1,28 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-import pylab as pb
-import GPy
-np.random.seed(1)
-print "GPLVM with RBF kernel"
-
-N = 100
-Q = 1
-D = 2
-X = np.random.rand(N, Q)
-k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001)
-K = k.K(X)
-Y = np.random.multivariate_normal(np.zeros(N),K,D).T
-
-m = GPy.models.GPLVM(Y, Q)
-m.constrain_positive('(rbf|bias|white)')
-
-pb.figure()
-m.plot()
-pb.title('PCA initialisation')
-pb.figure()
-m.optimize(messages = 1)
-m.plot()
-pb.title('After optimisation')
diff --git a/GPy/examples/GP_regression_demo.py b/GPy/examples/GP_regression_demo.py
deleted file mode 100644
index 0fe2fabd..00000000
--- a/GPy/examples/GP_regression_demo.py
+++ /dev/null
@@ -1,51 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-"""
-Simple Gaussian Processes regression with an RBF kernel
-"""
-import pylab as pb
-import numpy as np
-import GPy
-pb.ion()
-pb.close('all')
-
-
-######################################
-## 1 dimensional example
-
-# sample inputs and outputs
-X = np.random.uniform(-3.,3.,(20,1))
-Y = np.sin(X)+np.random.randn(20,1)*0.05
-
-# create simple GP model
-m = GPy.models.GP_regression(X,Y)
-
-# contrain all parameters to be positive
-m.constrain_positive('')
-
-# optimize and plot
-m.optimize('tnc', max_f_eval = 1000)
-m.plot()
-print(m)
-
-######################################
-## 2 dimensional example
-
-# sample inputs and outputs
-X = np.random.uniform(-3.,3.,(40,2))
-Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(40,1)*0.05
-
-# create simple GP model
-m = GPy.models.GP_regression(X,Y)
-
-# contrain all parameters to be positive
-m.constrain_positive('')
-# optimize and plot
-pb.figure()
-m.optimize('tnc', max_f_eval = 1000)
-m.plot()
-print(m)
-
-
diff --git a/GPy/examples/GP_regression_kern_demo.py b/GPy/examples/GP_regression_kern_demo.py
deleted file mode 100644
index 8f0b9226..00000000
--- a/GPy/examples/GP_regression_kern_demo.py
+++ /dev/null
@@ -1,33 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-"""
-Simple one-dimensional Gaussian Processes with assorted kernel functions
-"""
-import pylab as pb
-import numpy as np
-import GPy
-
-# sample inputs and outputs
-D = 1
-X = np.random.randn(10,D)*2
-X = np.linspace(-1.5,1.5,5)[:,None]
-X = np.append(X,[[5]],0)
-Y = np.sin(np.pi*X/2) #+np.random.randn(X.shape[0],1)*0.05
-
-models = [GPy.models.GP_regression(X,Y, k) for k in (GPy.kern.rbf(D), GPy.kern.Matern52(D), GPy.kern.Matern32(D), GPy.kern.exponential(D), GPy.kern.linear(D) + GPy.kern.white(D),  GPy.kern.bias(D) + GPy.kern.white(D))]
-
-pb.figure(figsize=(12,8))
-for i,m in enumerate(models):
-    m.constrain_positive('')
-    m.optimize()
-    pb.subplot(3,2,i+1)
-    m.plot()
-    #pb.title(m.kern.parts[0].name)
-
-GPy.util.plot.align_subplots(3,2,(-3,6),(-2.5,2.5))
-
-pb.show()
-
-
diff --git a/GPy/examples/classification.py b/GPy/examples/classification.py
index eb7a887e..989ed08a 100644
--- a/GPy/examples/classification.py
+++ b/GPy/examples/classification.py
@@ -8,8 +8,6 @@ Simple Gaussian Processes classification
 import pylab as pb
 import numpy as np
 import GPy
-pb.ion()
-pb.close('all')
 
 default_seed=10000
 ######################################
@@ -27,7 +25,7 @@ def crescent_data(model_type='Full', inducing=10, seed=default_seed):
     likelihood = GPy.inference.likelihoods.probit(data['Y'])
 
     if model_type=='Full':
-        m = GPy.models.simple_GP_EP(data['X'],likelihood)
+        m = GPy.models.GP_EP(data['X'],likelihood)
     else:
         # create sparse GP EP model
         m = GPy.models.sparse_GP_EP(data['X'],likelihood=likelihood,inducing=inducing,ep_proxy=model_type)
@@ -49,7 +47,7 @@ def oil():
     likelihood = GPy.inference.likelihoods.probit(data['Y'][:, 0:1])
 
     # create simple GP model
-    m = GPy.models.simple_GP_EP(data['X'],likelihood)
+    m = GPy.models.GP_EP(data['X'],likelihood)
 
     # contrain all parameters to be positive
     m.constrain_positive('')
diff --git a/GPy/examples/oil_flow_demo.py b/GPy/examples/oil_flow_demo.py
deleted file mode 100644
index eee74461..00000000
--- a/GPy/examples/oil_flow_demo.py
+++ /dev/null
@@ -1,53 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import cPickle as pickle
-import numpy as np
-import pylab as pb
-import GPy
-import pylab as plt
-np.random.seed(1)
-
-def plot_oil(X, theta, labels, label):
-    plt.figure()
-    X = X[:,np.argsort(theta)[:2]]
-    flow_type = (X[labels[:,0]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'rx')
-    flow_type = (X[labels[:,1]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'gx')
-    flow_type = (X[labels[:,2]==1])
-    plt.plot(flow_type[:,0], flow_type[:,1], 'bx')
-    plt.title(label)
-
-data = pickle.load(open('../util/datasets/oil_flow_3classes.pickle', 'r'))
-
-Y = data['DataTrn']
-N, D = Y.shape
-selected = np.random.permutation(N)[:200]
-labels = data['DataTrnLbls'][selected]
-Y = Y[selected]
-N, D = Y.shape
-Y -= Y.mean(axis=0)
-Y /= Y.std(axis=0)
-
-Q = 2
-m1 = GPy.models.sparse_GPLVM(Y, Q, M = 15)
-m1.constrain_positive('(rbf|bias|noise)')
-m1.constrain_bounded('white', 1e-6, 1.0)
-
-plot_oil(m1.X, np.array([1,1]), labels, 'PCA initialization')
-# m.optimize(messages = True)
-m1.optimize('bfgs', messages = True)
-plot_oil(m1.X, np.array([1,1]), labels, 'sparse GPLVM')
-# pb.figure()
-# m.plot()
-# pb.title('PCA initialisation')
-# pb.figure()
-# m.optimize(messages = 1)
-# m.plot()
-# pb.title('After optimisation')
-m = GPy.models.GPLVM(Y, Q)
-m.constrain_positive('(white|rbf|bias|noise)')
-m.optimize()
-plot_oil(m.X, np.array([1,1]), labels, 'GPLVM')
diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index d61c8b10..79763504 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -8,8 +8,6 @@ Gaussian Processes regression examples
 import pylab as pb
 import numpy as np
 import GPy
-pb.ion()
-pb.close('all')
 
 
 def toy_rbf_1d():
@@ -48,6 +46,10 @@ def rogers_girolami_olympics():
     print(m)
     return m
 
+def della_gatta_TRP63_gene_expression(number=942):
+    """Run a standard Gaussian process regression on the della Gatta et al TRP63 Gene Expression data set for a given gene number."""
+
+
 def toy_rbf_1d_50():
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
     data = GPy.util.datasets.toy_rbf_1d_50()
@@ -81,3 +83,94 @@ def silhouette():
 
     print(m)
     return m
+
+
+def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000):
+    """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""
+
+    # Contour over a range of length scales and signal/noise ratios.
+    length_scales = np.linspace(0.1, 60., resolution)
+    log_SNRs = np.linspace(-3., 4., resolution)
+
+    data = GPy.util.datasets.della_gatta_TRP63_gene_expression(gene_number)
+    # Sub sample the data to ensure multiple optima
+    #data['Y'] = data['Y'][0::2, :]
+    #data['X'] = data['X'][0::2, :]
+
+    # Remove the mean (no bias kernel to ensure signal/noise is in RBF/white)
+    data['Y'] = data['Y'] - np.mean(data['Y'])
+
+    lls = GPy.examples.regression.contour_data(data, length_scales, log_SNRs, GPy.kern.rbf)
+    pb.contour(length_scales, log_SNRs, np.exp(lls), 20)
+    ax = pb.gca()
+    pb.xlabel('length scale')
+    pb.ylabel('log_10 SNR')
+
+    xlim = ax.get_xlim()
+    ylim = ax.get_ylim()
+    
+    # Now run a few optimizations
+    models = []
+    optim_point_x = np.empty(2)
+    optim_point_y = np.empty(2)
+    np.random.seed(seed=seed)
+    for i in range(0, model_restarts):
+        kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.))
+        
+        m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern)
+        params = m._get_params()
+        optim_point_x[0] = params[1]
+        optim_point_y[0] = np.log10(params[0]) - np.log10(params[2]);
+        
+        # contrain all parameters to be positive
+        m.constrain_positive('')
+        
+        # optimize
+        m.optimize(xtol=1e-6,ftol=1e-6)
+
+        params = m._get_params()
+        optim_point_x[1] = params[1]
+        optim_point_y[1] = np.log10(params[0]) - np.log10(params[2]);
+        print(m)
+        
+        pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
+        models.append(m)
+
+    ax.set_xlim(xlim)
+    ax.set_ylim(ylim)
+    return (models, lls)
+
+def contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf):
+    """Evaluate the GP objective function for a given data set for a range of signal to noise ratios and a range of lengthscales.
+
+    :data_set: A data set from the utils.datasets director.
+    :length_scales: a list of length scales to explore for the contour plot.
+    :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot.
+    :signal_kernel: a kernel to use for the 'signal' portion of the data."""
+
+    lls = []
+    total_var = np.var(data['Y'])
+    for log_SNR in log_SNRs:
+        SNR = 10**log_SNR
+        length_scale_lls = []
+        for length_scale in length_scales:
+            noise_var = 1.
+            signal_var = SNR
+            noise_var = noise_var/(noise_var + signal_var)*total_var
+            signal_var = signal_var/(noise_var + signal_var)*total_var
+
+            signal_kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale)
+            noise_kernel = GPy.kern.white(1, variance=noise_var)
+            kernel = signal_kernel + noise_kernel
+            K = kernel.K(data['X'])
+            total_var = (np.dot(np.dot(data['Y'].T,GPy.util.linalg.pdinv(K)[0]), data['Y'])/data['Y'].shape[0])[0,0]
+            noise_var *= total_var
+            signal_var *= total_var
+            
+            kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale) + GPy.kern.white(1, variance=noise_var)
+
+            model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel)
+            model.constrain_positive('')
+            length_scale_lls.append(model.log_likelihood())
+        lls.append(length_scale_lls)
+    return np.array(lls)
diff --git a/GPy/examples/warped_GP_demo.py b/GPy/examples/warped_GP_demo.py
deleted file mode 100644
index 3b75694a..00000000
--- a/GPy/examples/warped_GP_demo.py
+++ /dev/null
@@ -1,45 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-
-import numpy as np
-import scipy as sp
-import pdb, sys, pickle
-import matplotlib.pylab as plt
-import GPy
-np.random.seed(1)
-
-N = 100
-# sample inputs and outputs
-X = np.random.uniform(-np.pi,np.pi,(N,1))
-Y = np.sin(X)+np.random.randn(N,1)*0.05
-# Y += np.abs(Y.min()) + 0.5
-Z = np.exp(Y)# Y**(1/3.0)
-
-# rescaling targets?
-Zmax = Z.max()
-Zmin = Z.min()
-Z = (Z-Zmin)/(Zmax-Zmin) - 0.5
-
-m = GPy.models.warpedGP(X, Z, warping_terms = 2)
-m.constrain_positive('(tanh_a|tanh_b|rbf|white|bias)')
-m.randomize()
-plt.figure()
-plt.xlabel('predicted f(Z)')
-plt.ylabel('actual f(Z)')
-plt.plot(m.Y, Y, 'o', alpha = 0.5, label = 'before training')
-m.optimize(messages = True)
-plt.plot(m.Y, Y, 'o', alpha = 0.5, label = 'after training')
-plt.legend(loc = 0)
-m.plot_warping()
-plt.figure()
-plt.title('warped GP fit')
-m.plot()
-
-m1 = GPy.models.GP_regression(X, Z)
-m1.constrain_positive('(rbf|white|bias)')
-m1.randomize()
-m1.optimize(messages = True)
-plt.figure()
-plt.title('GP fit')
-m1.plot()
diff --git a/GPy/inference/samplers.py b/GPy/inference/samplers.py
index 1216f1eb..da159c23 100644
--- a/GPy/inference/samplers.py
+++ b/GPy/inference/samplers.py
@@ -73,12 +73,12 @@ class Metropolis_Hastings:
 
     def predict(self,function,args):
         """Make a prediction for the function, to which we will pass the additional arguments"""
-        param = self.model.get_param()
+        param = self.model._get_params()
         fs = []
         for p in self.chain:
-            self.model.set_param(p)
+            self.model._set_params(p)
             fs.append(function(*args))
-        self.model.set_param(param)# reset model to starting state
+        self.model._set_params(param)# reset model to starting state
         return fs
 
 
diff --git a/GPy/kern/Brownian.py b/GPy/kern/Brownian.py
index 07d5fbd6..fd19137c 100644
--- a/GPy/kern/Brownian.py
+++ b/GPy/kern/Brownian.py
@@ -23,16 +23,16 @@ class Brownian(kernpart):
         assert self.D==1, "Brownian motion in 1D only"
         self.Nparam = 1.
         self.name = 'Brownian'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):
diff --git a/GPy/kern/Matern32.py b/GPy/kern/Matern32.py
index 377b687f..6517ac2c 100644
--- a/GPy/kern/Matern32.py
+++ b/GPy/kern/Matern32.py
@@ -34,19 +34,19 @@ class Matern32(kernpart):
             lengthscales = np.ones(self.D)
         self.Nparam = self.D + 1
         self.name = 'Mat32'
-        self.set_param(np.hstack((variance,lengthscales)))
+        self._set_params(np.hstack((variance,lengthscales)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscales))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==(self.D+1)
         self.variance = x[0]
         self.lengthscales = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
         if self.D==1:
             return ['variance','lengthscale']
diff --git a/GPy/kern/Matern52.py b/GPy/kern/Matern52.py
index 172c0117..056da32d 100644
--- a/GPy/kern/Matern52.py
+++ b/GPy/kern/Matern52.py
@@ -32,19 +32,19 @@ class Matern52(kernpart):
             lengthscales = np.ones(self.D)
         self.Nparam = self.D + 1
         self.name = 'Mat52'
-        self.set_param(np.hstack((variance,lengthscales)))
+        self._set_params(np.hstack((variance,lengthscales)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscales))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==(self.D+1)
         self.variance = x[0]
         self.lengthscales = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
         if self.D==1:
             return ['variance','lengthscale']
diff --git a/GPy/kern/bias.py b/GPy/kern/bias.py
index a89f56c5..0a6872a5 100644
--- a/GPy/kern/bias.py
+++ b/GPy/kern/bias.py
@@ -17,16 +17,16 @@ class bias(kernpart):
         self.D = D
         self.Nparam = 1
         self.name = 'bias'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):
diff --git a/GPy/kern/exponential.py b/GPy/kern/exponential.py
index 8ad5b813..2df6a958 100644
--- a/GPy/kern/exponential.py
+++ b/GPy/kern/exponential.py
@@ -32,19 +32,19 @@ class exponential(kernpart):
             lengthscales = np.ones(self.D)
         self.Nparam = self.D + 1
         self.name = 'exp'
-        self.set_param(np.hstack((variance,lengthscales)))
+        self._set_params(np.hstack((variance,lengthscales)))
 
-    def get_param(self):
+    def _get_params(self):
         """return the value of the parameters."""
         return np.hstack((self.variance,self.lengthscales))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         """set the value of the parameters."""
         assert x.size==(self.D+1)
         self.variance = x[0]
         self.lengthscales = x[1:]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         """return parameter names."""
         if self.D==1:
             return ['variance','lengthscale']
diff --git a/GPy/kern/finite_dimensional.py b/GPy/kern/finite_dimensional.py
index 98c99628..36afd1dc 100644
--- a/GPy/kern/finite_dimensional.py
+++ b/GPy/kern/finite_dimensional.py
@@ -27,15 +27,15 @@ class finite_dimensional(kernpart):
             weights = np.ones(self.n)
         self.Nparam = self.n + 1
         self.name = 'finite_dim'
-        self.set_param(np.hstack((variance,weights)))
+        self._set_params(np.hstack((variance,weights)))
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.variance,self.weights))
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.size == (self.Nparam)
         self.variance = x[0]
         self.weights = x[1:]
-    def get_param_names(self):
+    def _get_param_names(self):
         if self.n==1:
             return ['variance','weight']
         else:
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 3ba7d97b..5f259f55 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -133,20 +133,20 @@ class kern(parameterised):
         newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices]
         return newkern
 
-    def get_param(self):
-        return np.hstack([p.get_param() for p in self.parts])
+    def _get_params(self):
+        return np.hstack([p._get_params() for p in self.parts])
 
-    def set_param(self,x):
-        [p.set_param(x[s]) for p, s in zip(self.parts, self.param_slices)]
+    def _set_params(self,x):
+        [p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
 
-    def get_param_names(self):
+    def _get_param_names(self):
         #this is a bit nasty: we wat to distinguish between parts with the same name by appending a count
         part_names = np.array([k.name for k in self.parts],dtype=np.str)
         counts = [np.sum(part_names==ni) for i, ni in enumerate(part_names)]
         cum_counts = [np.sum(part_names[i:]==ni) for i, ni in enumerate(part_names)]
         names = [name+'_'+str(cum_count) if count>1 else name for name,count,cum_count in zip(part_names,counts,cum_counts)]
 
-        return sum([[name+'_'+n for n in k.get_param_names()] for name,k in zip(names,self.parts)],[])
+        return sum([[name+'_'+n for n in k._get_param_names()] for name,k in zip(names,self.parts)],[])
 
     def K(self,X,X2=None,slices1=None,slices2=None):
         assert X.shape[1]==self.D
diff --git a/GPy/kern/kernpart.py b/GPy/kern/kernpart.py
index d06749b8..3a5486de 100644
--- a/GPy/kern/kernpart.py
+++ b/GPy/kern/kernpart.py
@@ -16,11 +16,11 @@ class kernpart(object):
         self.Nparam = 1
         self.name = 'unnamed'
 
-    def get_param(self):
+    def _get_params(self):
         raise NotImplementedError
-    def set_param(self,x):
+    def _set_params(self,x):
         raise NotImplementedError
-    def get_param_names(self):
+    def _get_param_names(self):
         raise NotImplementedError
     def K(self,X,X2,target):
         raise NotImplementedError
diff --git a/GPy/kern/linear.py b/GPy/kern/linear.py
index 7246244e..14c0ece3 100644
--- a/GPy/kern/linear.py
+++ b/GPy/kern/linear.py
@@ -20,16 +20,16 @@ class linear(kernpart):
             variance = 1.0
         self.Nparam = 1
         self.name = 'linear'
-        self.set_param(variance)
+        self._set_params(variance)
         self._Xcache, self._X2cache = np.empty(shape=(2,))
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):
diff --git a/GPy/kern/linear_ARD.py b/GPy/kern/linear_ARD.py
index 149bfb4d..b9112044 100644
--- a/GPy/kern/linear_ARD.py
+++ b/GPy/kern/linear_ARD.py
@@ -23,16 +23,16 @@ class linear_ARD(kernpart):
             variances = np.ones(self.D)
         self.Nparam = int(self.D)
         self.name = 'linear'
-        self.set_param(variances)
+        self._set_params(variances)
 
-    def get_param(self):
+    def _get_params(self):
         return self.variances
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.size==(self.Nparam)
         self.variances = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         if self.D==1:
             return ['variance']
         else:
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index 9e2bb509..c929369f 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -30,23 +30,23 @@ class rbf(kernpart):
         self.D = D
         self.Nparam = 2
         self.name = 'rbf'
-        self.set_param(np.hstack((variance,lengthscale)))
+        self._set_params(np.hstack((variance,lengthscale)))
 
         #initialize cache
         self._Z, self._mu, self._S = np.empty(shape=(3,1))
         self._X, self._X2, self._params = np.empty(shape=(3,1))
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.variance,self.lengthscale))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.variance, self.lengthscale = x
         self.lengthscale2 = np.square(self.lengthscale)
         #reset cached results
         self._X, self._X2, self._params = np.empty(shape=(3,1))
         self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance','lengthscale']
 
     def K(self,X,X2,target):
diff --git a/GPy/kern/rbf_ARD.py b/GPy/kern/rbf_ARD.py
index 732be590..12924b3f 100644
--- a/GPy/kern/rbf_ARD.py
+++ b/GPy/kern/rbf_ARD.py
@@ -22,16 +22,16 @@ class rbf_ARD(kernpart):
             lengthscales = np.ones(self.D)
         self.Nparam = self.D + 1
         self.name = 'rbf_ARD'
-        self.set_param(np.hstack((variance,lengthscales)))
+        self._set_params(np.hstack((variance,lengthscales)))
 
         #initialize cache
         self._Z, self._mu, self._S = np.empty(shape=(3,1))
         self._X, self._X2, self._params = np.empty(shape=(3,1))
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.variance,self.lengthscales))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.size==(self.D+1)
         self.variance = x[0]
         self.lengthscales = x[1:]
@@ -39,7 +39,7 @@ class rbf_ARD(kernpart):
         #reset cached results
         self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
 
-    def get_param_names(self):
+    def _get_param_names(self):
         if self.D==1:
             return ['variance','lengthscale']
         else:
@@ -135,8 +135,8 @@ class rbf_ARD(kernpart):
             if X2 is None: X2 = X
             self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy
             self._params = np.empty(shape=(1,0))#ensure the next section gets called
-        if not np.all(self._params == self.get_param()):
-            self._params == self.get_param()
+        if not np.all(self._params == self._get_params()):
+            self._params == self._get_params()
             self._K_dist2 = np.square(self._K_dist/self.lengthscales)
             self._K_exponent = -0.5*self._K_dist2.sum(-1)
             self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1))
@@ -189,7 +189,7 @@ if __name__=='__main__':
     from checkgrad import checkgrad
 
     def k_theta_test(param,k):
-        k.set_param(param)
+        k._set_params(param)
         K = k.K(Z)
         dK_dtheta = k.dK_dtheta(Z)
         f = np.sum(K)
@@ -216,7 +216,7 @@ if __name__=='__main__':
     checkgrad(psi1_S_test,np.random.rand(N*Q),args=(k,))
 
     def psi1_theta_test(theta,k):
-        k.set_param(theta)
+        k._set_params(theta)
         f = np.sum(k.psi1(Z,mu,S))
         df = np.array([np.sum(grad) for grad in k.dpsi1_dtheta(Z,mu,S)])
         return f,df
@@ -241,7 +241,7 @@ if __name__=='__main__':
     checkgrad(psi2_S_test,np.random.rand(N*Q),args=(k,))
 
     def psi2_theta_test(theta,k):
-        k.set_param(theta)
+        k._set_params(theta)
         f = np.sum(k.psi2(Z,mu,S))
         df = np.array([np.sum(grad) for grad in k.dpsi2_dtheta(Z,mu,S)])
         return f,df
diff --git a/GPy/kern/spline.py b/GPy/kern/spline.py
index 44033eea..030b2f02 100644
--- a/GPy/kern/spline.py
+++ b/GPy/kern/spline.py
@@ -25,15 +25,15 @@ class spline(kernpart):
         assert self.D==1
         self.Nparam = 1
         self.name = 'spline'
-        self.set_param(np.squeeze(variance))
+        self._set_params(np.squeeze(variance))
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):
diff --git a/GPy/kern/sympykern.py b/GPy/kern/sympykern.py
index d9f89f5b..db3cc976 100644
--- a/GPy/kern/sympykern.py
+++ b/GPy/kern/sympykern.py
@@ -44,7 +44,7 @@ class spkern(kernpart):
         if param is None:
             param = np.ones(self.Nparam)
         assert param.size==self.Nparam
-        self.set_param(param)
+        self._set_params(param)
 
         #Differentiate!
         self._sp_dk_dtheta = [sp.diff(k,theta).simplify() for theta in self._sp_theta]
@@ -247,12 +247,12 @@ class spkern(kernpart):
         Z = X
         weave.inline(self._dKdiag_dX_code,arg_names=['target','X','Z','param','partial'],**self.weave_kwargs)
 
-    def set_param(self,param):
+    def _set_params(self,param):
         #print param.flags['C_CONTIGUOUS']
         self._param = param.copy()
 
-    def get_param(self):
+    def _get_params(self):
         return self._param
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return [x.name for x in self._sp_theta]
diff --git a/GPy/kern/white.py b/GPy/kern/white.py
index 587a2b4a..b3b00c48 100644
--- a/GPy/kern/white.py
+++ b/GPy/kern/white.py
@@ -17,16 +17,16 @@ class white(kernpart):
         self.D = D
         self.Nparam = 1
         self.name = 'white'
-        self.set_param(np.array([variance]).flatten())
+        self._set_params(np.array([variance]).flatten())
 
-    def get_param(self):
+    def _get_params(self):
         return self.variance
 
-    def set_param(self,x):
+    def _set_params(self,x):
         assert x.shape==(1,)
         self.variance = x
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return ['variance']
 
     def K(self,X,X2,target):
diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py
index 44147b73..0a168f25 100644
--- a/GPy/models/GPLVM.py
+++ b/GPy/models/GPLVM.py
@@ -33,16 +33,16 @@ class GPLVM(GP_regression):
         else:
             return np.random.randn(Y.shape[0], Q)
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return (sum([['X_%i_%i'%(n,q) for n in range(self.N)] for q in range(self.Q)],[])
                 + self.kern.extract_param_names())
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack((self.X.flatten(), self.kern.extract_param()))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
-        GP_regression.set_param(self, x[self.X.size:])
+        GP_regression._set_params(self, x[self.X.size:])
 
     def log_likelihood_gradients(self):
         dL_dK = self.dL_dK()
diff --git a/GPy/models/GP_EP.py b/GPy/models/GP_EP.py
index bb582674..5cb2ba45 100644
--- a/GPy/models/GP_EP.py
+++ b/GPy/models/GP_EP.py
@@ -41,13 +41,13 @@ class GP_EP(model):
         self.K = self.kernel.K(self.X)
         model.__init__(self)
 
-    def set_param(self,p):
+    def _set_params(self,p):
         self.kernel.expand_param(p)
 
-    def get_param(self):
+    def _get_params(self):
         return self.kernel.extract_param()
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return self.kernel.extract_param_names()
 
     def approximate_likelihood(self):
@@ -138,7 +138,7 @@ class GP_EP(model):
         """
         self.epsilon_em = epsilon
         log_likelihood_change = self.epsilon_em + 1.
-        self.parameters_path = [self.kernel.get_param()]
+        self.parameters_path = [self.kernel._get_params()]
         self.approximate_likelihood()
         self.site_approximations_path = [[self.ep_approx.tau_tilde,self.ep_approx.v_tilde]]
         self.log_likelihood_path = [self.log_likelihood()]
@@ -155,6 +155,6 @@ class GP_EP(model):
             else:
                 self.approximate_likelihood()
                 self.log_likelihood_path.append(self.log_likelihood())
-                self.parameters_path.append(self.kernel.get_param())
+                self.parameters_path.append(self.kernel._get_params())
                 self.site_approximations_path.append([self.ep_approx.tau_tilde,self.ep_approx.v_tilde])
             iteration += 1
diff --git a/GPy/models/GP_regression.py b/GPy/models/GP_regression.py
index ee2bb448..da657950 100644
--- a/GPy/models/GP_regression.py
+++ b/GPy/models/GP_regression.py
@@ -70,15 +70,15 @@ class GP_regression(model):
 
         model.__init__(self)
 
-    def set_param(self,p):
+    def _set_params(self,p):
         self.kern.expand_param(p)
         self.K = self.kern.K(self.X,slices1=self.Xslices)
         self.Ki, self.L, self.Li, self.K_logdet = pdinv(self.K)
 
-    def get_param(self):
+    def _get_params(self):
         return self.kern.extract_param()
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return self.kern.extract_param_names()
 
     def _model_fit_term(self):
diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py
index 4ef36f0e..40f18ff9 100644
--- a/GPy/models/generalized_FITC.py
+++ b/GPy/models/generalized_FITC.py
@@ -42,14 +42,14 @@ class generalized_FITC(model):
         self.jitter = 1e-12
         model.__init__(self)
 
-    def set_param(self,p):
+    def _set_params(self,p):
         self.kernel.expand_param(p[0:-self.Z.size])
         self.Z = p[-self.Z.size:].reshape(self.M,self.D)
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack([self.kernel.extract_param(),self.Z.flatten()])
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return self.kernel.extract_param_names()+['iip_%i'%i for i in range(self.Z.size)]
 
     def approximate_likelihood(self):
@@ -214,7 +214,7 @@ class generalized_FITC(model):
         """
         self.epsilon_em = epsilon
         log_likelihood_change = self.epsilon_em + 1.
-        self.parameters_path = [self.kernel.get_param()]
+        self.parameters_path = [self.kernel._get_params()]
         self.approximate_likelihood()
         self.site_approximations_path = [[self.ep_approx.tau_tilde,self.ep_approx.v_tilde]]
         self.inducing_inputs_path = [self.Z]
@@ -235,7 +235,7 @@ class generalized_FITC(model):
             else:
                 self.approximate_likelihood()
                 self.log_likelihood_path.append(self.log_likelihood())
-                self.parameters_path.append(self.kernel.get_param())
+                self.parameters_path.append(self.kernel._get_params())
                 self.site_approximations_path.append([self.ep_approx.tau_tilde,self.ep_approx.v_tilde])
                 self.inducing_inputs_path.append(self.Z)
             iteration += 1
diff --git a/GPy/models/sparse_GPLVM.py b/GPy/models/sparse_GPLVM.py
index c5125c85..1c6a03cf 100644
--- a/GPy/models/sparse_GPLVM.py
+++ b/GPy/models/sparse_GPLVM.py
@@ -27,16 +27,16 @@ class sparse_GPLVM(sparse_GP_regression, GPLVM):
         X = self.initialise_latent(init, Q, Y)
         sparse_GP_regression.__init__(self, X, Y, **kwargs)
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return (sum([['X_%i_%i'%(n,q) for n in range(self.N)] for q in range(self.Q)],[])
-                + sparse_GP_regression.get_param_names(self))
+                + sparse_GP_regression._get_param_names(self))
 
-    def get_param(self):
-        return np.hstack((self.X.flatten(), sparse_GP_regression.get_param(self)))
+    def _get_params(self):
+        return np.hstack((self.X.flatten(), sparse_GP_regression._get_params(self)))
 
-    def set_param(self,x):
+    def _set_params(self,x):
         self.X = x[:self.X.size].reshape(self.N,self.Q).copy()
-        sparse_GP_regression.set_param(self, x[self.X.size:])
+        sparse_GP_regression._set_params(self, x[self.X.size:])
 
     def log_likelihood(self):
         return sparse_GP_regression.log_likelihood(self)
diff --git a/GPy/models/sparse_GP_regression.py b/GPy/models/sparse_GP_regression.py
index fe5f7cc1..fc04d5ed 100644
--- a/GPy/models/sparse_GP_regression.py
+++ b/GPy/models/sparse_GP_regression.py
@@ -59,10 +59,10 @@ class sparse_GP_regression(GP_regression):
         if self.has_uncertain_inputs:
             self.X_uncertainty /= np.square(self._Xstd)
 
-    def set_param(self, p):
+    def _set_params(self, p):
         self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
         self.beta = p[self.M*self.Q]
-        self.kern.set_param(p[self.Z.size + 1:])
+        self.kern._set_params(p[self.Z.size + 1:])
         self.beta2 = self.beta**2
         self._compute_kernel_matrices()
         self._computations()
@@ -106,10 +106,10 @@ class sparse_GP_regression(GP_regression):
         self.dL_dKmm += -0.5 * self.D * (- self.LBL_inv - 2.*self.beta*mdot(self.LBL_inv, self.psi2, self.Kmmi) + self.Kmmi) # dC
         self.dL_dKmm +=  np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE
 
-    def get_param(self):
+    def _get_params(self):
         return np.hstack([self.Z.flatten(),self.beta,self.kern.extract_param()])
 
-    def get_param_names(self):
+    def _get_param_names(self):
         return sum([['iip_%i_%i'%(i,j) for i in range(self.Z.shape[0])] for j in range(self.Z.shape[1])],[]) + ['noise_precision']+self.kern.extract_param_names()
 
     def log_likelihood(self):
diff --git a/GPy/models/warped_GP.py b/GPy/models/warped_GP.py
index bf5af21f..ea3c44cc 100644
--- a/GPy/models/warped_GP.py
+++ b/GPy/models/warped_GP.py
@@ -13,7 +13,7 @@ from GP_regression import GP_regression
 
 class warpedGP(GP_regression):
     """
-    TODO: fucking docstrings!
+    TODO: fecking docstrings!
 
     @nfusi: I'#ve hacked a little on this, but no guarantees. J.
     """
@@ -30,17 +30,17 @@ class warpedGP(GP_regression):
         self.transform_data()
         GP_regression.__init__(self, X, self.Y, **kwargs)
 
-    def set_param(self, x):
+    def _set_params(self, x):
         self.warping_params = x[:self.warping_function.num_parameters].reshape(self.warp_params_shape).copy()
         self.transform_data()
-        GP_regression.set_param(self, x[self.warping_function.num_parameters:].copy())
+        GP_regression._set_params(self, x[self.warping_function.num_parameters:].copy())
 
-    def get_param(self):
-        return np.hstack((self.warping_params.flatten().copy(), GP_regression.get_param(self).copy()))
+    def _get_params(self):
+        return np.hstack((self.warping_params.flatten().copy(), GP_regression._get_params(self).copy()))
 
-    def get_param_names(self):
-        warping_names = self.warping_function.get_param_names()
-        param_names = GP_regression.get_param_names(self)
+    def _get_param_names(self):
+        warping_names = self.warping_function._get_param_names()
+        param_names = GP_regression._get_param_names(self)
         return warping_names + param_names
 
     def transform_data(self):
@@ -81,7 +81,7 @@ class warpedGP(GP_regression):
     def predict(self, X, in_unwarped_space = False, **kwargs):
         mu, var = GP_regression.predict(self, X, **kwargs)
 
-        # The plot() function calls set_param() before calling predict()
+        # The plot() function calls _set_params() before calling predict()
         # this is causing the observations to be plotted in the transformed
         # space (where Y lives), making the plot looks very wrong
         # if the predictions are made in the untransformed space
diff --git a/GPy/util/warping_functions.py b/GPy/util/warping_functions.py
index 5a6dc6b3..a87deb5b 100644
--- a/GPy/util/warping_functions.py
+++ b/GPy/util/warping_functions.py
@@ -33,7 +33,7 @@ class WarpingFunction(object):
         """inverse function transformation"""
         raise NotImplementedError
 
-    def get_param_names(self):
+    def _get_param_names(self):
         raise NotImplementedError
 
     def plot(self, psi, xmin, xmax):
@@ -151,7 +151,7 @@ class TanhWarpingFunction(WarpingFunction):
 
         return gradients
 
-    def get_param_names(self):
+    def _get_param_names(self):
         variables = ['a', 'b', 'c']
         names = sum([['warp_tanh_%s_t%i' % (variables[n],q) for n in range(3)] for q in range(self.n_terms)],[])
         return names