2

Estou usando o fitdist para estimar os parâmetros da distribuição beta e ajustá-los a um conjunto de dados da seguinte forma como apresenta os códigos. Eu gostaria de entender por que não está se ajustando aos dados. Os códigos estão anexados. A figura mostra a não-variação da curva no histograma que há de tras, o que pode ser?

Data:

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)


enter image description here

1 Resposta 1

2

O erro está em tentar ajustar uma beta a dados que não estão no intervalo [0, 1], o suporte da distribuição beta. Deve-se primeiro transformar os dados.

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)

inserir a descrição da imagem aqui

O ajuste da distribuição beta é feito com o método dos momentos porque há valores para os quais dbeta dá zero, e portanto loglik é igual a -Inf. Um modo de resolver isso será eliminar os valores em causa.

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

Os valores dos parâmetros são semelhantes mas não são exatamente iguais. Agora é só usar o método acima para traçar o gráfico.

Dados.

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
")
3
  • Olá Rui, tudo bem? Muito obrigado! Sabe me dizer por qual motivo o AIC está dando infinito?
    – user55546
    19/06/2019 às 11:44
  • @BrenoGabriel Dei a resposta logo a seguir ao gráfico. Espero que esteja clara. 19/06/2019 às 12:42
  • @BrenoGabriel Só mais uma nota: a pergunta está bastante completa e por isso votei para cima. Mas o exemplo não é mínimo, bastava carregar um ou dois pacotes (MASS e fitdistrplus) e dar um ou dois exemplos. O exemplo onde está o erro e talvez outro exemplo onde corre bem para vermos o que se pretende 19/06/2019 às 12:51

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .