Fit distributions with glm

I'm trying to fit different statistical distributions (Gamma, Poisson, normal, inverse Gaussian) to my data with a glm. An example could be like this:

data   <- rexp(500, 3)
model  <- glm(data~1, family=Gamma())
shape  <- 1 / gamma.dispersion(model)
rate   <- shape*model$coef[1] model2 <- glm(data~1, family=poisson()) lambda <- unique(model2$fitted.values)


So, here are the questions:

1. If I want to select between all of them the best distribution fit with a glm, can I use the models AIC? Or should I use another method, like MSE? It's because I have really different results from one distribution to another.
2. I already know that packages fitdistrplus or MASS include specific functions to do this, but it's too slow for my data. If I use those functions to fit the data, mostly of the time I get as a result an exponential distribution. Is there a way to fit it with the glm families?

Update: The data comes from sales orders, but it is always grater than 0, that's why I can use the exponential or gamma distributions. (Only include the normal for the cases where there are so many orders that the negative part under the curve is almost 0).

Also, I need the distribution for simulate from it, a new series. With the result showed by the AIC i can't make them compete to get the best fitting. I would like a formal test, but i also need it to be really fast, so I'll settle for an approximation or something exploratory.

• Where do these data come from? Why do you need to know what distribution they are (ie, what is your end goal)? What are the "really different results", & what is the impact of them? Do you need a formal test, or are you looking exploratatively for a better understanding? – gung - Reinstate Monica May 6 '14 at 15:24
• Thanks for the clarification. Why would it need to be really fast? Are you going to code this & have it run automatically at regular intervals? If you just need to figure out what the distribution is, presumably that could be done once & then you could simulate however much you want later. – gung - Reinstate Monica May 6 '14 at 16:21
• Yes, I will have to compute it for regular intervals. I have almost 500000 records, and for each one I have to compute it's distribution. With the other calculations of the code I reach the half a second for each record, so this makes little less than 3 days. I was thinking on keeping the distribution for each record, and this way save time, but I will have to compute the calculation for half of them periodically. – Tomás M. May 6 '14 at 16:40
• These data are "sales orders", are they counts (eg, 5 orders), or something like revenue (which can have fractional values)? – gung - Reinstate Monica May 6 '14 at 16:50
• Yes, they are the number of items ordered for each time. Why? – Tomás M. May 6 '14 at 16:56

Since your data are counts, they can only come in whole numbers. That is, you could get a count of $5$ orders or $6$ orders, but you could never get $5.5$ orders. Thus, your data cannot be distributed as Gamma, normal or inverse Gaussian, as these are continuous distributions. (It is always possible in some case that the fits will be good enough for your purposes, but your goal is to use the distribution to simulate new data, so we can just rule them out.)

For your case, the default first choice would be the Poisson distribution. However, the Poisson is very restrictive, in the sense that the variance must equal the mean. This can happen in the real world, but rarely does in practice and I highly doubt it would hold for sales. A distribution that relaxes this restriction / generalizes the Poisson is the negative binomial. This will allow you to estimate the mean and variance separately.

Thus, your first strategy for fitting these distributions and doing so quickly, is to only fit the negative binomial. The primary function for getting these parameters is ?fitdistr in the MASS package. Subsequently, you can use ?rnbinom for simulation.

If speed and size remain a problem, another possibility is to extract random subsamples of manageable size from your data and fit those. Even if these were a small proportion of the total, the fitted parameters should be very close to the true parameters if the sample is large enough.

• I think I can't use the subsamples, because the total data for each set is not really large. I will try only the negative binomial and check the results. Also, I tried the Poisson, that it was my first approach, but it fitted really bad the data. So my strategy was to add more distributions although they were for continuous data. Thanks a lot for the explanation. – Tomás M. May 12 '14 at 10:22

It is quite interesting to hear using glm to replace the fitdistr. But at least the code has some problems, which you may have found from the totally different AIC or likelihood.

The glm family distribution function will consider the link function by default. So in your code, gamma() means Gamma(link = "inverse") by default. And a log link for the poisson distribution. You may try gamma(link ="identity") and compare with the poisson(link = "identity")

• I already checked this, but I have one question. If I change the glm code like that, then, the parameters of the distributions have to be changed with the corresponding link function? Thanks! – Tomás M. May 6 '14 at 16:00
• @TomásM.Sorry for not understanding your exact meaning. But different link function will transform the underlying data first. So the result of fitting parameters are different and so are the AIC and llh. – Vincent May 6 '14 at 16:34
• Sorry for my bad english and not making me understanding. I was wondering if the parameters need to be changed in any way with the link function or they can be used to simulate a new dataset, similar to the fitted, as they are on the glm result. – Tomás M. May 6 '14 at 16:50
• @TomásM. Link function should be set as link ="identity". In this case, the parameter does not need transforming. Eg. lambda of Poisson, lambda from glm is the same as lambda from fitdistr and it can be used to simulate new distribution directly. But gamma distribution need to care about using the same pattern (shape/rate, or shape/scale). If the link function is not set to identity, some transformation has been done by glm by default, then the fitting parameters are not the same as fitdistr. This not only affects the value of parameter but also the standard error of the fitting. – Vincent May 7 '14 at 9:55