# Fitting a distribution to my timeseries: two R-packages, two contrary results

typical R-User here, applying a bunch of packages to my data, hoping for a convincing result although I understand only half of what I do at most (and then getting rejected by the reviewers, surprise). But bear with me, I am here to learn.

I have a time-series of monthly precipitation data, for which I'd like to find a suitable distribution. The Cullen and Frey graph from the fitdistrplus package indicates that a lognormal distribution should outperform a gamma distribution:

and the results of fitting a gamma, Weibull and lognormal distribution seem to say the same:

Goodness-of-fit statistics
Weibull  lognormal      gamma
Kolmogorov-Smirnov statistic 0.0263855 0.03245692 0.07115254
Cramer-von Mises statistic   0.1427163 0.11374430 1.13149886
Anderson-Darling statistic   1.0934641 0.62867455 6.73258465

Goodness-of-fit criteria
Weibull lognormal    gamma
Akaike's Information Criterion 7845.461  7842.426 7954.480
Bayesian Information Criterion 7854.940  7851.905 7963.959


This is of particular interest to me, as in my field of research a gamma distribution is often used for precipitation data, without considering that there might be better fitting distributions.

I was quite satisfied with my results until I stumbled upon the GAMLSS package, which offers an all-in-one solution for finding a suiting distribution to your data. So I applied the easy to use fitDist' function to see what's the outcome:

> fit <- fitDist(df$$prcp, k = round(log(length(df$$prcp))), type = "realplus", trace = F , try.gamlss = TRUE)
There were 11 warnings (use warnings() to see them)
> summary(fit)
*******************************************************************
Family:  c("BCCG", "Box-Cox-Cole-Green")

Call:  gamlssML(formula = y, family = DIST[i], data = sys.parent())

Fitting method: "nlminb"

Coefficient(s):
Estimate  Std. Error   t value   Pr(>|t|)
eta.mu    45.1741563   0.9560784  47.24943 < 2.22e-16 ***
eta.sigma -0.5522153   0.0265535 -20.79636 < 2.22e-16 ***
eta.nu     0.4250894   0.0439835   9.66474 < 2.22e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Degrees of Freedom for the fit: 3 Residual Deg. of Freedom   842
Global Deviance:     7827.23
AIC:     7833.23
SBC:     7847.45


First of all, the function offers me a distribution I've never heard of, but that's fair since (you might've already guessed it) I am quite new to this statistical topic it makes sense that I don't know "all" distributions. But, as I investigated the results of the other distributions, it got quite interesting (values are the BIC):

> fit\$fits
BCCG    BCCGo       GA       GG      BCT     BCTo     BCPE    BCPEo     WEI2      WEI     WEI3      GB2      GIG   exGAUS      LNO
7848.230 7848.230 7852.426 7853.005 7854.763 7854.763 7854.989 7854.989 7855.461 7855.461 7855.461 7857.338 7859.426 7871.134 7964.480
LOGNO2    LOGNO       IG      EXP PARETO2o  PARETO2       GP   IGAMMA
7964.480 7964.480 8260.617 8290.105 8297.106 8297.106 8297.112 8521.495


Right behind the mysterious BCCG-distributions follows the GA (gamma) distribution, quite closely followed by Weibull, and then, lognormal ranked even lower.

So here are my questions:

1. What's the reason for these contrary results? Is it due to different fitting methods?
2. What can I further do to find the most fitting distribution?

I know in general that all models are wrong and that it probably won't make a huge difference whether I pick the result of fitdistrplus or GAMLSS` for my further analysis (calculation of SPI). But still, I am wondering about the comparability of research:

1. In a publication, I would probably state the package that I used for the fitting and present some goodness-of-fit values (such as the BIC e.g.). But, obviously, there's more to that and the results depend on the package used. But in the publications that I investigated, there's rarely more information offered. Doesn't that hinder the comparability of research?
• At a guess, you have one large outlier. (I'm having to make some suppositions about what your graphics and output represent.) It's that data value exceeding 200. What results do you get when omitting that value? I'm also going to guess that much of the output will be qualitatively the same but that the new Cullen-Frey graph may give you a very different impression.
– whuber
Jun 20, 2020 at 20:07
• @whuber you are right and I thought about removing it before, but I have cross-checked this measurement and I can assume that the outlier is true. Nevertheless, I tested it without the outlier and you were right, the Cullen-Frey graph changes dramatically, but doesn't have much impact on the actual results. Jun 22, 2020 at 7:46