Hvordan bestemme hvilken distribusjon som passer best til dataene mine?

> 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 

Først er det noen raske kommentarer:

• $ p $ -verdiene til en Kolmovorov -Smirnov-Test (KS-Test) med estimerte parametere vil være ganske feil. Så dessverre kan du ikke bare passe på en fordeling og deretter bruke de estimerte parametrene i en Kolmogorov-Smirnov-test for å teste prøven din.
• Eksemplet ditt vil aldri følge en spesifikk distribusjon nøyaktig. Så selv om $ p $ -verdiene fra KS-testen ville være gyldige og $ > 0,05 $ , vil det bare bety at du ikke kan utelukke at dataene dine følger denne spesifikke distribusjonen. En annen formulering vil være at prøven din er kompatibel med en viss fordeling. Men svaret på spørsmålet » Følger dataene mine fordelingen xy nøyaktig? » er alltid nei.
• Målet her kan ikke være å fastslå med sikkerhet hvilken fordeling prøven din følger. Målet er det @whuber (i kommentarene) kaller parsimonious omtrentlige beskrivelser av dataene. Å ha en spesifikk parametrisk fordeling kan være nyttig som en modell for dataene.

---

Men la oss utforske litt. Jeg vil bruke den utmerkede fitdistrplus pakke som tilbyr noen fine funksjoner for distribusjonstilpasning. Vi vil bruke funksjonen descdist for å få noen 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) 

Nå kan vi bruke descdist:

descdist(x, discrete = FALSE) 

Kurtosen og den kvadratiske skjevheten til prøven din er plottet som et blått punkt med navnet » Observasjon «. Det ser ut til at mulige distribusjoner inkluderer Weibull, Lognormal og muligens gammafordeling.

La oss passe til Weibull-fordeling og en normalfordeling:

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

Kontroller nå om det passer til det normale:

plot(fit.norm) 

Og for Weibull-fit:

plot(fit.weibull) 

Begge ser bra ut, men bedømt av QQ-plottet, ser Weibull kanskje litt bedre ut, spesielt i halene. Tilsvarende er AIC for Weibull-passformen lavere sammenlignet med normal passform:

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

---

Kolmogorov-Smirnov testsimulering

Jeg vil bruke @Aksakals prosedyre forklart her for å 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 simulerte KS-statistikken ser slik ut:

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

Endelig, vår $ p $ -verdi ved hjelp av den simulerte nullfordelingen av 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 bekrefter vår grafiske konklusjon at prøven er kompatibel med en Weibull-fordeling.

Som forklart her , kan vi bruke bootstrapping for å legge til punktvise konfidensintervaller til den estimerte 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 distribusjonstilpasning med GAMLSS

gamlss pakke for R tilbyr muligheten til å prøve mange forskjellige distribusjoner og velge best » i henhold til GAIC (det generelle Akaike-informasjonskriteriet). Hovedfunksjonen er fitDist. Et viktig alternativ i denne funksjonen er typen distribusjoner som blir prøvd. For eksempel vil innstilling av type = "realline" prøve alle implementerte distribusjoner definert på hele den virkelige linjen, mens type = "realsplus" bare vil prøve distribusjoner definert på den virkelige positive linjen . Et annet viktig alternativ er parameteren $ k $ , som er straffen for GAIC. I eksemplet nedenfor setter jeg parameteren $ k = 2 $ som betyr at » best » distribusjon er valgt i henhold til den klassiske AIC. Du kan angi $ k $ til hva du vil, for eksempel $ \ 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 *** 

I følge AIC er Weibull-fordelingen (nærmere bestemt WEI2, en spesiell parametrisering av den ) passer best til dataene. Den nøyaktige parametriseringen av fordelingen WEI2 er detaljert i dette dokumentet på side 279. La oss inspisere passformen ved å ser på restene i et ormplott (i utgangspunktet et avviklet QQ-plot):

Vi forventer at restene ligger nær den midtre horisontale linjen og 95% av dem ligger mellom den øvre og nedre prikkede kurver, som fungerer som 95% konfidensintervaller i punkt. I dette tilfellet ser ormplottet bra ut for meg, noe som indikerer at Weibull-fordelingen passer tilstrekkelig.

Kommentarer

