Hvordan finder jeg ud af, hvilken distribution der passer bedst til mine data?

> 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 

For det første er der nogle hurtige kommentarer:

• $ p $ -værdierne for en Kolmovorov -Smirnov-Test (KS-Test) med estimerede parametre vil være ret forkert. Så desværre kan du ikke bare tilpasse en distribution og derefter bruge de estimerede parametre i en Kolmogorov-Smirnov-test til at teste din prøve.
• Din prøve vil aldrig følge en bestemt distribution nøjagtigt. Så selvom din $ p $ -værdier fra KS-testen ville være gyldig og $ > 0,05 $ , det ville bare betyde, at du ikke kan udelukke at dine data følger denne specifikke distribution. En anden formulering ville være, at din prøve er kompatibel med en vis fordeling. Men svaret på spørgsmålet ” Følger mine data fordelingen xy nøjagtigt? ” er altid nej.
• Målet her kan ikke være med sikkerhed at bestemme, hvilken distribution din prøve følger. Målet er, hvad @whuber (i kommentarerne) kalder parsimonious omtrentlige beskrivelser af dataene. At have en bestemt parametrisk fordeling kan være nyttig som model for dataene.

---

Men lad os udforske noget. Jeg vil bruge den fremragende fitdistrplus pakke, der tilbyder nogle gode funktioner til distributionstilpasning. Vi bruger funktionen descdist nogle ideer om mulige kandidatfordelinger.

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) 

Lad os nu bruge descdist:

descdist(x, discrete = FALSE) 

Kurtosen og den kvadratiske skævhed i din prøve er plottet som et blåt punkt med navnet ” Observation “. Det ser ud til, at mulige fordelinger inkluderer Weibull, Lognormal og muligvis gammafordeling.

Lad os passe til Weibull-distribution og en normalfordeling:

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

Undersøg nu pasformen til normal:

plot(fit.norm) 

Og for Weibull-fit:

plot(fit.weibull) 

Begge ser godt ud, men bedømt af QQ-plottet, ser Weibull måske lidt bedre ud, især i halerne. Tilsvarende er AIC for Weibull-pasformen lavere sammenlignet med den normale pasform:

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

---

Kolmogorov-Smirnov testsimulering

Jeg vil bruge @Aksakals procedure forklaret her til at simulere KS-statistikken under null.

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

ECDF for den simulerede KS-statistik ser således ud:

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

Endelig, vores $ p $ -værdi ved hjælp af den simulerede nulfordeling af KS-statistikken er:

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

Dette bekræfter vores grafiske konklusion, at prøven er kompatibel med en Weibull-distribution.

Som forklaret her kan vi bruge bootstrapping til at tilføje punktvise konfidensintervaller til den estimerede Weibull PDF eller 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") 

---

Automatisk fordelingstilpasning med GAMLSS

gamlss pakke til R giver mulighed for at prøve mange forskellige distributioner og vælge bedste ” i henhold til GAIC (det generelle Akaike-informationskriterium). Hovedfunktionen er fitDist. En vigtig mulighed i denne funktion er typen af distributioner, der er prøvet. For eksempel vil indstilling af type = "realline" prøve alle implementerede distributioner defineret på hele den reelle linje, mens type = "realsplus" kun vil prøve distributioner defineret på den reelle positive linje . En anden vigtig mulighed er parameteren $ k $ , hvilket er straffen for GAIC. I eksemplet nedenfor indstiller jeg parameteren $ k = 2 $ , hvilket betyder, at ” bedste ” distribution er valgt i henhold til den klassiske AIC. Du kan indstille $ k $ til alt hvad du kan lide, såsom $ \ log (n) $ for 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 *** 

Ifølge AIC er Weibull-fordelingen (mere specifikt WEI2, en særlig parametrisering af den ) passer bedst til dataene. Den nøjagtige parametrisering af fordelingen WEI2 er detaljeret i dette dokument på side 279. Lad os inspicere pasformen ved at ser på resterne i et ormplot (dybest set et udrettet QQ-plot):

Vi forventer, at resterne ligger tæt på den midterste vandrette linje og 95% af dem ligger mellem den øverste og nedre punkterede kurver, som fungerer som 95% punktvise konfidensintervaller. I dette tilfælde ser ormplottet fint ud for mig, hvilket indikerer, at Weibull-fordelingen er passende.

