Bad initial energy when trying to build a multivariate ExGaussian model

Questions
1

I am trying to sample from a multivariate exmodified gaussian but somehow I keep getting bad initial energy error. I tried checking FAQ’s and also the workflow suggested by @ junpenglao however, I am not able to figure out if there some misspecification in my model.

with pm.Model() as exp_mod_mdel: 
    
    # Priors
    ν = pm.InverseGamma('ν', alpha=5)
    β = pm.Normal('β', 0, sd=10, shape=len(X.columns))
    σ = pm.InverseGamma('σ', alpha=3)
    
    #Liklihood mean
    μ = pm.Deterministic('μ', pm.math.dot(X, β))
    
    # Likelihood
    likelihood = pm.ExGaussian('y', mu=μ, sigma=σ, nu=ν, observed=y)

    # Inference!
    #trace = pm.sample(2000, tune=500, cores=2)

Running logp, dlogp = step.integrator._logp_dlogp_func(q0) tells me that logp is infinite, however when I check the test values for example using np.sum(~np.isfinite(μ.tag.test_value)) I do not find any non-finite values (I check this for all of my variables ).

I find that ‘y’ is infinite by running, exp_mod_mdel.check_test_point(), however, as I said using values of parameters, μ = 0, σ = 0.26, ν = 0.16 found using tag.test_value, gives me finite log pdf through:

import scipy.stats as st
x = np.linspace(-6, 9, 200)
log_pdf = np.log(st.exponnorm.pdf(x, ν.tag.test_value, loc=0, scale=σ.tag.test_value))
2

Hi,
Have you checked if there are missing values in your observed data, here y? This sounds like it could be the case

3

Hi,

There are no invalid values in y. I can also reproduce the same error on simulated data through following snippet:

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)

with pm.Model() as exp_mod_mdel: 
    
    # Priors
    ν = pm.InverseGamma('ν', alpha=5)
    β = pm.Normal('β', sd=10)
    σ = pm.InverseGamma('σ', alpha=3)
    
    #Liklihood mean
    μ = pm.Deterministic('μ', x*β)
    
    # Likelihood
    likelihood = pm.ExGaussian('y', mu=μ, sigma=σ, nu=ν, observed=y)

    # Inference!
    trace = pm.sample(500, tune=50, cores=2)

Interestingly, sometimes these work if I pass ‘advi_map’ as init when sampling, i.e. trace = pm.sample(500, tune=50, cores=2, init='advi_map') but only sometimes.

1 Like
4

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

I opened https://github.com/pymc-devs/pymc3/issues/4045, for now, try using larger sigma (e.g., a truncated Normal prior to avoid small sigma)

1 Like
5

@non_convex, the issue is now fixed on master :tada:

2 Likes
6

Thanks a lot @AlexAndorra @junpenglao