Wie kann ich feststellen, welche Verteilung am besten zu meinen Daten passt?

> 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 

Zunächst einige kurze Kommentare:

• Die $ p $ -Werte eines Kolmovorov -Smirnov-Test (KS-Test) mit geschätzten Parametern ist völlig falsch. Leider können Sie nicht einfach eine Verteilung anpassen und dann die geschätzten Parameter in einem Kolmogorov-Smirnov-Test verwenden, um Ihre Probe zu testen.
• Ihre Probe wird niemals einer bestimmten folgen Verteilung genau. Selbst wenn Ihre $ p $ -Werte aus dem KS-Test gültig wären und $ > 0,05 $ würde nur bedeuten, dass Sie nicht ausschließen können , dass Ihre Daten dieser spezifischen Verteilung folgen. Eine andere Formulierung wäre, dass Ihre Probe mit einer bestimmten Verteilung kompatibel ist. Die Antwort auf die Frage “ Entspricht meine Daten genau der Verteilung xy? “ lautet immer no.
• Das Ziel hier kann nicht sein, mit Sicherheit zu bestimmen, welcher Verteilung Ihre Stichprobe folgt. Das Ziel ist, was @whuber (in den Kommentaren) sparsame ungefähre Beschreibungen der Daten nennt. Eine bestimmte parametrische Verteilung kann als Modell für die Daten nützlich sein.

---

Aber lassen Sie uns etwas untersuchen. Ich werde die ausgezeichnete fitdistrplus Paket, das einige nette Funktionen für die Verteilungsanpassung bietet. Wir werden die Funktion descdist verwenden, um zu gewinnen Einige Ideen zu möglichen Kandidatenverteilungen.

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) 

Verwenden wir jetzt descdist:

descdist(x, discrete = FALSE) 

Die Kurtosis und die quadratische Schiefe Ihrer Probe werden als blauer Punkt mit dem Namen “ Beobachtung „. Es scheint, dass mögliche Verteilungen die Weibull-, Lognormal- und möglicherweise die Gamma-Verteilung umfassen.

Lassen Sie uns a passen Weibull-Verteilung und Normalverteilung:

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

Überprüfen Sie nun die Passform für die Normalverteilung:

plot(fit.norm) 

Und für die Weibull-Anpassung:

plot(fit.weibull) 

Beide sehen gut aus, aber nach dem QQ-Plot sieht der Weibull vielleicht etwas besser aus, besonders an den Schwänzen. Entsprechend ist der AIC der Weibull-Anpassung im Vergleich zur normalen Anpassung niedriger:

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

---

Kolmogorov-Smirnov-Testsimulation

Ich werde @Aksakals Prozedur verwenden, die hier erklärt wurde, um die KS-Statistik unter der Null zu simulieren.

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

Das ECDF der simulierten KS-Statistik sieht folgendermaßen aus:

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

Schließlich unser $ p $ -Wert unter Verwendung der simulierten Nullverteilung der KS-Statistik lautet:

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

Dies bestätigt unsere grafische Schlussfolgerung, dass die Stichprobe mit einer Weibull-Verteilung kompatibel ist.

Wie erläutert hier können wir Bootstrapping verwenden, um dem geschätzten Weibull-PDF oder CDF punktweise Konfidenzintervalle hinzuzufügen:

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 Verteilungsanpassung mit GAMLSS

Die gamlss Paket für R bietet die Möglichkeit, viele verschiedene Distributionen auszuprobieren und best “ gemäß GAIC (dem verallgemeinerten Akaike-Informationskriterium). Die Hauptfunktion ist fitDist. Eine wichtige Option in dieser Funktion ist der Typ der Distributionen, die ausprobiert werden. Wenn Sie beispielsweise type = "realline" festlegen, werden alle implementierten Verteilungen ausprobiert, die in der gesamten realen Zeile definiert sind, während type = "realsplus" nur Verteilungen versucht, die in der realen positiven Zeile definiert sind . Eine weitere wichtige Option ist der Parameter $ k $ , der die Strafe für den GAIC darstellt. Im folgenden Beispiel habe ich den Parameter $ k = 2 $ gesetzt, was bedeutet, dass die “ beste wird gemäß dem klassischen AIC ausgewählt. Sie können $ k $ auf einen beliebigen Wert setzen, z. B. $ \ 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 *** 

Laut AIC die Weibull-Verteilung (genauer gesagt WEI2, eine spezielle Parametrisierung davon ) passt am besten zu den Daten. Die genaue Parametrisierung der Verteilung WEI2 wird in dieses Dokuments auf Seite 279 detailliert beschrieben. Lassen Sie uns die Anpassung durch überprüfen Betrachten der Residuen in einem -Wurmplot (im Grunde ein de-trendierter QQ-Plot):

Wir erwarten, dass die Residuen nahe an der mittleren horizontalen Linie liegen und 95% von ihnen zwischen der oberen liegen und niedrigere gepunktete Kurven, die als 95% punktweise Konfidenzintervalle dienen. In diesem Fall sieht das Wurmdiagramm für mich gut aus, was darauf hinweist, dass die Weibull-Verteilung eine angemessene Anpassung aufweist.

