diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 89313ae2..898c7b58 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -89,7 +89,6 @@ class GP(Model):
             assert mean_function.output_dim == self.output_dim
             self.link_parameter(mean_function)
 
-
         #find a sensible inference method
         logger.info("initializing inference method")
         if inference_method is None:
diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index 1bc6a29e..48c05a8d 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -430,23 +430,38 @@ class Indexable(Nameable, Updateable):
 
     def log_prior(self):
         """evaluate the prior"""
-        if self.priors.size > 0:
-            x = self.param_array
-            #py3 fix
-            #return reduce(lambda a, b: a + b, (p.lnpdf(x[ind]).sum() for p, ind in self.priors.iteritems()), 0)
-            return reduce(lambda a, b: a + b, (p.lnpdf(x[ind]).sum() for p, ind in self.priors.items()), 0)
-        return 0.
+        if self.priors.size == 0:
+            return 0.
+        x = self.param_array
+        #evaluate the prior log densities
+        log_p = reduce(lambda a, b: a + b, (p.lnpdf(x[ind]).sum() for p, ind in self.priors.items()), 0)
+
+        #account for the transformation by evaluating the log Jacobian (where things are transformed)
+        log_j = 0.
+        priored_indexes = np.hstack([i for p, i in self.priors.items()])
+        for c,j in self.constraints.items():
+            if not isinstance(c, Transformation):continue
+            for jj in j:
+                if jj in priored_indexes:
+                    log_j += c.log_jacobian(x[jj])
+        return log_p + log_j
 
     def _log_prior_gradients(self):
         """evaluate the gradients of the priors"""
-        if self.priors.size > 0:
-            x = self.param_array
-            ret = np.zeros(x.size)
-            #py3 fix
-            #[np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.iteritems()]
-            [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.items()]
-            return ret
-        return 0.
+        if self.priors.size == 0:
+            return 0.
+        x = self.param_array
+        ret = np.zeros(x.size)
+        #compute derivate of prior density
+        [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.items()]
+        #add in jacobian derivatives if transformed
+        priored_indexes = np.hstack([i for p, i in self.priors.items()])
+        for c,j in self.constraints.items():
+            if not isinstance(c, Transformation):continue
+            for jj in j:
+                if jj in priored_indexes:
+                    ret[jj] += c.log_jacobian_grad(x[jj])
+        return ret
 
     #===========================================================================
     # Tie parameters together
diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py
index 7e15cee9..d53cb6c8 100644
--- a/GPy/core/parameterization/transformations.py
+++ b/GPy/core/parameterization/transformations.py
@@ -31,6 +31,16 @@ class Transformation(object):
         raise NotImplementedError
     def finv(self, model_param):
         raise NotImplementedError
