Fix ExGaussian logp (#4050)

* Fix exgaussian logp

* Updated release notes

Author: AlexAndorra

Diff for: RELEASE-NOTES.md

@@ -4,7 +4,7 @@
 
 ### Maintenance
 - Mentioned the way to do any random walk with `theano.tensor.cumsum()` in `GaussianRandomWalk` docstrings (see [#4048](https://github.com/pymc-devs/pymc3/pull/4048)).
-
+- Fixed numerical instability in ExGaussian's logp by preventing `logpow` from returning `-inf` (see [#4050](https://github.com/pymc-devs/pymc3/pull/4050)).
 
 ### Documentation
 

Diff for: pymc3/distributions/continuous.py

@@ -18,6 +18,7 @@
 nodes in PyMC.
 """
 import numpy as np
+import theano
 import theano.tensor as tt
 from scipy import stats
 from scipy.special import expit
@@ -3268,13 +3269,21 @@ def logp(self, value):
         sigma = self.sigma
         nu = self.nu
 
-        # This condition suggested by exGAUS.R from gamlss
-        lp = tt.switch(tt.gt(nu,  0.05 * sigma),
-                       - tt.log(nu) + (mu - value) / nu + 0.5 * (sigma / nu)**2
-                       + logpow(std_cdf((value - mu) / sigma - sigma / nu), 1.),
-                       - tt.log(sigma * tt.sqrt(2 * np.pi))
-                       - 0.5 * ((value - mu) / sigma)**2)
-        return bound(lp, sigma > 0., nu > 0.)
+        standardized_val = (value - mu) / sigma
+        cdf_val = std_cdf(standardized_val - sigma / nu)
+        cdf_val_safe = tt.switch(tt.eq(cdf_val, 0), np.finfo(theano.config.floatX).eps, cdf_val)
+
+        # This condition is suggested by exGAUS.R from gamlss
+        lp = tt.switch(
+            tt.gt(nu, 0.05 * sigma),
+            -tt.log(nu)
+            + (mu - value) / nu
+            + 0.5 * (sigma / nu) ** 2
+            + logpow(cdf_val_safe, 1.0),
+            -tt.log(sigma * tt.sqrt(2 * np.pi)) - 0.5 * standardized_val ** 2,
+        )
+
+        return bound(lp, sigma > 0.0, nu > 0.0)
 
     def _repr_latex_(self, name=None, dist=None):
         if dist is None: