Negative binomial fit over millions of data

Thanks, I’ve considered that. Ultimately I want to build a framework to generate predictive distributions, but right now I’m testing the feasibility of using Pymc3.

So some background, Y is purchase count, X are some features. Thinking of GLM with NB distribution, I want to predict Y_hat for each person. After the model is built, I want to identify the people with most purchase, but with some uncertainty. Although person A is better than B in terms of MLE, but I want to know the probability for later usage.