Comment déterminer quelle distribution correspond le mieux à mes données?

> 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 

Tout dabord, voici quelques commentaires rapides:

• Les $ p $ -valeurs dun Kolmovorov -Smirnov-Test (KS-Test) avec des paramètres estimés sera tout à fait faux. Donc, malheureusement, vous ne pouvez pas simplement ajuster une distribution et ensuite utiliser les paramètres estimés dans un test Kolmogorov-Smirnov pour tester votre échantillon.
• Votre échantillon ne suivra jamais un distribution exacte. Ainsi, même si les valeurs de $ p $ du test KS seraient valides et $ > 0,05 $ , cela signifierait simplement que vous ne pouvez pas exclure que vos données suivent cette distribution spécifique. Une autre formulation serait que votre échantillon est compatible avec une certaine distribution. Mais la réponse à la question  » Mes données suivent-elles exactement la distribution xy?  » est toujours non.
• Le but ici ne peut pas être de déterminer avec certitude quelle distribution suit votre échantillon. Le but est ce que @whuber (dans les commentaires) appelle des descriptions approximatives parcimonieuses des données. Avoir une distribution paramétrique spécifique peut être utile comme modèle des données.

---

Mais allons faire un peu d’exploration. Je vais utiliser l’excellent fitdistrplus qui offre quelques fonctions intéressantes pour lajustement de la distribution. Nous utiliserons la fonction descdist pour gagner quelques idées sur déventuelles distributions candidates.

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) 

Utilisons maintenant descdist:

descdist(x, discrete = FALSE) 

Le kurtosis et lasymétrie au carré de votre échantillon sont tracés sous forme de point bleu nommé  » Observation « . Il semble que les distributions possibles incluent la distribution Weibull, Lognormal et éventuellement la distribution Gamma.

Soit a Distribution Weibull et une distribution normale:

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

Maintenant, inspectez lajustement pour la normale:

plot(fit.norm) 

Et pour le fit Weibull:

plot(fit.weibull) 

Les deux ont lair bien mais à en juger par le QQ-Plot, le Weibull est peut-être un peu meilleur, surtout au niveau des queues. En conséquence, lAIC de lajustement de Weibull est inférieur à celui de lajustement normal:

fit.weibull$aic [1] 519.8537 fit.norm$aic [1] 523.3079 

---

Simulation de test de Kolmogorov-Smirnov

Jutiliserai la procédure de @Aksakal « expliquée ici pour simuler la statistique KS sous la valeur nulle.

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

Le ECDF des statistiques KS simulées ressemble à ceci:

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

Enfin, notre valeur de $ p $ en utilisant la distribution nulle simulée des statistiques KS est:

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

Ceci confirme notre conclusion graphique que léchantillon est compatible avec une distribution Weibull.

Comme expliqué ici , nous pouvons utiliser le bootstrap pour ajouter des intervalles de confiance ponctuels au PDF ou CDF de Weibull estimé:

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

---

Ajustement de la distribution automatique avec GAMLSS

Le gamlss package pour R offre la possibilité dessayer de nombreuses distributions différentes et de sélectionner le meilleur  » selon le GAIC (le critère dinformation généralisé dAkaike). La fonction principale est fitDist. Une option importante dans cette fonction est le type des distributions qui sont essayées. Par exemple, définir type = "realline" essaiera toutes les distributions implémentées définies sur toute la ligne réelle tandis que type = "realsplus" essaiera uniquement les distributions définies sur la ligne positive réelle . Une autre option importante est le paramètre $ k $ , qui est la pénalité pour le GAIC. Dans lexemple ci-dessous, jai défini le paramètre $ k = 2 $ ce qui signifie que le  » best est sélectionnée selon lAIC classique. Vous pouvez définir $ k $ sur tout ce que vous voulez, par exemple $ \ log (n) $ pour le 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 *** 

Selon lAIC, la distribution de Weibull (plus précisément WEI2, une paramétrisation spéciale de celle-ci ) correspond le mieux aux données. Le paramétrage exact de la distribution WEI2 est détaillé dans ce document à la page 279. Inspectons lajustement en en regardant les résidus dans un graphique de vers (essentiellement un graphique QQ dépourvu de tendance):

Nous nous attendons à ce que les résidus soient proches de la ligne horizontale médiane et 95% dentre eux se situent entre la ligne supérieure et des courbes en pointillés inférieurs, qui agissent comme des intervalles de confiance à 95%. Dans ce cas, le tracé du ver me semble correct, indiquant que la distribution de Weibull est un ajustement adéquat.

