Como determinar qual distribuição se ajusta melhor aos meus dados?

133

Eu tenho um conjunto de dados e gostaria de descobrir qual distribuição se ajusta melhor aos meus dados.

Eu usei a fitdistr()função para estimar os parâmetros necessários para descrever a distribuição assumida (ou seja, Weibull, Cauchy, Normal). Usando esses parâmetros, posso realizar um Teste de Kolmogorov-Smirnov para estimar se meus dados de amostra são da mesma distribuição que minha distribuição assumida.

Se o valor-p for> 0,05, posso assumir que os dados da amostra são extraídos da mesma distribuição. Mas o valor-p não fornece nenhuma informação sobre a divindade do ajuste, não é?

Portanto, se o valor p dos meus dados de amostra for> 0,05 para uma distribuição normal e uma distribuição weibull, como posso saber qual distribuição se ajusta melhor aos meus dados?

Isto é basicamente o que eu fiz:

> 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

Os valores de p são 0,8669 para a distribuição Weibull e 0,5522 para a distribuição normal. Assim, posso assumir que meus dados seguem uma distribuição Weibull e uma distribuição normal. Mas qual função de distribuição descreve melhor meus dados?


Referindo- me a elevendollar , encontrei o seguinte código, mas não sei como interpretar os resultados:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
tobibo
fonte
5
Por que você gostaria de descobrir qual distribuição se ajusta melhor aos seus dados?
Roland
6
Porque eu quero gerar números pseudo-aleatórios seguindo a distribuição fornecida.
tobibo
6
Você não pode usar o KS para verificar se uma distribuição com parâmetros encontrados no conjunto de dados corresponde ao conjunto de dados. Veja o item 2 desta página, por exemplo, além de alternativas (e outras maneiras pelas quais o teste KS pode ser enganoso).
precisa saber é o seguinte
Outra discussão aqui com exemplos de código sobre como aplicar o teste KS quando parâmetros são estimados a partir da amostra.
Aksakal
1
I used the fitdistr() function ..... O que é fitdistrfunção? Algo do Excel? Ou algo que você escreveu em C?
wolfies

Respostas:

162

Primeiro, aqui estão alguns comentários rápidos:

  • Os valores de um teste de Kolmovorov-Smirnov (teste KS) com parâmetros estimados estarão completamente errados. Infelizmente, você não pode ajustar apenas uma distribuição e, em seguida, usar os parâmetros estimados em um teste de Kolmogorov-Smirnov para testar sua amostra.p
  • Sua amostra nunca seguirá exatamente uma distribuição específica. Portanto, mesmo que seus valores- do KS-Test sejam válidos e , isso significa apenas que você não pode descartar que seus dados seguem essa distribuição específica. Outra formulação seria que sua amostra é compatível com uma certa distribuição. Mas a resposta para a pergunta "Meus dados seguem exatamente a distribuição xy?" é sempre não.p>0.05
  • O objetivo aqui não pode ser determinar com certeza qual distribuição sua amostra segue. O objetivo é o que @whuber (nos comentários) chama de descrições aproximadas parcimoniosas dos dados. Ter uma distribuição paramétrica específica pode ser útil como modelo dos dados.

Mas vamos fazer alguma exploração. Usarei o excelente fitdistrpluspacote que oferece algumas funções interessantes para o ajuste da distribuição. Usaremos a função descdistpara obter algumas idéias sobre possíveis distribuições de candidatos.

library(fitdistrplus)
library(logspline)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

Agora vamos usar descdist:

descdist(x, discrete = FALSE)

Descdist

A curtose e assimetria quadrada da sua amostra são marcadas como um ponto azul chamado "Observação". Parece que as possíveis distribuições incluem a distribuição Weibull, Lognormal e possivelmente a gama.

Vamos ajustar uma distribuição Weibull e uma distribuição normal:

fit.weibull <- fitdist(x, "weibull")
fit.norm <- fitdist(x, "norm")