Kommentare

• +1 Gute Analyse. Eine Frage jedoch. Kann eine positive Schlussfolgerung zur Kompatibilität mit einer bestimmten Hauptverteilung (in diesem Fall Weibull) die Möglichkeit einer Mischungsverteilung ausschließen ‚ Anwesenheit? Oder wir müssen eine ordnungsgemäße Gemischanalyse durchführen und GoF auf überprüfen Diese Option ausschließen?
• @AleksandrBlekh Es ist unmöglich, genügend Leistung zu haben, um eine Mischung auszuschließen: Wenn die Mischung zwei nahezu identische Verteilungen aufweist, kann sie nicht erkannt werden, und wenn alle bis auf eine Komponente sehr kleine Anteile haben es kann auch nicht erkannt werden. Typischerweise (in Abwesenheit einer Theorie, die eine Verteilungsform vorschlagen könnte) passt man parametrische Verteilungen an, um sparsame ungefähre Beschreibungen von Daten zu erhalten. Mischungen sind keine davon: Sie erfordern zu viele Parameter und sind für diesen Zweck zu flexibel.
• @whuber: +1 Schätzen Sie Ihre ausgezeichnete Erklärung!
• @Lourenco Ich habe mir das Diagramm von Cullen und Fey angesehen. Der blaue Punkt kennzeichnet unsere Probe. Sie sehen, dass der Punkt nahe an den Linien von Weibull, Lognormal und Gamma liegt (zwischen Weibull und Gamma). Nachdem ich jede dieser Verteilungen angepasst hatte, verglich ich die Anpassungsgütestatistik mit der Funktion gofstat und dem AIC. Es gibt keinen ‚ Konsens darüber, wie die “ beste “ -Verteilung am besten bestimmt werden kann ist. Ich mag grafische Methoden und den AIC.
• @Lourenco Meinst du das Lognormal? Die logistische Verteilung (das Zeichen “ + „) ist ziemlich weit von den beobachteten Daten entfernt. Das Lognormal wäre auch ein Kandidat, den ich ‚ normalerweise betrachten würde. Für dieses Tutorial habe ich ‚ ausgewählt, es nicht anzuzeigen, um den Beitrag kurz zu halten. Das Lognormal zeigt eine schlechtere Anpassung im Vergleich zur Weibull- und Normalverteilung. Der AIC ist 537,59 und die Grafiken sehen auch ‚ nicht allzu gut aus.

Diagramme sind meistens eine gute Möglichkeit, um eine bessere Vorstellung davon zu bekommen, wie Ihre Daten aussehen. In Ihrem Fall würde ich empfehlen, die empirische kumulative Verteilungsfunktion (ecdf) gegen die theoretischen cdfs mit den von fitdistr erhaltenen Parametern zu zeichnen ().

Ich habe das einmal für meine Daten gemacht und auch die Konfidenzintervalle eingeschlossen. Hier ist das Bild, das ich mit ggplot2 () erhalten habe.

Die schwarze Linie ist die empirische kumulative Verteilungsfunktion und die farbigen Linien sind cdfs aus verschiedenen Verteilungen unter Verwendung von Parametern, die ich mit der Maximum Likelihood-Methode erhalten habe. Man kann leicht erkennen, dass die Exponential- und Normalverteilung nicht gut zu den Daten passt, da die Linien eine andere Form als das ecdf haben und die Linien ziemlich weit vom ecdf entfernt sind. Leider sind die anderen Verteilungen ziemlich nah. Aber ich würde sagen, dass die logNormal-Linie der schwarzen Linie am nächsten kommt. Mit einem Abstandsmaß (z. B. MSE) könnte die Annahme validiert werden.

Wenn Sie nur zwei konkurrierende Verteilungen haben (z. B. diejenigen auswählen, die am besten in die Darstellung passen), können Sie ein Likelihood-Ratio-Test , um zu testen, welche Verteilungen besser passen.

Kommentare

• Willkommen bei CrossValidated! Ihre Antwort könnte nützlicher sein, wenn Sie sie so bearbeiten könnten, dass sie (a) den Code enthält, mit dem Sie die Grafik erstellt haben, und (b) wie man die Grafik lesen würde.
• Was wird dort gezeichnet? Ist das eine Art Exponentialitätsdiagramm?
• Aber wie entscheiden Sie, welche Verteilung am besten zu Ihren Daten passt? Nur gemäß der Grafik konnte ich ‚ nicht sagen, ob logNormal oder weibull am besten zu Ihren Daten passt.
• Wenn Sie einen Pseudozufallszahlengenerator erstellen möchten, warum nicht das empirische cdf verwenden? Möchten Sie Zahlen zeichnen, die über Ihre beobachtete Verteilung hinausgehen?
• Wenn Sie Ihr Diagramm zum Nennwert nehmen, scheint es, dass keine Ihrer Kandidatenverteilungen überhaupt gut zu den Daten passt. Außerdem scheint Ihr ecdf eine horizontale Asymptote von weniger als 0,03 zu haben, was ‚ keinen Sinn ergibt, daher bin ich mir nicht sicher, ob ‚ Es ist wirklich in erster Linie ein Ecdf.