Step transformation added

GPy/likelihoods/noise_model_constructors.py

@@ -17,9 +17,16 @@ def binomial(gp_link=None):
 
     if isinstance(gp_link,noise_models.gp_transformations.Probit):
         analytical_mean = True
+        analytical_variance = False
+
+    elif isinstance(gp_link,noise_models.gp_transformations.Step):
+        analytical_mean = True
+        analytical_variance = True
+
     else:
         analytical_mean = False
-    analytical_variance = False
+        analytical_variance = False
+
     return noise_models.binomial_noise.Binomial(gp_link,analytical_mean,analytical_variance)
 
 def exponential(gp_link=None):

GPy/likelihoods/noise_models/binomial_noise.py

@@ -42,11 +42,18 @@ class Binomial(NoiseDistribution):
         :param tau_i: precision of the cavity distribution (float)
         :param v_i: mean/variance of the cavity distribution (float)
         """
-        z = data_i*v_i/np.sqrt(tau_i**2 + tau_i)
-        Z_hat = std_norm_cdf(z)
-        phi = std_norm_pdf(z)
-        mu_hat = v_i/tau_i + data_i*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
-        sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
+        if isinstance(self.gp_link,gp_transformations.Probit):
+            z = data_i*v_i/np.sqrt(tau_i**2 + tau_i)
+            Z_hat = std_norm_cdf(z)
+            phi = std_norm_pdf(z)
+            mu_hat = v_i/tau_i + data_i*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
+            sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
+
+        elif isinstance(self.gp_link,gp_transformations.Step):
+            Z_hat = None
+            mu_hat = None
+            sigma2_hat = None
+
         return Z_hat, mu_hat, sigma2_hat
 
     def _predictive_mean_analytical(self,mu,sigma):

GPy/likelihoods/noise_models/gp_transformations.py

@@ -55,6 +55,8 @@ class Probit(GPTransformation):
     def d2transf_df2(self,f):
         return -f * std_norm_pdf(f)
 
+    def 
+
 class Log(GPTransformation):
     """
     $$
@@ -108,4 +110,18 @@ class Reciprocal(GPTransformation):
     def d2transf_df2(self,f):
         return 2./f**3
 
+class Step(GPTransformation):
+    """
+    $$
+    g(f) = I_{x \in A}
+    $$
+    """
+    def transf(self,f):
+        #transformation goes here
+        return np.where(f>0, 1, 0)
 
+    def dtransf_df(self,f):
+        pass
+
+    def d2transf_df2(self,f):
+        pass