Fitdist problems in beta distribution

Question

I am using fitdist to estimate the parameters of the beta distribution and adjust them to a data set as follows as it presents the codes. I would like to understand why it is not adjusting to the data. The codes are attached. The figure shows the non-variation of the curve in the background histogram, which may be?

Date:

dados1 <- scan(text = "
36.6 33.8 37.8 34.0 32.7 31.8 31.7 36.5 32.7
  33.0 36.9 32.6 33.2 34.0 34.3 33.0 30.4 30.0
  30.6 30.2 33.6 35.8 35.6 36.8 36.8 33.3 33.2
  35.2 35.4 35.0 36.4 36.8 37.4 32.6 32.8 31.4
  31.5 34.4 35.2 38.5 38.3 35.9 37.7 34.0 35.6
  35.1 32.4 35.6 34.5 34.6 34.7 34.3 32.4 31.8
  31.8 36.4 34.1 35.7 34.4 37.1 35.0 31.0 36.5
  28.8 28.3 29.4 28.3 30.9 31.3 35.4 34.0 35.9
  33.4 33.6 33.3 30.0 32.6 27.0 27.0 26.4 23.8
  24.4 26.8 27.6 30.2 28.7 30.4 34.4 35.6 31.0
  33.2 36.6 37.9 34.5 35.0 31.5 37.9 36.5 31.0
 32.0 32.5 36.2 35.3 33.6 31.9 27.8 31.2 31.8
 35.1 36.6 36.8 31.7 30.0 31.5 32.2 34.9 35.7
 38.2 38.5 36.2 33.4 33.0 32.0 31.8
")

library(Hmisc) 
library(agricolae)
library(moments)
library(car)
library(MASS)
library(hnp)
library(fitdistrplus)
library(ggplot2)
library(grid)
library(fBasics)
library(VGAM)

dados1=dados$TempMaxima
dados1

####################Estimação dos Modelos########################
Gamm1 = fitdist(data = dados1, distr = "gamma")
summary(Gamm1)

Weibull1 = fitdist(data = dados1, distr = "weibull")
summary(Weibull1)

lnorm1 = fitdist(data = dados1, distr = "lnorm")
summary(lnorm1)

beta1 = fitdist((data=dados1)/40, distr="beta")
summary(beta1)

rm(dgumbel) ## get rid of previous definition
## hack behaviour of VGAM::pgumbel() a little bit
pgumbel <- function(x,...) {
  if (length(x)==0) numeric(0) else VGAM::pgumbel(x,...)
}

gumbel1 <- fitdist(dados1, "gumbel", 
                     start=list(location=10, scale=10))
summary(gumbel1)


norm1 = fitdist(data = dados1, distr = "norm")
summary(norm1) 

########################### Graphics ###########################
x11()
par(mfrow=c(2,3))
hist(dados1, probability = T, ylab = NULL, 
     main = "Distribuição Gamma", xlab = NULL, ylim = c(0,0.15),cex = 1.5)
curve(dgamma(x, shape=Gamm1$estimate[1], rate=Gamm1$estimate[2]), 
      add=T, lwd = 2, lty = 5, col ="red")

hist(dados1, probability = T, ylab = NULL, ylim = c(0,0.15),
     main = "Distribuição Weibull ", xlab = NULL, cex = 1.5)
curve(dweibull(x, shape=Weibull1$estimate[1], scale=Weibull1$estimate[2]), 
      add=T, lwd = 2, col ="red")

hist(dados1, probability = T, ylab = NULL, ylim = c(0,0.15),
     main = "Distribuição Log-Normal ", xlab = NULL, cex = 1.5)
curve(dlnorm(x, lnorm1$estimate[1], lnorm1$estimate[2]), 
      add=T, lwd = 2, lty = 3, col ="red")

hist(dados1, probability = T, ylab = NULL,  ylim = c(0,0.15),
     main = "Distribuição Gumbel I", xlab = NULL, cex = 1.5)
curve(dgumbel(x, gumbel1$estimate[1], gumbel1$estimate[2]), 
      add=T, lwd = 2, lty = 2, col ="red")

hist(dados1, probability = T, ylab = NULL,  ylim = c(0,0.15),
     main = "Distribuição Normal", xlab = NULL, cex = 1.5)
curve(dnorm(x, norm1$estimate[1], norm1$estimate[2]), 
      add=T, lwd = 2,lty = 4, col ="red")

hist(dados1, probability = T, ylab = NULL,
     main = "Distribuição Beta", xlab = NULL, cex = 1.5)
curve(dbeta(x, beta1$estimate[1], beta1$estimate[2]),add=T,lwd = 2,lty = 4)

Answer 1

The error is in trying to adjust a beta to data that is not in the range [0, 1], the beta distribution support. The data must first be transformed.

ScaleData <- function(x, na.rm = TRUE) {
  (x - min(x, na.rm = na.rm))/diff(range(x, na.rm = na.rm))
}

beta1 <- fitdistrplus::fitdist(ScaleData(dados1), 
                               distr = "beta",
                               method = "mme")
hist(ScaleData(dados1), prob = TRUE)
curve(dbeta(x, shape1 = beta1$estimate[1], shape2 = beta1$estimate[2]), 
      from = 0, to = 1,
      add = TRUE, lwd = 2, lty = 4)

The adjustment of the beta distribution is done with the moments method because there are values for which dbeta is zero, and therefore loglik is equal to -Inf. One way to resolve this is to eliminate the values in question.

i <- dbeta(ScaleData(dados1), 
           shape1 = beta1$estimate[1], 
           shape2 = beta1$estimate[2]) > 0

beta2 <- fitdistrplus::fitdist(ScaleData(dados1)[i], 
                               distr = "beta",
                               method = "mle")

beta2
#Fitting of the distribution ' beta ' by maximum likelihood 
#Parameters:
#       estimate Std. Error
#shape1 3.082147  0.3921171
#shape2 1.690168  0.2015881

The values of the parameters are similar but not exactly equal. Now just use the above method to plot the graph.

Dice.

dados1 <- scan(text = "
36.6 33.8 37.8 34.0 32.7 31.8 31.7 36.5 32.7
  33.0 36.9 32.6 33.2 34.0 34.3 33.0 30.4 30.0
  30.6 30.2 33.6 35.8 35.6 36.8 36.8 33.3 33.2
  35.2 35.4 35.0 36.4 36.8 37.4 32.6 32.8 31.4
  31.5 34.4 35.2 38.5 38.3 35.9 37.7 34.0 35.6
  35.1 32.4 35.6 34.5 34.6 34.7 34.3 32.4 31.8
  31.8 36.4 34.1 35.7 34.4 37.1 35.0 31.0 36.5
  28.8 28.3 29.4 28.3 30.9 31.3 35.4 34.0 35.9
  33.4 33.6 33.3 30.0 32.6 27.0 27.0 26.4 23.8
  24.4 26.8 27.6 30.2 28.7 30.4 34.4 35.6 31.0
  33.2 36.6 37.9 34.5 35.0 31.5 37.9 36.5 31.0
 32.0 32.5 36.2 35.3 33.6 31.9 27.8 31.2 31.8
 35.1 36.6 36.8 31.7 30.0 31.5 32.2 34.9 35.7
 38.2 38.5 36.2 33.4 33.0 32.0 31.8
")