Kommentarer

• +1 Flot analyse. Et spørgsmål, dog. Gør positiv konklusion om kompatibilitet med en bestemt større distribution (Weibull, i dette tilfælde) det muligt at udelukke en mulighed for en blandingsfordeling ‘ s tilstedeværelse? Eller vi skal udføre en ordentlig blandingsanalyse og kontrollere GoF til udelukke denne mulighed?
• @AleksandrBlekh Det er umuligt at have nok kraft til at udelukke en blanding: når blandingen har to næsten identiske fordelinger, kan den ikke detekteres, og når alle undtagen en komponent har meget små proportioner det kan heller ikke detekteres. Typisk (i mangel af en teori, der kan antyde en distributionsform), passer man til parametriske fordelinger for at opnå parsimonious omtrentlige beskrivelser af data. Blandinger er ingen af dem: de kræver for mange parametre og er for fleksible til formålet.
• @whuber: +1 Værdsat din fremragende forklaring!
• @Lourenco Jeg kiggede på Cullen og Fey grafen. Det blå punkt angiver vores prøve. Du ser, at punktet er tæt på linjerne i Weibull, Lognormal og Gamma (som er mellem Weibull og Gamma). Efter tilpasning af hver af disse distributioner sammenlignede jeg statistikken over godhed ved hjælp af funktionen gofstat og AIC. Der er ikke ‘ en konsensus om, hvad den bedste måde at bestemme ” bedste ” fordeling på er. Jeg kan godt lide grafiske metoder og AIC.
• @Lourenco Mener du det lognormale? Den logistiske fordeling (” + ” -tegnet) er ret lidt væk fra de observerede data. Det lognormale ville også være en kandidat, jeg ‘ normalt ser på. Til denne vejledning har jeg ‘ valgt at ikke vise det for at holde posten kort. Det lognormale viser en dårligere pasform sammenlignet med både Weibull- og Normalfordelingen. AIC er 537,59, og graferne ser heller ikke ‘ ud for godt.

Plotter er for det meste en god måde at få en bedre idé om, hvordan dine data ser ud. I dit tilfælde vil jeg anbefale at plotte empirisk kumulativ fordelingsfunktion (ecdf) mod de teoretiske cdfs med de parametre, du fik fra fitdistr ().

Det gjorde jeg en gang for mine data og inkluderede også tillidsintervaller. Her er det billede, jeg fik ved hjælp af ggplot2 ().

Den sorte linje er den empiriske kumulative fordelingsfunktion og de farvede linjer er cdfs fra forskellige distributioner ved hjælp af parametre, jeg fik ved hjælp af metoden Maximum Likelihood. Man kan let se, at den eksponentielle og normale fordeling ikke passer godt til dataene, fordi linjerne har en anden form end ecdf, og linjer er ret langt væk fra ecdf. Desværre er de andre fordele ret tætte. Men jeg vil sige, at logNormal-linjen er tættest på den sorte linje. Ved hjælp af et mål for afstand (for eksempel MSE) kunne man validere antagelsen.

Hvis du kun har to konkurrerende distributioner (for eksempel at vælge dem, der synes at passe bedst i plottet), kan du bruge en Likelihood-Ratio-Test for at teste, hvilke distributioner der passer bedre.

Kommentarer

• Velkommen til CrossValidated! Dit svar kan være mere nyttigt, hvis du kunne redigere det til at inkludere (a) den kode, du brugte til at producere grafikken, og (b) hvordan man ville læse grafikken.
• Hvad er der plottet der? Er det en slags eksponentialitetsplot?
• Men hvordan beslutter du, hvilken distribution der passer bedst til dine data? Kun ifølge grafikken kunne jeg ikke ‘ ikke fortælle dig, om logNormal eller weibull passer bedst til dine data.
• Hvis du vil oprette en pseudo-tilfældig talgenerator, hvorfor ikke bruge den empiriske cdf? Vil du tegne tal, der går ud over din observerede fordeling?
• Når du tager din graf til pålydende værdi, ser det ud til, at ingen af dine kandidatfordelinger overhovedet passer til dataene. Også din ecdf ser ud til at have en vandret asymptote på mindre end 0,03, som ikke ‘ ikke giver mening, så jeg ‘ er ikke sikker på, at det er virkelig en ecdf i første omgang.