From f2fa9bd74d0141c2869d2090ed28fa4f9ffa49e4 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Mon, 16 Sep 2013 16:20:26 +0100
Subject: [PATCH] addedHeraviside functionality to EP

---
 GPy/likelihoods/noise_model_constructors.py   |  2 +-
 .../noise_models/binomial_noise.py            | 27 +++++++++++++++----
 .../noise_models/gp_transformations.py        |  6 ++---
 3 files changed, 26 insertions(+), 9 deletions(-)

diff --git a/GPy/likelihoods/noise_model_constructors.py b/GPy/likelihoods/noise_model_constructors.py
index 158f1674..b3f4bdaf 100644
--- a/GPy/likelihoods/noise_model_constructors.py
+++ b/GPy/likelihoods/noise_model_constructors.py
@@ -19,7 +19,7 @@ def binomial(gp_link=None):
         analytical_mean = True
         analytical_variance = False
 
-    elif isinstance(gp_link,noise_models.gp_transformations.Step):
+    elif isinstance(gp_link,noise_models.gp_transformations.Heaviside):
         analytical_mean = True
         analytical_variance = True
 
diff --git a/GPy/likelihoods/noise_models/binomial_noise.py b/GPy/likelihoods/noise_models/binomial_noise.py
index 9ceacf79..abad4822 100644
--- a/GPy/likelihoods/noise_models/binomial_noise.py
+++ b/GPy/likelihoods/noise_models/binomial_noise.py
@@ -49,15 +49,32 @@ class Binomial(NoiseDistribution):
             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
+        elif isinstance(self.gp_link,gp_transformations.Heaviside):
+            a = data_i*v_i/np.sqrt(tau_i)
+            Z_hat = std_norm_cdf(a)
+            N = std_norm_pdf(a)
+            mu_hat = v_i/tau_i + data_i*N/Z_hat/np.sqrt(tau_i)
+            sigma2_hat = (1. - a*N/Z_hat - np.square(N/Z_hat))/tau_i
+            if np.any(np.isnan([Z_hat, mu_hat, sigma2_hat])):
+                stop
 
         return Z_hat, mu_hat, sigma2_hat
 
     def _predictive_mean_analytical(self,mu,sigma):
-        return stats.norm.cdf(mu/np.sqrt(1+sigma**2))
+        if isinstance(self.gp_link,gp_transformations.Probit):
+            return stats.norm.cdf(mu/np.sqrt(1+sigma**2))
+        elif isinstance(self.gp_link,gp_transformations.Heaviside):
+            return stats.norm.cdf(mu/sigma) 
+        else:
+            raise NotImplementedError
+
+    def _predictive_variance_analytical(self,mu,sigma, pred_mean):
+        if isinstance(self.gp_link,gp_transformations.Heaviside):
+            return 0.
+        else:
+            raise NotImplementedError
+
+
 
     def _mass(self,gp,obs):
         #NOTE obs must be in {0,1}
diff --git a/GPy/likelihoods/noise_models/gp_transformations.py b/GPy/likelihoods/noise_models/gp_transformations.py
index 5f337aa4..0f8cf18c 100644
--- a/GPy/likelihoods/noise_models/gp_transformations.py
+++ b/GPy/likelihoods/noise_models/gp_transformations.py
@@ -108,7 +108,7 @@ class Reciprocal(GPTransformation):
     def d2transf_df2(self,f):
         return 2./f**3
 
-class Step(GPTransformation):
+class Heaviside(GPTransformation):
     """
     $$
     g(f) = I_{x \in A}
@@ -119,7 +119,7 @@ class Step(GPTransformation):
         return np.where(f>0, 1, -1)
 
     def dtransf_df(self,f):
-        pass
+        raise NotImplementedError, "this function is not differentiable!"
 
     def d2transf_df2(self,f):
-        pass
+        raise NotImplementedError, "this function is not differentiable!"