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

133

J'ai un jeu de données et j'aimerais savoir quelle distribution correspond le mieux à mes données.

J'ai utilisé le fitdistr() fonction pour estimer les paramètres nécessaires pour décrire la distribution supposée (c.-à-d. Weibull, Cauchy, Normal). En utilisant ces paramètres, je peux effectuer un test de Kolmogorov-Smirnov pour estimer si les données de mon échantillon proviennent de la même distribution que ma distribution supposée.

Si la valeur p est> 0,05, je peux supposer que les données de l'échantillon proviennent de la même distribution. Mais la valeur p ne fournit aucune information sur la divinité de l'ajustement, n'est-ce pas?

Donc, si la valeur p de mes données d'échantillon est> 0,05 pour une distribution normale ainsi que pour une distribution de Weibull, comment puis-je savoir quelle distribution correspond le mieux à mes données?

C'est fondamentalement ce que j'ai fait:

> 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

Les valeurs de p sont 0,8669 pour la distribution de Weibull et 0,5522 pour la distribution normale. Je peux donc supposer que mes données suivent une distribution de Weibull aussi bien qu'une distribution normale. Mais quelle fonction de distribution décrit le mieux mes données?


En se référant à elevendollar, j'ai trouvé le code suivant, mais je ne sais pas comment interpréter les résultats:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
Tobibo
la source
5
Pourquoi voudriez-vous savoir quelle distribution correspond le mieux à vos données?
Roland
6
Parce que je veux générer des nombres pseudo-aléatoires à la suite de la distribution donnée.
lundi
6
Vous ne pouvez pas utiliser KS pour vérifier si une distribution avec des paramètres trouvés dans l'ensemble de données correspond à l'ensemble de données. Voir le n ° 2 sur cette page par exemple, plus les alternatives (et les autres manières par lesquelles le test de KS peut être trompeur).
tpg2114
Une autre discussion ici avec des exemples de code sur la façon d'appliquer le test KS lorsque des paramètres sont estimés à partir de l'échantillon.
Aksakal
1
I used the fitdistr() function ..... Quelle est la fitdistrfonction? Quelque chose d'Excel? Ou quelque chose que vous avez écrit en C?
loups

Réponses:

163

Tout d'abord, voici quelques commentaires rapides:

  • p
  • p>0.05
  • L’objectif ici ne peut pas être de déterminer avec certitude la répartition de votre échantillon. L'objectif 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 en tant que modèle des données.

Mais faisons de l'exploration. Je vais utiliser l'excellent fitdistrpluspaquet qui offre quelques fonctions intéressantes pour l'ajustement de la distribution. Nous utiliserons cette fonction descdistpour obtenir des idées sur les distributions de candidats possibles.

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)

Maintenant, utilisons descdist:

descdist(x, discrete = FALSE)

Descdiste

Le kurtosis et l'asymétrie au carré de votre échantillon sont représentés par un point bleu appelé "Observation". Il semble que les distributions possibles incluent la distribution de Weibull, Lognormal et éventuellement la distribution Gamma.

Adaptons une distribution de Weibull et une distribution normale:

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

Maintenant, vérifiez l'ajustement pour la normale:

plot(fit.norm)

Ajustement normal

Et pour la coupe Weibull:

plot(fit.weibull)

Ajustement de Weibull

Les deux ont l'air bien, mais à en juger par le QQ-Plot, le Weibull est peut-être un peu mieux, surtout au niveau des queues. De même, l’AIC de l’ajustement de Weibull est inférieur à l’ajustement normal:

fit.weibull$aic
[1] 519.8537

fit.norm$aic
[1] 523.3079

Simulation de test de Kolmogorov-Smirnov

Je vais utiliser la procédure de @ Aksakal expliquée ici pour simuler la statistique KS sous le zéro.

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 fichier ECDF des statistiques KS simulées se présente comme suit:

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

KS-statistiques simulées

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

Cela confirme notre conclusion graphique selon laquelle l'échantillon est compatible avec une distribution de Weibull.

Comme expliqué ici , nous pouvons utiliser l’amorçage pour ajouter des intervalles de confiance ponctuels au fichier PDF ou CDF estimé de Weibull:

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


