Hoe bepaal ik welke distributie het beste bij mijn gegevens past?

> 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 

Allereerst zijn hier enkele korte opmerkingen:

• De $ p $ -waarden van een Kolmovorov -Smirnov-Test (KS-Test) met geschatte parameters zal behoorlijk fout zijn. Dus helaas kunt u “niet zomaar een distributie passen en vervolgens de geschatte parameters in een Kolmogorov-Smirnov-Test gebruiken om uw steekproef te testen.
• Uw steekproef zal nooit een specifieke exacte distributie. Dus zelfs als uw $ p $ -waarden uit de KS-Test geldig zijn en $ > 0,05 $ , zou het gewoon betekenen dat u niet kunt uitsluiten dat uw gegevens deze specifieke distributie volgen. Een andere formulering zou zijn dat uw monster compatibel is met een bepaalde distributie. Maar het antwoord op de vraag ” Volgen mijn gegevens precies de distributie xy? ” is altijd nee.
• Het doel hier kan niet zijn om met zekerheid te bepalen welke verdeling uw steekproef volgt. Het doel is wat @whuber (in de commentaren) spaarzame geschatte beschrijvingen van de gegevens noemt. Het hebben van een specifieke parametrische verdeling kan handig zijn als een model van de gegevens.

---

Maar laten we wat onderzoeken. Ik zal de uitstekende fitdistrplus pakket dat een aantal leuke functies biedt voor distributie-aanpassing. We zullen de functie descdist gebruiken om enkele ideeën over mogelijke kandidaat-distributies.

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) 

Laten we nu descdist gebruiken:

descdist(x, discrete = FALSE) 

De kurtosis en de gekwadrateerde scheefheid van uw steekproef wordt uitgezet als een blauw punt genaamd ” Observatie “. Het lijkt erop dat mogelijke distributies de Weibull, Lognormal en mogelijk de Gamma-distributie omvatten.

Laten we een Weibull-distributie en een normale distributie:

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

Inspecteer nu de fit voor de normale:

plot(fit.norm) 

En voor de Weibull fit:

plot(fit.weibull) 

Beide zien er goed uit, maar beoordeeld door de QQ-Plot, ziet de Weibull er misschien wat beter uit, vooral aan de staart. Dienovereenkomstig is de AIC van de Weibull-fit lager in vergelijking met de normale fit:

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

---

Kolmogorov-Smirnov testsimulatie

Ik zal de procedure van @Aksakal gebruiken die hier wordt uitgelegd om de KS-statistiek onder de nul te simuleren.

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

De ECDF van de gesimuleerde KS-statistieken ziet er als volgt uit:

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

Eindelijk, onze $ p $ -waarde met behulp van de gesimuleerde nulverdeling van de KS-statistieken is:

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

Dit bevestigt onze grafische conclusie dat de sample compatibel is met een Weibull-distributie.

Zoals uitgelegd hier kunnen we bootstrapping gebruiken om puntsgewijze betrouwbaarheidsintervallen toe te voegen aan de geschatte Weibull PDF of CDF:

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

---

Automatische distributie-aanpassing met GAMLSS

De gamlss pakket voor R biedt de mogelijkheid om veel verschillende distributies te proberen en het best ” volgens de GAIC (het gegeneraliseerde Akaike-informatiecriterium). De belangrijkste functie is fitDist. Een belangrijke optie in deze functie is het type distributies dat wordt geprobeerd. Als u bijvoorbeeld type = "realline" instelt, worden alle geïmplementeerde distributies geprobeerd die zijn gedefinieerd op de hele echte lijn, terwijl type = "realsplus" alleen distributies probeert die zijn gedefinieerd op de echte positieve lijn . Een andere belangrijke optie is de parameter $ k $ , wat de straf is voor de GAIC. In het onderstaande voorbeeld heb ik de parameter $ k = 2 $ ingesteld, wat betekent dat de ” beste ” distributie is geselecteerd volgens de klassieke AIC. U kunt $ k $ instellen op alles wat u maar wilt, zoals $ \ log (n) $ voor de 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 *** 

Volgens de AIC is de Weibull-distributie (meer specifiek WEI2, een speciale parametrisering ervan ) past het beste bij de gegevens. De exacte parametrisering van de distributie WEI2 wordt gedetailleerd beschreven in dit document op pagina 279. Laten we de aanpassing bekijken door kijken naar de residuen in een wormplot (in feite een gedateerde QQ-plot):

We verwachten dat de residuen zich dicht bij de middelste horizontale lijn bevinden en dat 95% ervan tussen de bovenste en lagere gestippelde curves, die fungeren als 95% puntsgewijze betrouwbaarheidsintervallen. In dit geval ziet de wormplot er volgens mij goed uit, wat aangeeft dat de Weibull-verdeling geschikt is.

