Initial evaluation of model at starting point failed, when attempting to use truncated model

Questions version agnostic
1

I am new to PyMC, and I am attempting to estimate the parameters of a Generalized Extreme Value distribution for my dataset. My data has been truncated artificially, and so I would like to use a truncated version of this distribution. I am able to estimate parameters for the GEV distribution for my data without truncation, but I am having trouble making the distribution truncated.

The following DOES work, without truncation:

data = read_in_data() # a numpy array of floats
basic_model = pm.Model()
with basic_model: 
    # Priors 
    mu = pm.Uniform("mu", lower = 1300, upper = 1500)
    sigma = pm.Uniform("sigma", lower = 90, upper = 150)
    xi = pm.TruncatedNormal("chi", mu=0, sigma=0.4, lower=-0.7, upper=0.7)   
    # Estimation
    gev = pmx.GenExtreme("gev", mu=mu, sigma=sigma, xi=xi, observed=data)
    gev_idata = pm.sample()

However, I believe my data comes from a truncated GEV distribution. My attempt to use a truncated GEV distribution, after reading the documentation here is as follows:

data = read_in_data() # a numpy array of floats
basic_model = pm.Model()
with basic_model: 
    # Priors 
    mu = pm.Uniform("mu", lower = 1300, upper = 1500)
    sigma = pm.Uniform("sigma", lower = 90, upper = 150)
    xi = pm.TruncatedNormal("xi", mu=0, sigma=0.4, lower=-0.7, upper=0.7)   
    # Estimation
    gev_full = pmx.GenExtreme.dist(mu=mu, sigma=sigma, xi=xi)
    gev = pm.Truncated("gev", gev_full, upper=2000, lower = None, observed=z_ions)
    gev_idata = pm.sample()

I get the error:

SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'mu_interval__': array(-0.53549444), 'sigma_interval__': array(-0.41259799), 'xi_interval__': array(0.48288129)}

Logp initial evaluation results:
{'mu': -1.46, 'sigma': -1.43, 'xi': -1.18, 'gev': -inf}
You can call `model.debug()` for more details.

I’ve got this error many times before, without using the truncated distribution, and never understood why it was trying to evaluate starting values of the parameters outside the distributions I specified. I see that it’s giving an infinite value of the gev distribution which is why it is erroring out, but I don’t understand why. In the past, just playing around with the priors for some reason fixed this error, but I can’t seem to find a set of priors that works in this case.

I assume I’m making a rookie mistake in implementing the truncated distribution.
How can I implement a truncated GEV distribution? Thank you so much in advance!

2

Welcome!

You can simply specify valid starting values for your model parameters. This can either be done on a per-parameter basis:

with pm.Model() as basic_model: 
    # Priors 
    mu = pm.Uniform("mu", lower = 1300, upper = 1500, initval=1400)
    sigma = pm.Uniform("sigma", lower = 90, upper = 150, initval=120)
    xi = pm.TruncatedNormal("xi", mu=0, sigma=0.4, lower=-0.7, upper=0.7, initval=0)
    # Estimation
    gev_full = pmx.GenExtreme.dist(mu=mu, sigma=sigma, xi=xi)
    gev = pm.Truncated("gev", gev_full, upper=2000, lower = None, observed=[1400])
    gev_idata = pm.sample(init="adapt_diag")

Or you can do it all in one go when calling pm.sample():

with basic_model:
    gev_idata = pm.sample(init="adapt_diag", initvals={"mu":1400, "sigma":120, "xi":0})

In both cases note I have changed the default initialization routine to init="adapt_diag". This is because the default (init="jitter+adapt_diag") will jitter the starting values and thus can still end up in invalid parts of the parameter space.

3

Thank you so much for the help!
I attempted to do as you suggested, and replaced the sampling line with (using values that I know these parameters should be very close to):

gev_idata = pm.sample(init="adapt_diag", initvals={"mu":1440, "sigma":125, "xi":0.03})

However, I still get the same error:

SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'mu_interval__': array(0.84729786), 'sigma_interval__': array(0.), 'xi_interval__': array(-0.28185115)}

Logp initial evaluation results:
{'mu': -1.56, 'sigma': -1.39, 'xi': -1.16, 'gev': -inf}
You can call `model.debug()` for more details.

I don’t understand how it can be picking initial values that aren’t allowed given the priors. For instance, my prior on mu is a uniform distribution between 1300 and 1500, yet if I am interpreting the output correctly it is trying to use 0.84729786 as the initial value for mu.

4

The value you’re seeing is on the unconstrained transformed scale, note the __interval__ suffix

5

Ok, I see, thank you. I guess regardless of that though, it’s still giving me the SamplingError. Do you have any thoughts as to why?

6

“It” hasn’t picked any parameter values. It is (now) using the parameter values you have supplied. and the model parameter values are not the problem. The error message tells you the (log) probability of each model parameter and they are all fine. However, the probability of gev is -inf (i.e., impossible according to your model/parameter values). Because gev is an observed variable (rather than a model parameter), this suggests that one or more observations are inconsistent with your model and the starting parameter values. This could happen, for example, if any of your observations are >2000 (straying outside the interval specified for the truncation transform). But there could be other reasons as well.

7

Thank you so much! You’re right, the issue is that my dataset does have a few values past 2000. I didn’t realize that would give inf values for the truncated distribution. I still want the distribution they are theoretically drawn from to cap at 2000, but that’s an issue with my data not PyMC. Thanks again for your help in pointing this out!

8

This has nothing to do with the trucation transformation per se. The probability of any “impossible” observation is zero, so logp==ln(0)==-inf. This is most commonly seen when feeding parameters that try to move model parameters outside where the associated distributions have support. This is a bit of a contrived example, but it illustrates that you don’t need any explicit transformations:


with pm.Model() as model: 
    # hyper parameter will almost certainly be negative
    sigma = pm.Normal("sigma", mu=-1000, sigma=1)
    # half normal has support on R+ and the observed data
    # is consistent with this support
    # but the negative sigma makes the positive observation "impossible"
    obs = pm.HalfNormal("obs", sigma=sigma, observed=[1])
    pm.sample()