Commentaires

• +1 Belle analyse. Une question, cependant. Une conclusion positive sur la compatibilité avec une distribution majeure particulière (Weibull, dans ce cas) permet dexclure la possibilité dune distribution mixte ‘ s? Ou nous devons effectuer une analyse de mélange appropriée et vérifier GoF pour exclure cette option?
• @AleksandrBlekh Il est impossible davoir assez de puissance pour exclure un mélange: lorsque le mélange est de deux distributions presque identiques, il ne peut pas être détecté et lorsque tous les composants sauf un ont de très petites proportions il ne peut pas non plus être détecté. Typiquement (en labsence dune théorie qui pourrait suggérer une forme distributionnelle), on ajuste des distributions paramétriques afin dobtenir des descriptions approchées parcimonieuses des données. Les mélanges nen font pas partie: ils nécessitent trop de paramètres et sont trop flexibles à cet effet.
• @whuber: +1 Appréciez votre excellente explication!
• @Lourenco Jai regardé le graphique de Cullen et Fey. Le point bleu désigne notre échantillon. Vous voyez que le point est proche des lignes de Weibull, Lognormal et Gamma (qui est entre Weibull et Gamma). Après avoir ajusté chacune de ces distributions, jai comparé les statistiques de qualité dajustement en utilisant la fonction gofstat et lAIC. Il ny a pas ‘ de consensus sur la meilleure façon de déterminer la  » meilleure  » distribution est. Jaime les méthodes graphiques et lAIC.
• @Lourenco Voulez-vous dire le lognormal? La distribution logistique (le signe  » + « ) est assez éloignée des données observées. Le lognormal serait également un candidat que je ‘ regarder normalement. Pour ce didacticiel, jai ‘ choisi de ne pas lafficher afin de garder le message court. La log-normale montre un ajustement moins bon que la distribution de Weibull et la distribution normale. LAIC est 537,59 et les graphiques ne sont pas ‘ trop beaux.

Les graphiques sont généralement un bon moyen davoir une meilleure idée de lapparence de vos données. Dans votre cas, je recommanderais de tracer la fonction de distribution cumulative empirique (ecdf) par rapport aux cdfs théoriques avec les paramètres que vous avez obtenus de fitdistr ().

Je lai fait une fois pour mes données et jai également inclus les intervalles de confiance. Voici limage que jai obtenue en utilisant ggplot2 ().

La ligne noire est la fonction de distribution cumulative empirique et les lignes colorées sont des cdfs de différentes distributions utilisant des paramètres que jai obtenus en utilisant la méthode du maximum de vraisemblance. On peut facilement voir que les distributions exponentielle et normale ne correspondent pas bien aux données, car les lignes ont une forme différente de lecdf et les lignes sont assez éloignées de lecdf. Malheureusement, les autres distributions sont assez proches. Mais je dirais que la ligne logNormal est la plus proche de la ligne noire. En utilisant une mesure de distance (par exemple MSE), on pourrait valider lhypothèse.

Si vous navez que deux distributions concurrentes (par exemple en choisissant celles qui semblent correspondre le mieux au graphique), vous pouvez utiliser un Test du rapport de vraisemblance pour tester quelles distributions correspondent le mieux.

Commentaires

Bienvenue dans CrossValidated! Votre réponse pourrait être plus utile si vous pouviez la modifier pour inclure (a) le code que vous avez utilisé pour produire le graphique, et (b) comment on lirait le graphique.

Quest-ce qui y est tracé? Sagit-il dune sorte de graphique dexponentialité?

Mais comment décidez-vous quelle distribution correspond le mieux à vos données? Ce nest que daprès le graphique que je nai pas pu ‘ vous dire si logNormal ou weibull correspond le mieux à vos données.

Si vous voulez créer un générateur de nombres pseudo-aléatoires pourquoi pas utiliser le CDF empirique? Voulez-vous dessiner des nombres qui vont au-delà de votre distribution observée?

En prenant votre graphique à sa valeur nominale, il semblerait que aucune de vos distributions candidates ne corresponde bien aux données. De plus, votre ecdf semble avoir une asymptote horizontale inférieure à 0,03, ce qui ‘ na aucun sens, donc je ‘ ne suis pas sûr que cest vraiment un ecdf en premier lieu.