diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py
index 44b48ccb..a209cb39 100644
--- a/GPy/core/parameterization/variational.py
+++ b/GPy/core/parameterization/variational.py
@@ -37,7 +37,7 @@ class SpikeAndSlabPrior(VariationalPrior):
     def __init__(self, pi, variance = 1.0, name='SpikeAndSlabPrior', **kw):
         super(VariationalPrior, self).__init__(name=name, **kw)
         assert variance==1.0, "Not Implemented!"
-        self.pi = Param('pi', pi, Logistic(1e-10,1-1e-10))
+        self.pi = Param('pi', pi, Logistic(1e-10,1.-1e-10))
         self.variance = Param('variance',variance)
         self.add_parameters(self.pi)
         
@@ -105,7 +105,7 @@ class SpikeAndSlabPosterior(VariationalPosterior):
         binary_prob : the probability of the distribution on the slab part.
         """
         super(SpikeAndSlabPosterior, self).__init__(means, variances, name)
-        self.gamma = Param("binary_prob",binary_prob, Logistic(1e-10,1-1e-10))
+        self.gamma = Param("binary_prob",binary_prob, Logistic(1e-10,1.-1e-10))
         self.add_parameter(self.gamma)
 
     def plot(self, *args):
diff --git a/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py b/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py
index 13b3f65f..d8414cfb 100644
--- a/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py
+++ b/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py
@@ -42,7 +42,8 @@ def _psi1computations(variance, lengthscale, Z, mu, S, gamma):
     _psi1_common = gamma[:,None,:] / (lengthscale2*_psi1_denom*_psi1_denom_sqrt) #Nx1xQ
     _psi1_exponent1 = np.log(gamma[:,None,:]) -0.5 * (_psi1_dist_sq + np.log(_psi1_denom)) # NxMxQ
     _psi1_exponent2 = np.log(1.-gamma[:,None,:]) -0.5 * (np.square(Z[None,:,:])/lengthscale2) # NxMxQ
-    _psi1_exponent = np.log(np.exp(_psi1_exponent1) + np.exp(_psi1_exponent2)) #NxMxQ
+    _psi1_exponent_max = np.maximum(_psi1_exponent1,_psi1_exponent2)
+    _psi1_exponent = _psi1_exponent_max+np.log(np.exp(_psi1_exponent1-_psi1_exponent_max) + np.exp(_psi1_exponent2-_psi1_exponent_max)) #NxMxQ
     _psi1_exp_sum = _psi1_exponent.sum(axis=-1) #NxM
     _psi1_exp_dist_sq = np.exp(-0.5*_psi1_dist_sq) # NxMxQ
     _psi1_exp_Z = np.exp(-0.5*np.square(Z[None,:,:])/lengthscale2) # 1xMxQ
@@ -89,7 +90,8 @@ def _psi2computations(variance, lengthscale, Z, mu, S, gamma):
     _psi2_common = gamma[:,None,None,:]/(lengthscale2 * _psi2_denom * _psi2_denom_sqrt) # Nx1x1xQ
     _psi2_exponent1 = -_psi2_Zdist_sq -_psi2_mudist_sq -0.5*np.log(_psi2_denom)+np.log(gamma[:,None,None,:]) #N,M,M,Q
     _psi2_exponent2 = np.log(1.-gamma[:,None,None,:]) - 0.5*(_psi2_Z_sq_sum) # NxMxMxQ
-    _psi2_exponent = np.log(np.exp(_psi2_exponent1) + np.exp(_psi2_exponent2))
+    _psi2_exponent_max = np.maximum(_psi2_exponent1, _psi2_exponent2)
+    _psi2_exponent = _psi2_exponent_max+np.log(np.exp(_psi2_exponent1-_psi2_exponent_max) + np.exp(_psi2_exponent2-_psi2_exponent_max))
     _psi2_exp_sum = _psi2_exponent.sum(axis=-1) #NxM
     _psi2_q = np.square(variance) * np.exp(_psi2_exp_sum[:,:,:,None]-_psi2_exponent) # NxMxMxQ 
     _psi2_exp_dist_sq = np.exp(-_psi2_Zdist_sq -_psi2_mudist_sq) # NxMxMxQ
diff --git a/GPy/models/ss_gplvm.py b/GPy/models/ss_gplvm.py
index 8c9d6f23..37309c94 100644
--- a/GPy/models/ss_gplvm.py
+++ b/GPy/models/ss_gplvm.py
@@ -36,7 +36,7 @@ class SSGPLVM(SparseGP):
             X_variance = np.random.uniform(0,.1,X.shape)
             
         gamma = np.empty_like(X) # The posterior probabilities of the binary variable in the variational approximation
-        gamma[:] = 0.5
+        gamma[:] = 0.5 + 0.01 * np.random.randn(X.shape[0], input_dim)
         
         if Z is None:
             Z = np.random.permutation(X.copy())[:num_inducing]