• +1 Fin analyse. Et spørsmål, skjønt. Har positiv konklusjon om kompatibilitet med en bestemt hovedfordeling (Weibull, i dette tilfellet) tillatelse til å utelukke muligheten for en blandingsfordeling ‘ s tilstedeværelse? Eller vi må utføre en ordentlig blandingsanalyse og sjekke GoF til utelukke det alternativet?
• @AleksandrBlekh Det er umulig å ha nok kraft til å utelukke en blanding: når blandingen har to nesten identiske fordelinger, kan den ikke oppdages, og når alle bortsett fra en komponent har veldig små proporsjoner det kan heller ikke oppdages. Vanligvis (i fravær av en teori som kan antyde en distribusjonsform), passer man til parametriske fordelinger for å oppnå parsimonious omtrentlige beskrivelser av data. Blandinger er ingen av disse: de krever for mange parametere og er for fleksible for formålet.
• @whuber: +1 Setter pris på din utmerkede forklaring!
• @Lourenco Jeg så på Cullen og Fey-grafen. Det blå punktet betegner prøven vår. Du ser at poenget er nær linjene til Weibull, Lognormal og Gamma (som er mellom Weibull og Gamma). Etter å ha tilpasset hver av disse distribusjonene, sammenlignet jeg godhetsstatistikken ved hjelp av funksjonen gofstat og AIC. Det er ikke ‘ en konsensus om hva den beste måten å bestemme » beste » er. Jeg liker grafiske metoder og AIC.
• @Lourenco Mener du det lognormale? Den logistiske fordelingen (» + » -tegnet) er ganske mye borte fra de observerte dataene. Det lognormale ville også være en kandidat jeg ‘ normalt ser på. For denne veiledningen har jeg ‘ valgt å ikke vise den for å holde innlegget kort. Det lognormale viser en dårligere passform sammenlignet med både Weibull- og Normal-fordelingen. AIC er 537,59, og grafene ser ikke ‘ t ut for bra.

Plott er for det meste en god måte å få et bedre inntrykk av hvordan dataene dine ser ut. I ditt tilfelle vil jeg anbefale å plotte empirisk kumulativ distribusjonsfunksjon (ecdf) mot de teoretiske cdfs med parametrene du fikk fra fitdistr ().

Jeg gjorde det en gang for dataene mine og inkluderte også konfidensintervallene. Her er bildet jeg fikk med ggplot2 ().

Den svarte linjen er den empiriske kumulative fordelingsfunksjonen og de fargede linjene er cdfs fra forskjellige distribusjoner ved hjelp av parametere jeg fikk ved hjelp av Maximum Likelihood-metoden. Man kan lett se at den eksponensielle og normalfordelingen ikke passer godt til dataene, fordi linjene har en annen form enn ecdf og linjer er ganske langt borte fra ecdf. Dessverre er de andre distribusjonene ganske tette. Men jeg vil si at logNormal-linjen er nærmest den svarte linjen. Ved å bruke et mål på avstand (for eksempel MSE) kan man validere antagelsen.

Hvis du bare har to konkurrerende distribusjoner (for eksempel å velge de som ser ut til å passe best i plottet), kan du bruke en Likelihood-Ratio-Test for å teste hvilke distribusjoner som passer bedre.

Kommentarer

• Velkommen til CrossValidated! Svaret ditt kan være mer nyttig hvis du kan redigere det for å inkludere (a) koden du brukte til å produsere grafikken, og (b) hvordan man ville lese grafikken.
• Hva blir tegnet der? Er det en slags eksponensialitetsplott?
• Men hvordan bestemmer du hvilken distribusjon som passer best til dine data? Bare i henhold til grafikken kunne jeg ikke ‘ ikke fortelle deg om logNormal eller weibull passer best til dine data.
• Hvis du vil lage en pseudo-tilfeldig tallgenerator, hvorfor ikke bruke empirisk cdf? Vil du tegne tall som går utover den observerte fordelingen din?
• Hvis du tar grafen til pålydende, ser det ut til at ingen av kandidatfordelingene dine passer i det hele tatt. Også, ecdf ser ut til å ha en horisontal asymptote på mindre enn 0,03 som ikke ‘ ikke gir mening, så jeg ‘ er ikke sikker på at det er virkelig en ecdf i utgangspunktet.