I have a dataset and would like to figure out which distribution fits my data best.

I used the fitdistr() function to estimate the necessary parameters to describe the assumed distribution (i.e. Weibull, Cauchy, Normal). Using those parameters I can conduct a Kolmogorov-Smirnov Test to estimate whether my sample data is from the same distribution as my assumed distribution.

If the p-value is > 0.05 I can assume that the sample data is drawn from the same distribution. But the p-value doesn't provide any information about the godness of fit, isn't it?

So in case the p-value of my sample data is > 0.05 for a normal distribution as well as a weibull distribution, how can I know which distribution fits my data better?

This is basically the what I have done:

> mydata
 [1] 37.50 46.79 48.30 46.04 43.40 39.25 38.49 49.51 40.38 36.98 40.00
[12] 38.49 37.74 47.92 44.53 44.91 44.91 40.00 41.51 47.92 36.98 43.40
[23] 42.26 41.89 38.87 43.02 39.25 40.38 42.64 36.98 44.15 44.91 43.40
[34] 49.81 38.87 40.00 52.45 53.13 47.92 52.45 44.91 29.54 27.13 35.60
[45] 45.34 43.37 54.15 42.77 42.88 44.26 27.14 39.31 24.80 16.62 30.30
[56] 36.39 28.60 28.53 35.84 31.10 34.55 52.65 48.81 43.42 52.49 38.00
[67] 38.65 34.54 37.70 38.11 43.05 29.95 32.48 24.63 35.33 41.34

# estimate shape and scale to perform KS-test for weibull distribution
> fitdistr(mydata, "weibull")
     shape        scale   
   6.4632971   43.2474500 
 ( 0.5800149) ( 0.8073102)

# KS-test for weibull distribution
> ks.test(mydata, "pweibull", scale=43.2474500, shape=6.4632971)

        One-sample Kolmogorov-Smirnov test

data:  mydata
D = 0.0686, p-value = 0.8669
alternative hypothesis: two-sided

# KS-test for normal distribution
> ks.test(mydata, "pnorm", mean=mean(mydata), sd=sd(mydata))

        One-sample Kolmogorov-Smirnov test

data:  mydata
D = 0.0912, p-value = 0.5522
alternative hypothesis: two-sided

The p-values are 0.8669 for the Weibull distribution, and 0.5522 for the normal distribution. Thus I can assume that my data follows a Weibull as well as a normal distribution. But which distribution function describes my data better?


Referring to elevendollar I found the following code, but don't know how to interpret the results:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
  • 5
    Why would you like to figure out which distribution fits your data best? – Roland Jan 8 '15 at 10:40
  • 5
    Because I want to generate pseudo-random numbers following the given distribution. – tobibo Jan 8 '15 at 11:21
  • 5
    You can't use KS to check whether a distribution with parameters found from the dataset matches the dataset. See #2 on this page for example, plus alternatives (and other ways the KS test can be misleading). – tpg2114 Jan 8 '15 at 13:08
  • Another discussion here with code samples on how to apply KS test when parameters are estimated from the sample. – Aksakal Jan 8 '15 at 18:42
  • I used the fitdistr() function .....What's fitdistr function? Something from Excel? Or something you wrote yourself in C? – wolfies Nov 13 '15 at 14:02
up vote 130 down vote accepted
+50

First, here are some quick comments:

  • The $p$-values of a Kolmovorov-Smirnov-Test (KS-Test) with estimated parameters will be quite wrong. So unfortunately, you can't just fit a distribution and then use the estimated parameters in a Kolmogorov-Smirnov-Test to test your sample.
  • Your sample will never follow a specific distribution exactly. So even if your $p$-values from the KS-Test would be valid and $>0.05$, it would just mean that you can't rule out that your data follow this specific distribution. Another formulation would be that your sample is compatible with a certain distribution. But the answer to the question "Does my data follow the distribution xy exactly?" is always no.
  • The goal here cannot be to determine with certainty what distribution your sample follows. The goal is what @whuber (in the comments) calls parsimonious approximate descriptions of the data. Having a specific parametric distribution can be useful as a model of the data.