Agora inspecione o ajuste para o normal:

plot(fit.norm)

Ajuste normal

E para o ajuste Weibull:

plot(fit.weibull)

Weibull fit

Ambos parecem bons, mas julgados pelo QQ-Plot, o Weibull talvez pareça um pouco melhor, especialmente nas caudas. Do mesmo modo, o AIC do ajuste Weibull é menor em comparação ao ajuste normal:

fit.weibull$aic
[1] 519.8537

fit.norm$aic
[1] 523.3079

Simulação de teste de Kolmogorov-Smirnov

Usarei o procedimento de @ Aksakal explicado aqui para simular a estatística KS sob o valor nulo.

n.sims <- 5e4

stats <- replicate(n.sims, {      
  r <- rweibull(n = length(x)
                , shape= fit.weibull$estimate["shape"]
                , scale = fit.weibull$estimate["scale"]
  )
  estfit.weibull <- fitdist(r, "weibull") # added to account for the estimated parameters
  as.numeric(ks.test(r
                     , "pweibull"
                     , shape= estfit.weibull$estimate["shape"]
                     , scale = estfit.weibull$estimate["scale"])$statistic
  )      
})

O ECDF das estatísticas KS simuladas tem a seguinte aparência:

plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7)
grid()

Estatísticas KS simuladas

Finalmente, nosso valor- usando a distribuição nula simulada das estatísticas KS é:p

fit <- logspline(stats)

1 - plogspline(ks.test(x
                       , "pweibull"
                       , shape= fit.weibull$estimate["shape"]
                       , scale = fit.weibull$estimate["scale"])$statistic
               , fit
)

[1] 0.4889511

Isso confirma nossa conclusão gráfica de que a amostra é compatível com uma distribuição Weibull.

Conforme explicado aqui , podemos usar o bootstrapping para adicionar intervalos de confiança pontuais ao PDF ou CDF Weibull estimado:

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

# Add pointwise confidence bands

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)

CI_Density

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

# Add pointwise confidence bands

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

CI_CDF


Acessório de distribuição automática com GAMLSS

gamlssRfitDisttype = "realline"type = "realsplus"kk=2klog(n)

library(gamlss)
library(gamlss.dist)
library(gamlss.add)

x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00,
       38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40,
       42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40,
       49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60,
       45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30,
       36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00,
       38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34)

fit <- fitDist(x, k = 2, type = "realplus", trace = FALSE, try.gamlss = TRUE)

summary(fit)

*******************************************************************
Family:  c("WEI2", "Weibull type 2") 

Call:  gamlssML(formula = y, family = DIST[i], data = sys.parent()) 

Fitting method: "nlminb" 


Coefficient(s):
             Estimate  Std. Error  t value   Pr(>|t|)    
eta.mu    -24.3468041   2.2141197 -10.9962 < 2.22e-16 ***
eta.sigma   1.8661380   0.0892799  20.9021 < 2.22e-16 ***

Segundo a AIC, a distribuição Weibull (mais especificamente WEI2, uma parametrização especial dela) se ajusta melhor aos dados. A parametrização exata da distribuição WEI2é detalhada neste documento na página 279. Vamos inspecionar o ajuste observando os resíduos em um gráfico de sem - fim (basicamente um gráfico de QQ descendente):

WormPlot

Esperamos que os resíduos fiquem próximos da linha horizontal média e 95% deles fiquem entre as curvas pontilhadas superior e inferior, que atuam como intervalos de confiança de 95% no sentido dos pontos. Nesse caso, o gráfico de worms parece bom para mim, indicando que a distribuição Weibull é um ajuste adequado.

