From 68f493b86cd152096852548f9ef891db95bea4e0 Mon Sep 17 00:00:00 2001
From: Ricardo <acq11ra@sheffield.ac.uk>
Date: Mon, 13 May 2013 11:10:27 +0100
Subject: [PATCH] New functions for EP-matching moments

---
 GPy/likelihoods/likelihood_functions.py |  7 ++---
 GPy/util/univariate_Gaussian.py         | 35 +++++++++++++++++++++++++
 2 files changed, 39 insertions(+), 3 deletions(-)
 create mode 100644 GPy/util/univariate_Gaussian.py

diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py
index 1196d88d..337c40fd 100644
--- a/GPy/likelihoods/likelihood_functions.py
+++ b/GPy/likelihoods/likelihood_functions.py
@@ -7,6 +7,7 @@ from scipy import stats
 import scipy as sp
 import pylab as pb
 from ..util.plot import gpplot
+from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
 
 class likelihood_function:
     """
@@ -37,11 +38,11 @@ class probit(likelihood_function):
         :param tau_i: precision of the cavity distribution (float)
         :param v_i: mean/variance of the cavity distribution (float)
         """
-        if data_i == 0: data_i = -1 #NOTE Binary classification algorithm works better with classes {-1,1}, 1D-plotting works better with classes {0,1}.
+        #if data_i == 0: data_i = -1 #NOTE Binary classification algorithm works better with classes {-1,1}, 1D-plotting works better with classes {0,1}.
         # TODO: some version of assert
         z = data_i*v_i/np.sqrt(tau_i**2 + tau_i)
-        Z_hat = stats.norm.cdf(z)
-        phi = stats.norm.pdf(z)
+        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)
         return Z_hat, mu_hat, sigma2_hat
diff --git a/GPy/util/univariate_Gaussian.py b/GPy/util/univariate_Gaussian.py
new file mode 100644
index 00000000..28946894
--- /dev/null
+++ b/GPy/util/univariate_Gaussian.py
@@ -0,0 +1,35 @@
+# Copyright (c) 2012, 2013 Ricardo Andrade
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from scipy import weave
+
+def std_norm_pdf(x):
+    """Standard Gaussian density function"""
+    return 1./np.sqrt(2.*np.pi)*np.exp(-.5*x**2)
+
+def std_norm_cdf(x):
+    """
+    Cumulative standard Gaussian distribution
+    Based on Abramowitz, M. and Stegun, I. (1970)
+    """
+    support_code = "#include <math.h>"
+    code = """
+
+    double sign = 1.0;
+    if (x < 0.0){
+        sign = -1.0;
+        x = -x;
+    }
+    x = x/sqrt(2.0);
+
+    double t = 1.0/(1.0 +  0.3275911*x);
+
+    double erf = 1. - exp(-x*x)*t*(0.254829592 + t*(-0.284496736 + t*(1.421413741 + t*(-1.453152027 + t*(1.061405429)))));
+
+    return_val = 0.5*(1.0 + sign*erf);
+    """
+    x = float(x)
+    return weave.inline(code,arg_names=['x'],support_code=support_code)
+
+