But let's do some exploration. I will use the excellent fitdistrplus package which offers some nice functions for distribution fitting. We will use the functiondescdist to gain some ideas about possible candidate distributions.

library(fitdistrplus)
library(logspline)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

Now lets use descdist:

descdist(x, discrete = FALSE)

Descdist

The kurtosis and squared skewness of your sample is plottet as a blue point named "Observation". It seems that possible distributions include the Weibull, Lognormal and possibly the Gamma distribution.

Let's fit a Weibull distribution and a normal distribution:

fit.weibull <- fitdist(x, "weibull")
fit.norm <- fitdist(x, "norm")

Now inspect the fit for the normal:

plot(fit.norm)

Normal fit

And for the Weibull fit:

plot(fit.weibull)

Weibull fit

Both look good but judged by the QQ-Plot, the Weibull maybe looks a bit better, especially at the tails. Correspondingly, the AIC of the Weibull fit is lower compared to the normal fit:

fit.weibull$aic
[1] 519.8537

fit.norm$aic
[1] 523.3079

Kolmogorov-Smirnov test simulation

I will use @Aksakal's procedure explained here to simulate the KS-statistic under the null.

n.sims <- 5e4

stats <- replicate(n.sims, {      
  r <- rweibull(n = length(x)
                , shape= fit.weibull$estimate["shape"]
                , scale = fit.weibull$estimate["scale"]
  )
  as.numeric(ks.test(r
                     , "pweibull"
                     , shape= fit.weibull$estimate["shape"]
                     , scale = fit.weibull$estimate["scale"])$statistic
  )      
})

The ECDF of the simulated KS-statistics looks like follows:

plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7)
grid()

Simulated KS-statistics

Finally, our $p$-value using the simulated null distribution of the KS-statistics is:

fit <- logspline(stats)

1 - plogspline(ks.test(x
                       , "pweibull"
                       , shape= fit.weibull$estimate["shape"]
                       , scale = fit.weibull$estimate["scale"])$statistic
               , fit
)

[1] 0.8428425

This confirms our graphical conclusion that the sample is compatible with a Weibull distribution.

As explained here, we can use bootstrapping to add pointwise confidence intervals to the estimated Weibull PDF or CDF:

xs <- seq(10, 65, len=500)

true.weibull <- rweibull(1e6, shape= fit.weibull$estimate["shape"]
                         , scale = fit.weibull$estimate["scale"])

boot.pdf <- sapply(1:1000, function(i) {
  xi <- sample(x, size=length(x), replace=TRUE)
  MLE.est <- suppressWarnings(fitdist(xi, distr="weibull"))  
  dweibull(xs, shape=MLE.est$estimate["shape"],  scale = MLE.est$estimate["scale"])
}
)

boot.cdf <- sapply(1:1000, function(i) {
  xi <- sample(x, size=length(x), replace=TRUE)
  MLE.est <- suppressWarnings(fitdist(xi, distr="weibull"))  
  pweibull(xs, shape= MLE.est$estimate["shape"],  scale = MLE.est$estimate["scale"])
}
)   

#-----------------------------------------------------------------------------
# Plot PDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf),
     xlab="x", ylab="Probability density")
for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.pdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.pdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)

CI_Density

#-----------------------------------------------------------------------------
# Plot CDF
#-----------------------------------------------------------------------------

par(bg="white", las=1, cex=1.2)
plot(xs, boot.cdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.cdf),
     xlab="x", ylab="F(x)")
for(i in 2:ncol(boot.cdf)) lines(xs, boot.cdf[, i], col=rgb(.6, .6, .6, .1))

# Add pointwise confidence bands