COOLSerdash
fonte
1
+1 boa análise. Uma pergunta, no entanto. A conclusão positiva sobre a compatibilidade com uma distribuição principal específica (Weibull, neste caso) permite descartar a possibilidade da presença de uma distribuição de mistura? Ou precisamos realizar uma análise adequada da mistura e verificar o GoF para descartar essa opção?
Aleksandr Blekh
18
@AleksandrBlekh É impossível ter poder suficiente para descartar uma mistura: quando a mistura é de duas distribuições quase idênticas, ela não pode ser detectada e quando todos, exceto um componente, têm proporções muito pequenas, também não podem ser detectados. Tipicamente (na ausência de uma teoria que possa sugerir uma forma distributiva), ajustamos as distribuições paramétricas para obter descrições aproximadas parcimoniosas dos dados. As misturas não são nenhuma delas: elas exigem muitos parâmetros e são muito flexíveis para a finalidade.
whuber
4
@whuber: +1 Aprecie sua excelente explicação!
Aleksandr Blekh
1
@Lourenco Olhei para o gráfico de Cullen e Fey. O ponto azul denota nossa amostra. Você vê que o ponto está próximo das linhas de Weibull, Lognormal e Gamma (que fica entre Weibull e Gamma). Depois de ajustar cada uma dessas distribuições, comparei as estatísticas de ajuste usando a função gofstate o AIC. Não há consenso sobre qual é a melhor maneira de determinar a "melhor" distribuição. Eu gosto de métodos gráficos e da AIC.
COOLSerdash
1
@Lourenco Você quer dizer o lognormal? A distribuição logística (o sinal "+") está um pouco distante dos dados observados. O lognormal também seria um candidato que eu normalmente consideraria. Para este tutorial, decidi não mostrá-lo para manter o post curto. O lognormal mostra um ajuste pior comparado à distribuição Weibull e Normal. A AIC é 537,59 e os gráficos também não parecem muito bons.
COOLSerdash
15

Os gráficos são principalmente uma boa maneira de ter uma idéia melhor da aparência dos seus dados. No seu caso, eu recomendaria plotar a função de distribuição cumulativa empírica (ecdf) contra os cdfs teóricos com os parâmetros obtidos em fitdistr ().

Fiz isso uma vez para os meus dados e também incluí os intervalos de confiança. Aqui está a foto que obtive usando ggplot2 ().

insira a descrição da imagem aqui

A linha preta é a função de distribuição cumulativa empírica e as linhas coloridas são cdfs de diferentes distribuições usando parâmetros que obtive usando o método de máxima verossimilhança. Pode-se ver facilmente que a distribuição exponencial e normal não se ajusta aos dados, porque as linhas têm uma forma diferente da ecdf e as linhas estão bem distantes da ecdf. Infelizmente, as outras distribuições são bastante próximas. Mas eu diria que a linha logNormal é a mais próxima da linha preta. Usando uma medida de distância (por exemplo, MSE), pode-se validar a suposição.

Se você tiver apenas duas distribuições concorrentes (por exemplo, escolhendo as que parecem se encaixar melhor na plotagem), poderá usar um Teste de Razão de Verossimilhança para testar quais distribuições se encaixam melhor.

elevendollar
fonte
20
Bem-vindo ao CrossValidated! Sua resposta pode ser mais útil se você puder editá-la para incluir (a) o código usado para produzir o gráfico e (b) como alguém leria o gráfico.
precisa saber é o seguinte
2
O que está sendo plotado lá? Isso é algum tipo de trama de exponencialidade?
Glen_b
1
Mas como você decide qual distribuição se ajusta melhor aos seus dados? Somente de acordo com o gráfico, não sei se o logNormal ou o weibull se encaixa melhor nos seus dados.
tobibo
4
Se você deseja criar um gerador de números pseudo-aleatórios, por que não usar o cdf empírico? Deseja desenhar números que vão além da distribuição observada?
Elevendollar
6
Tomando o seu gráfico pelo valor nominal, parece que nenhuma das suas distribuições candidatas se encaixa bem nos dados. Além disso, seu ecdf parece ter uma assíntota horizontal menor que 0,03, o que não faz sentido, então não tenho certeza de que seja realmente um ecdf em primeiro lugar.
Hong Ooi