Hur avgör jag vilken distribution som passar mina data bäst?

> 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, här är några snabba kommentarer:

• $ p $ -värden för en Kolmovorov -Smirnov-Test (KS-Test) med uppskattade parametrar kommer att vara helt fel. Tyvärr kan du inte bara passa en distribution och sedan använda de uppskattade parametrarna i ett Kolmogorov-Smirnov-test för att testa ditt prov.
• Ditt prov kommer aldrig att följa en specifik distribution exakt. Så även om dina $ p $ -värden från KS-testet skulle vara giltiga och $ > 0,05 $ , det skulle bara betyda att du inte kan utesluta att dina data följer denna specifika distribution. En annan formulering skulle vara att ditt prov är kompatibelt med en viss distribution. Men svaret på frågan ” Följer mina uppgifter fördelningen xy exakt? ” är alltid nej.
• Målet här kan inte vara att med säkerhet bestämma vilken distribution ditt prov följer. Målet är vad @whuber (i kommentarerna) kallar parsimonious approximativ beskrivning av data. Att ha en specifik parametrisk fördelning kan vara användbar som modell för datan.

---

Men låt oss göra lite utforskning. Jag kommer att använda den utmärkta fitdistrplus paket som erbjuder några fina funktioner för distributionstillpassning. Vi använder funktionen descdist för att få några idéer om möjliga kandidatfördelningar.

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) 

Nu kan vi använda descdist:

descdist(x, discrete = FALSE) 

Provets kurtos och kvadratiska snedhet är plottet som en blå punkt med namnet ” Observation ”. Det verkar som om möjliga fördelningar inkluderar Weibull, Lognormal och möjligen gammafördelningen.

Låt oss passa en Weibull-distribution och en normalfördelning:

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

Inspektera nu anpassningen för normal:

plot(fit.norm) 

Och för Weibull-fit:

plot(fit.weibull) 

Båda ser bra ut men bedöms av QQ-plot, Weibull kanske ser lite bättre ut, särskilt i svansarna. På motsvarande sätt är AIC för Weibull-passformen lägre jämfört med normal passform:

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

---

Kolmogorov-Smirnov testsimulering

Jag kommer att använda @Aksakals procedur förklarad här för att simulera KS-statistiken under noll.

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 för den simulerade KS-statistiken ser ut enligt följande:

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

Slutligen, vår $ p $ -värde med den simulerade nollfördelningen av KS-statistiken är:

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

Detta bekräftar vår grafiska slutsats att provet är kompatibelt med en Weibull-fördelning.

Som förklaras här kan vi använda bootstrapping för att lägga till punktvis konfidensintervall till den uppskattade 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 distributionspassning med GAMLSS

gamlss paket för R erbjuder möjligheten att prova många olika distributioner och välja bästa ” enligt GAIC (det allmänna Akaike-informationskriteriet). Huvudfunktionen är fitDist. Ett viktigt alternativ i den här funktionen är typen av distributioner som provas. Att till exempel ställa in type = "realline" försöker alla implementerade distributioner definierade på hela den verkliga raden medan type = "realsplus" bara försöker distributioner definierade på den verkliga positiva raden . Ett annat viktigt alternativ är parametern $ k $ , vilket är straffet för GAIC. I exemplet nedan ställer jag in parametern $ k = 2 $ vilket innebär att ” bäst ” distribution väljs enligt den klassiska AIC. Du kan ställa in $ k $ till vad du vill, till exempel $ \ log (n) $ för 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 *** 

Enligt AIC är Weibull-fördelningen (mer specifikt WEI2, en speciell parametrisering av den ) passar data bäst. Den exakta parametriseringen av distributionen WEI2 detailleras i detta dokument på sidan 279. Låt oss inspektera passformen genom att tittar på resterna i en maskritning (i grund och botten en utvecklad QQ-plot):

Vi förväntar oss att resterna ligger nära den mittre horisontella linjen och 95% av dem ligger mellan den övre och nedre prickade kurvor, som fungerar som 95% punktvis konfidensintervall. I det här fallet ser maskritningen mig bra ut, vilket indikerar att Weibull-fördelningen passar.

