diff --git a/GPy/util/univariate_Gaussian.py b/GPy/util/univariate_Gaussian.py
index 702ab25c..b5472e0a 100644
--- a/GPy/util/univariate_Gaussian.py
+++ b/GPy/util/univariate_Gaussian.py
@@ -40,6 +40,37 @@ def std_norm_cdf(x):
     weave.inline(code, arg_names=['x', 'cdf_x', 'N'], support_code=support_code)
     return cdf_x
 
+def std_norm_cdf_np(x):
+    """
+    Cumulative standard Gaussian distribution
+    Based on Abramowitz, M. and Stegun, I. (1970)
+    Around 3 times slower when x is a scalar otherwise quite a lot slower
+    """
+    x_shape = np.asarray(x).shape
+
+    if len(x_shape) == 0 or x_shape[0] == 1:
+        sign = np.sign(x)
+        x *= sign
+        x /= np.sqrt(2.)
+        t = 1.0/(1.0 +  0.3275911*x)
+        erf = 1. - np.exp(-x**2)*t*(0.254829592 + t*(-0.284496736 + t*(1.421413741 + t*(-1.453152027 + t*(1.061405429)))))
+        cdf_x = 0.5*(1.0 + sign*erf)
+        return cdf_x
+    else:
+        x = np.atleast_1d(x).copy()
+        cdf_x = np.zeros_like(x)
+        sign = np.ones_like(x)
+        neg_x_ind = x<0
+        sign[neg_x_ind] = -1.0
+        x[neg_x_ind] = -x[neg_x_ind]
+        x /= np.sqrt(2.)
+        t = 1.0/(1.0 +  0.3275911*x)
+        erf = 1. - np.exp(-x**2)*t*(0.254829592 + t*(-0.284496736 + t*(1.421413741 + t*(-1.453152027 + t*(1.061405429)))))
+        cdf_x = 0.5*(1.0 + sign*erf)
+        cdf_x = cdf_x.reshape(x_shape)
+    return cdf_x
+
+
 def inv_std_norm_cdf(x):
     """
     Inverse cumulative standard Gaussian distribution