0
$\begingroup$

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: enter image description here

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:

enter image description here

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?

$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$ – MrFlick Feb 1 at 20:21

Your Answer

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

Browse other questions tagged or ask your own question.