Raccord de distribution automatique avec 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 ***

Selon l'AIC, la distribution de Weibull (plus précisément WEI2, une paramétrisation spéciale de celle-ci) correspond le mieux aux données. La paramétrisation exacte de la distribution WEI2est détaillée dans ce document à la page 279. Inspectons l'ajustement en examinant les résidus dans un graphe de vers (en gros, un graphe QQ hors tension):

WormPlot

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

COOLSerdash
la source
1
+1 belle analyse. Une question cependant. Une conclusion positive sur la compatibilité avec une distribution majeure particulière (Weibull, dans ce cas) permet-elle d'exclure une possibilité de présence d'une distribution de mélange? Ou nous devons effectuer une analyse appropriée du mélange et vérifier que le GoF exclut cette option?
Aleksandr Blekh
18
@AleksandrBlekh Il est impossible d'avoir assez de puissance pour exclure un mélange: lorsque le mélange est composé 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, ils ne le sont pas non plus. Typiquement (en l'absence d'une théorie qui pourrait suggérer une forme de distribution), on adapte les distributions paramétriques afin d'obtenir des descriptions approximatives parcimonieuses des données. Les mélanges n'en font pas partie: ils nécessitent trop de paramètres et sont trop souples pour le but recherché.
whuber
4
@whuber: +1 Appréciez votre excellente explication!
Aleksandr Blekh
1
@Lourenco J'ai 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 (entre Weibull et Gamma). Après avoir ajusté chacune de ces distributions, j'ai comparé les statistiques sur la qualité de l'ajustement à l'aide de la fonction gofstatet de l'AIC. Il n'y a pas de consensus sur le meilleur moyen de déterminer la "meilleure" distribution. J'aime les méthodes graphiques et l'AIC.
COOLSerdash
1
@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 regarderais normalement. Pour ce tutoriel, j'ai choisi de ne pas l'afficher afin de garder le message court. La log-normale montre un meilleur ajustement par rapport à la distribution de Weibull et à la distribution normale. L’AIC porte le numéro 537.59 et les graphiques ne sont pas très beaux.
COOLSerdash
15

Les graphiques sont généralement un bon moyen d’avoir une meilleure idée de ce à quoi ressemblent 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 l'ai fait une fois pour mes données et j'ai également inclus les intervalles de confiance. Voici l'image que j'ai obtenue en utilisant ggplot2 ().

entrez la description de l'image ici

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 j'ai obtenus en utilisant la méthode du maximum de vraisemblance. On peut facilement voir que les distributions exponentielle et normale ne correspondent pas aux données, car les lignes ont une forme différente de celle du fichier ecdf et les lignes sont assez éloignées du fichier ecdf. 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 l'hypothèse.

Si vous ne disposez que de deux distributions concurrentes (par exemple, en choisissant celles qui semblent correspondre le mieux à l'intrigue), vous pouvez utiliser un test de vraisemblance-ratio pour tester les distributions qui conviennent le mieux.

elevendollar
la source
20
Bienvenue sur CrossValidated! Votre réponse pourrait être plus utile si vous pouviez l'éditer pour inclure (a) le code que vous avez utilisé pour produire le graphique et (b) comment lire le graphique.
Stephan Kolassa
2
Qu'est-ce qui est tracé ici? Est-ce une sorte de tracé d'exponentialité?
Glen_b
1
Mais comment décidez-vous quelle distribution correspond le mieux à vos données? Selon le graphique, je ne saurais vous dire si logNormal ou weibull correspond le mieux à vos données.
lundi
4
Si vous voulez créer un générateur de nombres pseudo-aléatoires, pourquoi ne pas utiliser le fichier cdf empirique? Voulez-vous dessiner des nombres qui vont au-delà de votre distribution observée?
Elevendollar
6
En prenant votre graphique à sa valeur nominale, il semblerait qu'aucune de vos distributions candidates ne corresponde bien aux données. En outre, votre fichier ecdf semble présenter une asymptote horizontale inférieure à 0,03, ce qui n’a aucun sens. Je ne suis donc pas sûr qu’il s’agisse vraiment d’un fichier ecdf.
Hong Ooi