# How to choose and plot the most appropriate distribution in R?

I need to choose the distribution that best fits my data for different datasets. There are similar discussions here and here, but I am still struggling to find proper solution.

My first attempt was simple density plot:

#  Creating data set
TLT <- c(rep(0,32), rep(1,120), rep(2,10), rep(3,67), rep(4,14),  rep(5,7), 6)
DLT_Code <- c(rep('DLT_Code',251))

data_to_plot <- data.frame(cbind(DLT_Code,TLT))
data_to_plot$TLT <- as.numeric(as.character(data_to_plot$TLT ))

#  Plot data
ggplot(data_to_plot, aes(x = TLT)) +
geom_histogram(aes(y = ..density..), binwidth = 1, colour= "black", fill = "white") +
geom_density( aes(y=..density..), fill="blue", alpha = .25)


But it does not look good - too redundant to the histogram:

My next attempt was to create a function, calculating the best fit:

library(plyr)
library(dplyr)
library(fitdistrplus)
library(evd)
library(gamlss)
library(purrr)

fdistr <- function(d) {

TLT <- d$TLT if (sum(TLT<=0)) {TLT[TLT <= 0] <- 0.001} # removing value < 0 for log calculation gev <- fgev(TLT, std.err=FALSE) distr <- c('norm', 'lnorm', 'weibull', 'gamma') fit <- lapply(X=distr, FUN=fitdist, data=TLT) fit[[5]] <- gev distr[5] <- 'gev' names(fit) <- distr Loglike <- sapply(X=fit, FUN=logLik) Loglike_Best <- which(Loglike == max(Loglike)) x_data <- max(d$TLT)

mynorm <- function(x) {
dnorm(x,
mean=fit[['norm']]$estimate[1], sd=fit[['norm']]$estimate[2])
}

mylnorm <- function(x){
dlnorm(x,
meanlog=fit[['lnorm']]$estimate[1], sdlog=fit[['lnorm']]$estimate[2])
}

myweibull <- function(x) {
dweibull(x,
shape=fit[['weibull']]$estimate[1], scale=fit[['weibull']]$estimate[2])
}

mygamma <- function(x) {
dgamma(x,
shape=fit[['gamma']]$estimate[1], rate=fit[['gamma']]$estimate[2])
}

mygev <- function(x){
dgev(x,
loc=fit[['gev']]$estimate[1], scale=fit[['gev']]$estimate[2],
shape=fit[['gev']]$estimate[3]) } distributions <- c(mynorm, mylnorm, myweibull, mygamma, mygev) # get the max of each density y <- purrr::map_dbl(distributions, ~ optimize(., interval = c(0, x_data), maximum = T)$objective)

# find the max (excluding infinity)
ymax <- max(y[abs(y) < Inf])

hist(TLT, prob=TRUE, breaks= x_data,
main=paste(d$DLT_Code[1], '- best :', names(Loglike[Loglike_Best])), sub = 'Total Lead Times', col='lightgrey', border='white', ylim= c(0, ymax) ) lines(density(TLT), col='darkgrey', lty=2, lwd=2) grid(nx = NA, ny = NULL, col = "gray", lty = "dotted", lwd = .5, equilogs = TRUE) curve(mynorm, add=TRUE, col='blue', lwd=2, n = 1E5) curve(mylnorm, add=TRUE, col='darkgreen', lwd=2, n = 1E5) curve(myweibull, add=TRUE, col='purple', lwd=2, n = 1E5) curve(mygamma, add=TRUE, col='Gold', lwd=2, n = 1E5) curve(mygev, add=TRUE, col='red', lwd=2, n = 1E5) legend_loglik <- paste(c('Norm', 'LogNorm', 'Weibull', 'Gamma','GEV'), c(':'), round(Loglike, digits=2)) legend("topright", legend=legend_loglik, col=c('blue', 'darkgreen', 'purple', 'gold', 'red'), lty=1, lwd=2, bty='o', bg='white', box.lty=2, box.lwd = 1, box.col='white') return(data.frame(DLT_Code = d$DLT_Code[1],
n = length(d$TLT), Best = names(Loglike[Loglike_Best]), lnorm = Loglike[1], norm = Loglike[2], weibul = Loglike[3], gamma = Loglike[4], GEV = Loglike[5])) } # Creating data set TLT <- c(rep(0,32), rep(1,120), rep(2,10), rep(3,67), rep(4,14), rep(5,7), 6) DLT_Code <- c(rep('DLT_Code',251)) data_to_plot <- data.frame(cbind(DLT_Code,TLT)) data_to_plot$TLT <- as.numeric(as.character(data_to_plot$TLT )) DLT_Distr <- do.call(rbind, by(data = data_to_plot, INDICES = data_to_plot$DLT_Code, FUN=fdistr))


But then my best distribution fit will be looking like that:

But now distribution curve does not look like distribution (I was expecting something similar to bell shape). Why the curve is so high? Are there better ways to find and plot fitting distribution?

• There's a reason that those other conversations happened on Cross Validated. This is not really a specific programming question that's appropriate for Stack Overflow. Choosing how to model your data doesn't really depend on the code. You need to make modeling assumptions with your knowledge about the data and how it was generated and what sort of statistical tests you want to perform. Code can't tell you what the "correct" way to model your data is. You simply have to accept the benefits and drawbacks of whatever method you choose. – MrFlick Feb 1 at 20:21