1
$\begingroup$

I have some data on the duration of several activities (rounded to the nearest half hour). I'm trying to add up these random variables (one per activity) so that I can calculate the total duration of a project, as well as extract some summary statistics from it.

In order to do that, I'm trying to determine which distribution fits it reasonably well. This is for two reasons:

  • if I can model these RVs by using a "named" distribution, it becomes simpler to combine them
  • by fitting a "named" distibution, I assume I'll be able to infer the generator behind the data. For instance, if a lognormal is a good fit for the data, then it might be generated by some sort of random walk process.

I've selected five candidate distributions: lognorm, exponweib, norm, t and dweibull. Three of them are unbounded; I've chosen them just to see if they fit my data reasonably well, even though time durations can't possibly be negative.

I've also selected two criteria with which to judge goodness-of-fit, K-S and AIC. I wanted something that would apply to all of the distributions I selected, in an automated manner. I also wanted a criterion that would penalize more parameters.

However, something apparently strange happened. Here are the distributions' CDFs and the duration data empirical CDF (dashed line):

enter image description here

Judging by this plot, the norm, dweibull and t are all reasonably good fits, which is confimed by their K-S scores.

Now here are the same distributions in PDF form (as well as the data histogram):

enter image description here

Judging by this plot, the lognorm and exponweib are the clear winners, with dweibull a distant third. This is also confirmed by their AIC scores.

Here is the code in Python that calculates the AIC:

def aic(dist, dados, second_order = True):
    fit = dist.fit(dados)
    k = len(fit)
    lnL = dist(*fit).logpdf(dados).sum()
    aic = 2 * k - 2 * lnL
    if second_order:
        n = len(dados)
        aicc_p = 2 * ((k ** 2 + k) / (n - k - 1))
        aic += aicc_p

    return aic

In light of this, I have a few questions:

  1. Is this a good strategy (fitting a distribution to the data) for my problem (adding RVs up and extracting summary statistics)?
  2. Is it possible/surprising/normal that the K-S and AIC statistics give opposing results when used as goodness-of-fit criteria?
  3. Which distribution should I pick? Am I misusing these statistics? Is there some other consideration I should be making before choosing?
$\endgroup$
  • 2
    $\begingroup$ Why is it important to have a 'name' for the distribution of your data. // You really seen to have four methods of matching AIC, KS; ECDF, Histogram--your interpretation of plots does not match mine. // Weibull often 'wins' contests such as this because that family includes fundamentally different shapes. $\endgroup$ – BruceET Sep 29 at 19:03
  • 2
    $\begingroup$ This sounds like an XY problem. You want to solve problem X: how to generate new samples. To do that, you propose Y: fit a parametric distribution and sample from it. However, you are unsure how to do Y, so you ask about Y, even though what you really want to do is X. $\endgroup$ – Dave Oct 2 at 16:38
  • 1
    $\begingroup$ I still do not understand what your ultimate aim is. Do you want to (a) learn something about the data (b) be able to generate data similar to the observed data (c) find a function that somewhat resembles the empirical distribution $\endgroup$ – Aleksejs Fomins Oct 4 at 15:46
  • 1
    $\begingroup$ What do you mean by "K-S score"? The p-value of the K-S test? This is not a "score" and isn't meant to be a model selection criterion. Actually if it's p-values, the results suggest that all these distributions are compatible with the data, and can all be used to generate "similar" data. "My thinking is, if I can "match" a named distribution to my data, that is evidence that the generator for the distribution is the same (or close enough) as the generator behind my data." - Well, surely not in the sense that you could exclude any of the others, which you can't, according to K-S. $\endgroup$ – Lewian Oct 8 at 21:42
  • 1
    $\begingroup$ @LmniCE: To me it would make more sense to look at the p-value, because that will actually tell you if and which of these distributions are compatible with the data. It could be they all are, or it could be that none is, but it may also be that one or two are and the others not. To me this seems to address your question better than looking at the distance, although there may not be one clear winner. But that may just be the honest result here. $\endgroup$ – Lewian Oct 9 at 8:22