quants <- apply(boot.cdf, 1, quantile, c(0.025, 0.5, 0.975))
min.point <- apply(boot.cdf, 1, min, na.rm=TRUE)
max.point <- apply(boot.cdf, 1, max, na.rm=TRUE)
lines(xs, quants[1, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[3, ], col="red", lwd=1.5, lty=2)
lines(xs, quants[2, ], col="darkred", lwd=2)
#lines(xs, min.point, col="purple")
#lines(xs, max.point, col="purple")

CI_CDF

  • 1
    +1 Nice analysis. One question, though. Does positive conclusion on compatibility with a particular major distribution (Weibull, in this case) allows to rule out a possibility of a mixture distribution's presence? Or we need to perform a proper mixture analysis and check GoF to rule out that option? – Aleksandr Blekh Jan 8 '15 at 19:00
  • 16
    @AleksandrBlekh It is impossible to have enough power to rule out a mixture: when the mixture is of two almost identical distributions it cannot be detected and when all but one component have very small proportions it cannot be detected, either. Typically (in the absence of a theory which might suggest a distributional form), one fits parametric distributions in order to achieve parsimonious approximate descriptions of data. Mixtures are none of those: they require too many parameters and are too flexible for the purpose. – whuber Jan 8 '15 at 19:04
  • 4
    @whuber: +1 Appreciate your excellent explanation! – Aleksandr Blekh Jan 8 '15 at 19:25
  • 1
    @Lourenco I looked at the Cullen and Fey graph. The blue point denotes our sample. You see that the point is close to the lines of the Weibull, Lognormal and Gamma (which is between Weibull and Gamma). After fitting each of those distributions, I compared the goodness-of-fit statistics using the function gofstat and the AIC. There isn't a consensus about what the best way to determine the "best" distribution is. I like graphical methods and the AIC. – COOLSerdash Jun 23 '17 at 17:33
  • 1
    @Lourenco Do you mean the lognormal? The logistic distribution (the "+" sign) is quite a bit away from the observed data. The lognormal would also be a candidate I'd normally look at. For this tutorial, I've chosen to not show it in order to keep the post short. The lognormal shows a worse fit compared to both the Weibull and Normal distribution. The AIC is 537.59 and the graphs also don't look too good. – COOLSerdash Jun 26 '17 at 19:19

Plots are mostly a good way to get a better idea of what your data looks like. In your case I would recommend plotting the empirical cumulative distribution function (ecdf) against the theoretical cdfs with the parameters you got from fitdistr().

I did that once for my data and also included the confidence intervals. Here is the picture I got using ggplot2().

enter image description here

The black line is the empirical cumulative distribution function and the colored lines are cdfs from different distributions using parameters I got using the Maximum Likelihood method. One can easily see that the exponential and normal distribution are not a good fit to the data, because the lines have a different form than the ecdf and lines are quite far away from the ecdf. Unfortunately the other distribtions are quite close. But I would say that the logNormal line is the closest to the black line. Using a measure of distance (for example MSE) one could validate the assumption.

If you only have two competing distributions (for example picking the ones that seem to fit best in the plot) you could use a Likelihood-Ratio-Test to test which distributions fits better.

  • 17
    Welcome to CrossValidated! Your answer might be more useful if you could edit it to include (a) the code you used to produce the graphic, and (b) how one would read the graphic. – Stephan Kolassa Jan 8 '15 at 10:49
  • 2
    What is being plotted there? Is that some sort of exponentiality-plot? – Glen_b Jan 8 '15 at 11:00
  • 1
    But how do you decide which distribution fits your data best? Only according to the graphic I couldn't tell you whether logNormal or weibull fits your data best. – tobibo Jan 8 '15 at 11:27
  • 4
    If you want to create a pseudo-random numbers generator why not use the empirical cdf? Do you want to draw numbers that go beyond your observed distribution? – elevendollar Jan 8 '15 at 12:14
  • 6
    Taking your graph at face value, it would appear that none of your candidate distributions fit the data well at all. Also, your ecdf appears to have a horizontal asymptote at less than 0.03 which doesn't make sense, so I'm not sure that it really is an ecdf in the first place. – Hong Ooi Jan 8 '15 at 15:56

protected by Nick Cox May 12 '17 at 13:28

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

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