ExGaussian logp is numerical unstable

[junpenglao]

Reported by user on discourse: https://discourse.pymc.io/t/bad-initial-energy-when-trying-to-build-a-multivariate-exgaussian-model/5591/3

There is some numerical problem of the ExGaussian log_prob. Specifically, the std_cdf in the line below returns 0., which resulting in logpow returns -inf

https://github.com/pymc-devs/pymc3/blob/63eba59fe2c42c936136b24babf6ca1e3a777d74/pymc3/distributions/continuous.py#L3274

minimal reproducible example:

import numpy as np
import scipy.stats as st
import pymc3 as pm

beta = 0.5
sigma = 1.0
nu = 1.2

x = np.linspace(-2, 4, 200)
mu = beta*x
y = st.exponnorm.rvs(nu, loc=mu, scale=sigma)

pm.ExGaussian.dist(0., .25, 1./6).logp(y).tag.test_value  # <== contains -inf

A tt.switch in ExGaussian logp to replace 0 with epsilon should be sufficient to fix this.

[AlexAndorra]

Just tried this code to see if I can reproduce but it returns an error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-1-d4f49d90ca18> in <module>
     11 y = st.exponnorm.rvs(nu, loc=mu, scale=sigma)
     12 
---> 13 pm.ExGaussian.dist(0., .25, 1./6).logp(y).tag.test_value

AttributeError: 'scratchpad' object has no attribute 'test_value'

[junpenglao]

pm.ExGaussian.dist(0., .25, 1./6).logp(y).eval() should work

[AlexAndorra]

It does, thanks! And I can reproduce the problem:

array([ 1.09121345e-01, -2.75359986e+00, -1.01118924e+01,  2.39749482e-01,
       -3.03631310e-01,            -inf, -1.20340373e+00,  2.73649008e-01,
       -5.53724783e-01,  1.92089061e-01, -7.24266195e-01, -2.05114903e+00, ...]