+    def log_jacobian(self, model_param):
+        """
+        compute the log of the jacobian of f, evaluated at f(x)= model_param
+        """
+	raise NotImplementedError
+    def log_jacobian_grad(self, model_param):
+        """
+        compute the drivative of the log of the jacobian of f, evaluated at f(x)= model_param
+        """
+	raise NotImplementedError
     def gradfactor(self, model_param, dL_dmodel_param):
         """ df(opt_param)_dopt_param evaluated at self.f(opt_param)=model_param, times the gradient dL_dmodel_param,
 
@@ -74,9 +84,33 @@ class Logexp(Transformation):
         if np.any(f < 0.):
             print("Warning: changing parameters to satisfy constraints")
         return np.abs(f)
+    def log_jacobian(self, model_param):
+        return np.where(model_param>_lim_val, model_param, np.log(np.exp(model_param+1e-20) - 1.)) - model_param
+    def log_jacobian_grad(self, model_param):
+        return 1./(np.exp(model_param)-1.)
     def __str__(self):
         return '+ve'
 
+class Exponent(Transformation):
+    domain = _POSITIVE
+    def f(self, x):
+        return np.where(x<_lim_val, np.where(x>-_lim_val, np.exp(x), np.exp(-_lim_val)), np.exp(_lim_val))
+    def finv(self, x):
+        return np.log(x)
+    def gradfactor(self, f, df):
+        return np.einsum('i,i->i', df, f)
+    def initialize(self, f):
+        if np.any(f < 0.):
+            print("Warning: changing parameters to satisfy constraints")
+        return np.abs(f)
+    def log_jacobian(self, model_param):
+        return np.log(model_param)
+    def log_jacobian_grad(self, model_param):
+        return 1./model_param
+    def __str__(self):
+        return '+ve'
+
+
 
 class NormalTheta(Transformation):
     "Do not use, not officially supported!"
@@ -417,22 +451,6 @@ class LogexpClipped(Logexp):
     def __str__(self):
         return '+ve_c'
 
-class Exponent(Transformation):
-    # TODO: can't allow this to go to zero, need to set a lower bound. Similar with negative Exponent below. See old MATLAB code.
-    domain = _POSITIVE
-    def f(self, x):
-        return np.where(x<_lim_val, np.where(x>-_lim_val, np.exp(x), np.exp(-_lim_val)), np.exp(_lim_val))
-    def finv(self, x):
-        return np.log(x)
-    def gradfactor(self, f, df):
-        return np.einsum('i,i->i', df, f)
-    def initialize(self, f):
-        if np.any(f < 0.):
-            print("Warning: changing parameters to satisfy constraints")
-        return np.abs(f)
-    def __str__(self):
-        return '+ve'
-
 class NegativeExponent(Exponent):
     domain = _NEGATIVE
     def f(self, x):
diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py
index 7a58ac5e..f33e7715 100644
--- a/GPy/core/parameterization/variational.py
+++ b/GPy/core/parameterization/variational.py
@@ -176,11 +176,11 @@ class SpikeAndSlabPosterior(VariationalPosterior):
             self.mean.fix(warning=False)
             self.variance.fix(warning=False)
         if group_spike:
-            self.gamma_group = Param("binary_prob_group",binary_prob.mean(axis=0),Logistic(1e-6,1.-1e-6))
+            self.gamma_group = Param("binary_prob_group",binary_prob.mean(axis=0),Logistic(1e-10,1.-1e-10))
             self.gamma = Param("binary_prob",binary_prob, __fixed__)
             self.link_parameters(self.gamma_group,self.gamma)
         else:
-            self.gamma = Param("binary_prob",binary_prob,Logistic(1e-6,1.-1e-6))
+            self.gamma = Param("binary_prob",binary_prob,Logistic(1e-10,1.-1e-10))
             self.link_parameter(self.gamma)
             
     def propogate_val(self):
diff --git a/GPy/plotting/matplot_dep/visualize.py b/GPy/plotting/matplot_dep/visualize.py
index b57b7dfc..6903c2d2 100644
--- a/GPy/plotting/matplot_dep/visualize.py
+++ b/GPy/plotting/matplot_dep/visualize.py
@@ -408,12 +408,13 @@ class mocap_data_show_vpython(vpython_show):
 class mocap_data_show(matplotlib_show):
     """Base class for visualizing motion capture data."""
 
-    def __init__(self, vals, axes=None, connect=None):
+    def __init__(self, vals, axes=None, connect=None, color='b'):
         if axes==None:
             fig = plt.figure()
             axes = fig.add_subplot(111, projection='3d', aspect='equal')
         matplotlib_show.__init__(self, vals, axes)
 
+        self.color = color
         self.connect = connect
         self.process_values()
         self.initialize_axes()
@@ -423,7 +424,7 @@ class mocap_data_show(matplotlib_show):
         self.axes.figure.canvas.draw()
 
     def draw_vertices(self):
-        self.points_handle = self.axes.scatter(self.vals[:, 0], self.vals[:, 1], self.vals[:, 2])
+        self.points_handle = self.axes.scatter(self.vals[:, 0], self.vals[:, 1], self.vals[:, 2], color=self.color)
 
     def draw_edges(self):
         self.line_handle = []
@@ -442,7 +443,7 @@ class mocap_data_show(matplotlib_show):
                 z.append(self.vals[i, 2])
                 z.append(self.vals[j, 2])
                 z.append(np.NaN)
-            self.line_handle = self.axes.plot(np.array(x), np.array(y), np.array(z), 'b-')
+            self.line_handle = self.axes.plot(np.array(x), np.array(y), np.array(z), '-', color=self.color)
 
     def modify(self, vals):
         self.vals = vals.copy()
@@ -494,7 +495,7 @@ class stick_show(mocap_data_show):
 
 class skeleton_show(mocap_data_show):
     """data_show class for visualizing motion capture data encoded as a skeleton with angles."""
-    def __init__(self, vals, skel, axes=None, padding=0):
+    def __init__(self, vals, skel, axes=None, padding=0, color='b'):
         """data_show class for visualizing motion capture data encoded as a skeleton with angles.
         :param vals: set of modeled angles to use for printing in the axis when it's first created.
         :type vals: np.array
@@ -506,7 +507,7 @@ class skeleton_show(mocap_data_show):
         self.skel = skel
         self.padding = padding
         connect = skel.connection_matrix()
-        mocap_data_show.__init__(self, vals, axes=axes, connect=connect)
+        mocap_data_show.__init__(self, vals, axes=axes, connect=connect, color=color)
     def process_values(self):
         """Takes a set of angles and converts them to the x,y,z coordinates in the internal prepresentation of the class, ready for plotting.