Como determinar qual distribuição se encaixa melhor em meus dados?

> 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 

fits <- list(no = fitdistr(mydata, "normal"), we = fitdistr(mydata, "weibull")) sapply(fits, function(i) i$loglik) no we -259.6540 -257.9268 

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

• Os $ p $ -valores de um Kolmovorov -Smirnov-Test (KS-Test) com parâmetros estimados estará muito errado. Então, infelizmente, você não pode apenas ajustar uma distribuição e então usar os parâmetros estimados em um teste de Kolmogorov-Smirnov para testar sua amostra.
• Sua amostra nunca seguirá um teste específico distribuição exata. Portanto, mesmo se seus $ p $ -valores do teste KS fossem válidos e $ > 0,05 $ , significaria apenas que você não pode “descartar que seus dados seguem esta distribuição específica. Outra formulação seria que sua amostra é compatível com uma determinada distribuição. Mas a resposta à pergunta ” Meus dados seguem a distribuição xy exatamente? ” é sempre não.
• 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 explorar. Usarei o excelente fitdistrplus pacote que oferece algumas funções interessantes para ajuste de distribuição. Usaremos a função descdist para ganhar algumas idéias sobre possíveis distribuições candidatas.

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) 

A curtose e a assimetria quadrada de sua amostra são representadas como um ponto azul chamado ” Observação “. Parece que as distribuições possíveis incluem a distribuição Weibull, Lognormal e possivelmente a distribuição Gama.

Vamos ajustar um 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) 

E para o ajuste Weibull:

plot(fit.weibull) 

Ambos parecem bons, mas a julgar pelo QQ-Plot, o Weibull talvez pareça um pouco melhor, especialmente nas caudas. Correspondentemente, o AIC do ajuste de Weibull é menor em comparação com o 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() 

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

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 bootstrapping para adicionar intervalos de confiança pontuais para o Weibull PDF ou CDF 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) 

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

---

Ajuste de distribuição automática com GAMLSS

O gamlss pacote para R oferece a capacidade de experimentar muitas distribuições diferentes e selecionar o melhor ” de acordo com o GAIC (o critério de informação Akaike generalizado). A função principal é fitDist. Uma opção importante nesta função é o tipo de distribuições que são tentadas. Por exemplo, definir type = "realline" tentará todas as distribuições implementadas definidas em toda a linha real, enquanto type = "realsplus" tentará apenas distribuições definidas na linha real positiva . Outra opção importante é o parâmetro $ k $ , que é a penalidade para o GAIC. No exemplo abaixo, defino o parâmetro $ k = 2 $ , o que significa que o ” best ” distribuição é selecionada de acordo com o clássico AIC. Você pode definir $ k $ como qualquer coisa que desejar, como $ \ log (n) $ para o BIC.

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

De acordo com o 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 por olhando para os resíduos em um gráfico de worm (basicamente um gráfico QQ de tendência):

Esperamos que os resíduos estejam próximos da linha horizontal do meio e 95% deles estejam entre a linha superior e curvas pontilhadas inferiores, que atuam como intervalos de confiança pontuais de 95%. Neste caso, o gráfico do worm parece bom para mim, indicando que a distribuição Weibull é um ajuste adequado. = “comments”>

+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 de uma distribuição de mistura ‘ s presença? Ou precisamos realizar uma análise de mistura adequada e verificar GoF para descarta essa opção?

@AleksandrBlekh É impossível ter energia suficiente para descartar uma mistura: quando a mistura é de duas distribuições quase idênticas, ela não pode ser detectada e quando todos os componentes, exceto um, têm proporções muito pequenas ele também não pode ser detectado. Normalmente (na ausência de uma teoria que possa sugerir uma forma de distribuição), ajusta-se as distribuições paramétricas para obter descrições aproximadas parcimoniosas dos dados. As misturas não são nada disso: elas requerem muitos parâmetros e são muito flexíveis para o propósito.

@whuber: +1 Agradeço sua excelente explicação!

@Lourenco Eu olhei para o gráfico de Cullen e Fey. O ponto azul denota nossa amostra. Você vê que o ponto está próximo às linhas de Weibull, Lognormal e Gamma (que está entre Weibull e Gamma). Depois de ajustar cada uma dessas distribuições, comparei as estatísticas de adequação usando a função gofstat e o AIC. Não existe um ‘ consenso sobre qual é a melhor maneira de determinar a ” melhor ” distribuição é. Gosto de métodos gráficos e do AIC.

@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 ‘ d normalmente consideraria. Para este tutorial, eu ‘ optei por não mostrá-lo para manter a postagem curta. O lognormal mostra um ajuste pior em comparação com a distribuição Weibull e Normal. O AIC é 537,59 e os gráficos também ‘ não parecem muito bons.

Os gráficos são, em geral, uma boa maneira de ter uma ideia melhor de como seus dados se parecem. No seu caso, eu recomendaria traçar a função de distribuição cumulativa empírica (ecdf) contra os cdfs teóricos com os parâmetros que você obteve de fitdistr ().

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

A linha preta é a função de distribuição cumulativa empírica e as linhas coloridas são cdfs de distribuições diferentes usando parâmetros que obtive usando o método de máxima verossimilhança. Pode-se ver facilmente que as distribuições exponencial e normal não são um bom ajuste para os dados, porque as linhas têm uma forma diferente da ecdf e as linhas estão bem longe 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, escolher aquelas que parecem se encaixar melhor no gráfico), você pode usar um Teste de razão de verossimilhança para testar quais distribuições se adaptam melhor.

Comentários

• Bem-vindo ao CrossValidated! Sua resposta poderia ser mais útil se você pudesse editá-la para incluir (a) o código usado para produzir o gráfico e (b) como se leria o gráfico.
• O que está sendo plotado lá? Isso é algum tipo de gráfico de exponencialidade?
• Mas como você decide qual distribuição se ajusta melhor aos seus dados? Apenas de acordo com o gráfico, eu não poderia ‘ dizer se logNormal ou weibull se ajusta melhor aos seus dados.
• Se você deseja criar um gerador de números pseudo-aleatórios, por que não usa o cdf empírico? Você deseja desenhar números que vão além de sua distribuição observada?
• Tomando seu gráfico pelo valor de face, pareceria que nenhuma de suas distribuições candidatas se ajustam bem aos dados. Além disso, seu ecdf parece ter uma assíntota horizontal inferior a 0,03, o que não ‘ não faz sentido, então ‘ não tenho certeza de que é realmente um ecdf em primeiro lugar.