Kommentarer

+1 Trevlig analys. En fråga. Men gör positiv slutsats om kompatibilitet med en viss större fördelning (Weibull, i det här fallet) att utesluta en möjlighet till en blandningsfördelning ’ närvaro? Eller vi måste utföra en ordentlig blandningsanalys och kontrollera GoF till utesluta det alternativet?

@AleksandrBlekh Det är omöjligt att ha tillräckligt med kraft för att utesluta en blandning: när blandningen har två nästan identiska fördelningar kan det inte detekteras och när alla utom en komponent har mycket små proportioner det kan inte heller detekteras. Vanligtvis (i avsaknad av en teori som kan föreslå en distributionsform) passar man parametriska fördelningar för att uppnå parsimonious ungefärliga beskrivningar av data. Blandningar är ingen av dessa: de kräver för många parametrar och är för flexibla för ändamålet.

@whuber: +1 Uppskattar din utmärkta förklaring!

@Lourenco Jag tittade på Cullen och Fey-diagrammet. Den blå punkten betecknar vårt urval. Du ser att punkten ligger nära linjerna i Weibull, Lognormal och Gamma (som ligger mellan Weibull och Gamma). Efter att ha anpassat var och en av dessa distributioner jämförde jag statistik över godhet med hjälp av funktionen gofstat och AIC. Det finns ’ en konsensus om vad det bästa sättet att bestämma ” bästa ” distribution är. Jag gillar grafiska metoder och AIC.

@Lourenco Menar du det lognormala? Den logistiska fördelningen (tecknet ” + ”) ligger ganska långt ifrån de observerade uppgifterna. Det lognormala skulle också vara en kandidat som jag ’ normalt tittar på. För den här självstudien har jag ’ valt att inte visa det för att hålla inlägget kort. Det lognormala visar en sämre passform jämfört med både Weibull och Normal distribution. AIC är 537,59 och graferna ser inte heller ’ ut.

Plottar är oftast ett bra sätt att få en bättre uppfattning om hur dina data ser ut. I ditt fall skulle jag rekommendera att plotta empirisk kumulativ fördelningsfunktion (ecdf) mot de teoretiska cdfs med de parametrar du fick från fitdistr ().

Jag gjorde det en gång för mina uppgifter och inkluderade också konfidensintervall. Här är bilden som jag fick med ggplot2 ().

Den svarta linjen är den empiriska kumulativa fördelningsfunktionen och de färgade linjerna är cdfs från olika distributioner med hjälp av parametrar som jag fick med metoden Maximum Likelihood. Man kan lätt se att den exponentiella och normala fördelningen inte passar bra till data, eftersom linjerna har en annan form än ecdf och linjer är ganska långt ifrån ecdf. Tyvärr är de andra fördelningarna ganska nära. Men jag skulle säga att logNormal-linjen är närmast den svarta linjen. Med hjälp av ett mått på avstånd (till exempel MSE) kan man validera antagandet.

Om du bara har två konkurrerande fördelningar (till exempel att välja de som verkar passa bäst i plottet) kan du använda en Likelihood-Ratio-Test för att testa vilka distributioner som passar bättre.

Kommentarer

• Välkommen till CrossValidated! Ditt svar kan vara mer användbart om du kan redigera det för att inkludera (a) koden du använde för att producera grafiken och (b) hur man skulle läsa bilden.
• Vad planeras där? Är det någon form av exponentiality-plot?
• Men hur bestämmer du vilken distribution som passar dina data bäst? Endast enligt bilden kunde jag inte ’ inte berätta om logNormal eller weibull passar bäst i dina data.
• Om du vill skapa en pseudo-slumpmässig talgenerator varför inte använda den empiriska cdf? Vill du rita siffror som går utöver din observerade fördelning?
• Om du tar grafen till nominellt värde verkar det som om ingen av dina kandidatfördelningar passar data alls. Din ecdf verkar också ha en horisontell asymptot på mindre än 0,03 vilket inte ’ inte är vettigt, så jag ’ är inte säker på att det är verkligen en ecdf i första hand.