Normal CDF and test (does not yet pass)

Author: fonnesbeck

pymc3/distributions/continuous.py

@@ -15,7 +15,8 @@
 from pymc3.theanof import floatX
 from . import transforms
 
-from .dist_math import bound, logpow, gammaln, betaln, std_cdf, i0, i1, alltrue_elemwise
+from .dist_math import (bound, logpow, gammaln, betaln, std_cdf, i0, i1, 
+                        alltrue_elemwise, zvalue)
 from .distribution import Continuous, draw_values, generate_samples, Bound
 
 __all__ = ['Uniform', 'Flat', 'Normal', 'Beta', 'Exponential', 'Laplace',
@@ -231,6 +232,18 @@ def logp(self, value):
         return bound((-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2.,
                      sd > 0)
 
+    def logcdf(self, value):
+        mu = self.mu
+        sd = self.sd
+        z = zvalue(value, mu=mu, sd=sd)
+
+        return tt.switch(
+            tt.lt(z, -1.0),
+            tt.log(tt.erfcx(-z / tt.sqrt(2.)) / 2.) -
+            tt.sqr(tt.abs_(z)) / 2,
+            tt.log1p(-tt.erfc(z / tt.sqrt(2.)) / 2.)
+        )
+        
 
 class HalfNormal(PositiveContinuous):
     R"""

pymc3/distributions/dist_math.py

@@ -96,3 +96,9 @@ def i1(x):
                      x**9 / 1474560 + x**11 / 176947200 + x**13 / 29727129600,
                      np.e**x / (2 * np.pi * x)**0.5 * (1 - 3 / (8 * x) + 15 / (128 * x**2) + 315 / (3072 * x**3)
                                                        + 14175 / (98304 * x**4)))
+
+def zvalue(value, sd=1, mu=0):
+    """
+    Calculate the z-value for a normal distribution. By default standard normal.
+    """
+    return (value - mu) / tt.sqrt(2 * sd**2)

pymc3/tests/test_distributions.py

@@ -296,8 +296,8 @@ def pymc3_matches_scipy(self, pymc3_dist, domain, paramdomains, scipy_dist, extr
         model = build_model(pymc3_dist, domain, paramdomains, extra_args)
         value = model.named_vars['value']
 
-        def logp(args):
-            return scipy_dist(**args)
+        def logp(kwargs):
+            return scipy_dist(**kwargs)
         self.check_logp(model, value, domain, paramdomains, logp)
 
     def check_logp(self, model, value, domain, paramdomains, logp_reference):
@@ -307,6 +307,17 @@ def check_logp(self, model, value, domain, paramdomains, logp_reference):
         for pt in product(domains, n_samples=100):
             pt = Point(pt, model=model)
             assert_almost_equal(logp(pt), logp_reference(pt), decimal=6, err_msg=str(pt))
+            
+    def check_logcdf(self, pymc3_dist, domain, paramdomains, scipy_logcdf):
+        domains = paramdomains.copy()
+        domains['value'] = domain
+        for pt in product(domains, n_samples=100):
+            params = dict(pt)
+            scipy_cdf = scipy_logcdf(**params)
+            value = params.pop('value')
+            dist = pymc3_dist.dist(**params)
+            assert_almost_equal(dist.logcdf(value).tag.test_value, scipy_cdf, 
+                                decimal=6, err_msg=str(pt))
 
     def check_int_to_1(self, model, value, domain, paramdomains):
         pdf = model.fastfn(exp(model.logpt))
@@ -379,6 +390,8 @@ def test_flat(self):
     def test_normal(self):
         self.pymc3_matches_scipy(Normal, R, {'mu': R, 'sd': Rplus},
                                  lambda value, mu, sd: sp.norm.logpdf(value, mu, sd))
+        self.check_logcdf(Normal, R, {'mu': R, 'sd': Rplus}, 
+                                lambda value, mu, sd: sp.norm.logcdf(value, mu, sd))
 
     def test_half_normal(self):
         self.pymc3_matches_scipy(HalfNormal, Rplus, {'sd': Rplus},