1
+50
$\begingroup$
  1. The AIC is normally a criterion that is used to compare models with different numbers of parameters to fit. It is not really clear exactly how many parameters you have fitted here (chances are 2-3 for most if not all of these, but for example it is not clear where you got the degrees of freedom of the t from; also it seems you have shifted or cut off distributions so that they start at 1/2 where the data start, rather than at zero (like, e.g., lognormal) or at $-\infty$ (like t or normal), so you used these in a nonstandard way and it may even debatable whether the shift to 1/2 should count as an additional parameter in the sense of AIC or not. In fact AIC differences are so striking that they cannot be explained in my opinion by different numbers of parameters. I wonder even whether the AIC was computed correctly (or whether the images reflect precisely what you did), because I haven't seen any such calculations for a truncated normal or t-distribution yet (no idea what dist.fit and len(fit) exactly do in your code, as I don't use python).

  2. In any case, even assuming that AIC values are correct, I wouldn't put much trust in them. This is because (as was stated already in the other answer) the data are highly discrete, and density (pdf) is a somewhat unstable concept in the sense that discretisation and changes of low amounts of probability can change densities a lot. The good results of lognorm and exponweib in terms of the density seem to be entirely due to the fact that they correctly capture the peak at 1/2; the quality of their fit at higher values (which is clearly worse than Weibull and somewhat worse than the other two) is in the density-based AIC-computation totally dominated by what goes on at 1/2. And this only works for those two distributions because you let them start at 1/2 rather than zero, which they theoretically do. (Of course this is not wrong in the sense that shifting them to 1/2 still leaves them as well defined distributions, but it would be hard to argue that this captures a real process). In fact it may be that the underlying more continuous real data that probably start at zero do not have a density peak at zero but somewhat higher, which would make a Weibull shape clearly more appropriate compared to lognormal or exponweib, but as the rounded data don't show anything below 1/2, the visible density peak is at 1/2, favouring in somewhat artificial manner the distributions that peak at their origin and then are shifted to 1/2.

  3. I recommend to look at K-S p-values, because these give you an indication if and whether these distributions are compatible with the data, which seems to be what you are interested in. This may rule out some (or all) of these distributions, or it may tell you that all are compatible, or it may tell you that some work and some don't, so it wouldn't necessarily pick a clear winner, but that will just be an honest result; it may be that the data cannot distinguish between all or some of these distributions apart from random variation. If you need to pick one that is best, you can still go for the one with lowest K-S distance. Note though that a proper computation of the K-S p-value needs to take into account the number of fitted parameters, and not all existing software does this. This of course brings in the difficulty of parameter counting mentioned for the AIC above (and in fact honest results here should depend on the number of degrees of freedom used in fitting), but I'd think even ignoring this, K-S p-values at least give some orientation, interpreted in a cautious "exploratory" way.

  4. An alternative approach would be to not shift or truncate distributions to 1/2, but rather to use them to generate continuous data, and then model applying the rounding/binning procedure as is done for the real data. This will give you probability mass functions starting at 1/2, which you can compare with the one from your real data in a more appropriate way than comparing the data pmf to continuous densities. (Also K-S distances can be computed in this way, and may then give a result more compatible to comparing based on densities/pmfs.)

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for the answer! dist.fit is a MLE fit for each distribution. len(fit) is number of fitted parameters. $\endgroup$ – LmnICE Oct 9 at 11:54
  • 1
    $\begingroup$ But what you show in your first plot is truncated normals and t-distributions, not just the standard form of normal and t. Did you run dist.fit with truncated distributions? Also I'd have guessed what len(fit) is, however this doesn't make exactly clear how this is counted - i.e., are the df of the t-distribution counted? Is the shift of the lognormal counted? Are the truncation parameters counted (in case you used them)? $\endgroup$ – Lewian Oct 9 at 13:40
  • $\begingroup$ I plot the fitted normal and t distributions in the data range. For instance If the data values range from 0.5 to 4.0, I plot the distributions using as x-axis values that go from 0.5 to 4.0. But neither of the distributions are truncated; it's just how they were plotted. $\endgroup$ – LmnICE Oct 9 at 18:29
  • $\begingroup$ As for the parameters: the df of the t is counted, as is the lognormal shift. In general, the fit vector contains all the parameters that are necessary to make a family of distributions unique. $\endgroup$ – LmnICE Oct 9 at 18:31
1
$\begingroup$

This seems to be an example of where binning a continuous variable is leading to problems, as it often does. With half-hour bins, over 3/4 of your observations are restricted to the first 3 bins. It's going to be hard to gauge any model fit very well, whether visually or by a statistic related to a fit, with data values restricted this way.

I'm also a bit worried that there seem to be no observed values in the time bin between 0 and 30 minutes in these histograms. If there were any actual times less than 15 minutes, rounding "to the nearest half hour" should lead to times of 0 for those cases.

And if you did round to the nearest half hour, it's not clear whether what you show in the bin from 1 to 1.5 hours, for example, represents actual task durations between 0.75 and 1.25 hours or between 1.25 and 1.75 hours. At least, for visual comparisons against parametric fits, the histogram bars should be shifted so that they are centered on the correct "nearest half hour" values rather than starting or ending at them.

It would be best to get actual durations instead of these highly binned data. If that's not possible, as this is a time-to-event (end-of-task event) analysis, you might consider parametric survival models with the times properly treated as interval-censored. If you do use an interval-censoring approach, make sure to use the correct boundary times for intervals if the current data display is for the "closest" half hour. For example, cases in the bin shown here from 1 to 1.5 hours would be coded to represent durations either between 0.75 and 1.25 hours or 1.25 to 1.75 hours, depending on how you have done the time-rounding.

In a situation like this I would tend to trust the integration provided by the cumulative distribution over whatever you would get from the highly discretized "pdf" (more like a probability mass function here). From that cumulative plot I suspect that the Weibull, with its relationship to extreme-value distributions, will work pretty well as @BruceET suggests in a comment.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for answer! Just a note on the binning: even though the software for inputting these data allows the beginning and end times to be recorded down to the second, most operators discretize them in 20 or 30 minute intervals because the underlying events are not very easily delimited. So these highly binned data are what actually gets recorded in the database. $\endgroup$ – LmnICE Oct 9 at 0:15
  • $\begingroup$ Also most operators tend to round up to the nearest half hour, so I would expect to not find any “0 hours” instances. $\endgroup$ – LmnICE Oct 9 at 0:17
  • 1
    $\begingroup$ @LmnICE if the binning is unavoidable and rounding is up rather than to closest half hour, then you should at least specify the time durations more toward the midpoint of the ranges. For example, times now specified as 1 hour should be specified as 45 minutes. It would be better to use ranges and treat the durations as interval-censored; for example, with rounding up to the next half hour, times now specified as 1 hour should be coded as 30-60 minutes. That mis-coding of time is probably more of a problem than any AIC/K-S discrepancy for finding a useful parametric form. $\endgroup$ – EdM Oct 9 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.