# How to determine which distribution fits my data best?

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"))
[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()


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))

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)


#-----------------------------------------------------------------------------
# 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))

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")


• +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
• @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
• @whuber: +1 Appreciate your excellent explanation! – Aleksandr Blekh Jan 8 '15 at 19:25
• @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
• @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().

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.

• 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
• What is being plotted there? Is that some sort of exponentiality-plot? – Glen_b Jan 8 '15 at 11:00
• 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
• 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
• 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 CoxMay 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).