Opmerkingen

• +1 Mooie analyse. Eén vraag echter. Geeft een positieve conclusie over compatibiliteit met een bepaalde grote distributie (in dit geval Weibull) de mogelijkheid om de mogelijkheid van een mengselverdeling uit te sluiten ‘ s aanwezigheid? Of we moeten een goede mengselanalyse uitvoeren en GoF controleren om sluit die optie uit?
• @AleksandrBlekh Het is onmogelijk om voldoende kracht te hebben om een mengsel uit te sluiten: wanneer het mengsel uit twee bijna identieke distributies bestaat, kan het niet worden gedetecteerd en wanneer alle componenten op één na zeer kleine proporties hebben het kan ook niet worden gedetecteerd. Typisch (bij gebrek aan een theorie die een distributievorm zou kunnen suggereren), past men parametrische distributies toe om spaarzame geschatte beschrijvingen van gegevens te verkrijgen. Mengsels zijn geen van deze: ze vereisen te veel parameters en zijn te flexibel voor het doel.
• @whuber: +1 Waardeer je uitstekende uitleg!
• @Lourenco Ik keek naar de Cullen en Fey-grafiek. Het blauwe punt geeft ons monster aan. Je ziet dat het punt dicht bij de lijnen van de Weibull, Lognormal en Gamma ligt (wat tussen Weibull en Gamma ligt). Nadat ik elk van deze distributies had aangepast, vergeleek ik de goodness-of-fit-statistieken met behulp van de functie gofstat en de AIC. Er is geen ‘ een consensus over wat de beste manier is om de ” beste ” distributie te bepalen is. Ik hou van grafische methoden en de AIC.
• @Lourenco Bedoel je de lognormale? De logistieke distributie (het ” + ” teken) is nogal wat verwijderd van de waargenomen gegevens. De lognormaal zou ook een kandidaat zijn waar ik ‘ normaal naar kijk. Voor deze tutorial heb ik ‘ ervoor gekozen om het niet weer te geven om het bericht kort te houden. De lognormaal laat een slechtere pasvorm zien in vergelijking met zowel de Weibull- als de normale distributie. De AIC is 537,59 en de grafieken ‘ zien er ook niet al te goed uit.

Plots zijn meestal een goede manier om een beter idee te krijgen van hoe uw gegevens eruit zien. In jouw geval zou ik aanraden om de empirische cumulatieve verdelingsfunctie (ecdf) uit te zetten tegen de theoretische cdfs met de parameters die je hebt gekregen van fitdistr ().

Ik heb dat een keer gedaan voor mijn gegevens en ook de betrouwbaarheidsintervallen opgenomen. Hier is de afbeelding die ik heb gekregen met ggplot2 ().

De zwarte lijn is de empirische cumulatieve verdelingsfunctie en de rassenbarrières zijn cdfs van verschillende distributies met behulp van parameters die ik heb gekregen met behulp van de Maximum Likelihood-methode. Men kan gemakkelijk zien dat de exponentiële en normale verdeling niet goed passen bij de data, omdat de lijnen een andere vorm hebben dan de ecdf en lijnen vrij ver van de ecdf verwijderd zijn. Helaas zijn de andere distributies redelijk dichtbij. Maar ik zou zeggen dat de lijn logNormal het dichtst bij de zwarte lijn ligt. Met behulp van een afstandsmaat (bijvoorbeeld MSE) zou je de aanname kunnen valideren.

Als je maar twee concurrerende distributies hebt (bijvoorbeeld als je degene kiest die het beste in de plot lijken te passen), kun je een Likelihood-Ratio-Test om te testen welke distributies beter passen.

Opmerkingen

• Welkom bij CrossValidated! Uw antwoord kan nuttiger zijn als u het zou kunnen bewerken om (a) de code op te nemen die u gebruikte om de afbeelding te produceren, en (b) hoe men de afbeelding zou lezen.
• Wat wordt daar geplot? Is dat een soort exponentialiteitsgrafiek?
• Maar hoe bepaal je welke distributie het beste bij je gegevens past? Alleen volgens de afbeelding kon ik ‘ niet vertellen of logNormal of weibull het beste bij uw gegevens passen.
• Als u een generator voor pseudo-willekeurige getallen wilt maken, waarom gebruik je de empirische cdf niet? Wilt u getallen trekken die verder gaan dan uw waargenomen verdeling?
• Als u uw grafiek voor de nominale waarde neemt, lijkt het erop dat geen van uw kandidaatverdelingen helemaal goed bij de gegevens past. Bovendien lijkt uw ecdf een horizontale asymptoot te hebben van minder dan 0,03 die ‘ niet logisch is, dus ik ‘ m niet zeker of het is in de eerste plaats echt een ecdf.