I use the pm.HurdleGamma (predominantly in BAMBI) often in my work. The native pm.HurdleGamma distribution has some issues . Jax won’t work. Nutpie will work but is frightfully slow because it drops into object mode. For your 2,000 observations, see how the sample with the following argument in pm.sample: compile_kwargs={‘mode’:‘FAST_RUN’}. How does that compare to your 2.5 min run? Millions of records? Even if it ran ok would be a challenge scaling up to your target data. I would recommend ADVI --, especially ADVI mini-batch. Another option would be to fire up PYMC Experimental and use pm.fit(method=‘pathfinder’) instead of pm.sample.
Finally, you could just do a custom HurdleGamma. I have done that and gotten good results. That also then opens up JAX, Blackjax, Nutpie and Bayeux and its samplers to you.
Do you have some example data that I could try a custom model with?