Default NUTS initialization fails on GLM

[kyleabeauchamp]

Description of your problem

I'm building a simple normal GLM with (roughly) continuous y and categorical features. I find that the default NUTS initializer (jitter+adapt_diag) fails (exception below), while init="advi" runs fine. While this is not a dealbreaker, it would be great if the defaults were robust to this sort of thing. In particular, I was hoping to hand off my work to some co-workers who are not pymc developers, but I think it's too much for me to ask if the modeling requires intervening on the sampler.

with pm.Model() as model:
    pm.glm.GLM.from_formula("y ~ x", x)
    trace = pm.sample(3000, chains=1)

Please provide the full traceback.

lib/python3.6/site-packages/pymc3/step_methods/hmc/quadpotential.py in raise_ok(self)
    188     def raise_ok(self):
    189         if np.any(self._stds == 0):
--> 190             raise ValueError('Mass matrix contains zeros on the diagonal. '
    191                              'Some derivatives might always be zero.')
    192         if np.any(self._stds < 0):

ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero.

Please provide any additional information below.

Versions and main components

• PyMC3 Version: 3.3
• Theano Version: 1.0.1
• Python Version: 3.6
• Operating system: OSX
• How did you install PyMC3: (conda/pip) conda

[junpenglao]

Thanks for reporting. Yes, this issue has been noticed a couple of times (eg. https://discourse.pymc.io/t/limit-or-prevent-unrealistic-output-of-neural-network/871/4?u=junpenglao). The reason is that jitter sometimes makes the initial value goes out of the support of the logp.

We will make this more robust after the PR #2876

[AlexAndorra]

I think this is fixed now on PyMC 3.9.1 -- closing