Wie bestimme ich, welche Distribution am besten zu meinen Daten passt?

133

Ich habe einen Datensatz und möchte herausfinden, welche Verteilung am besten zu meinen Daten passt.

Ich habe die fitdistr()Funktion verwendet, um die notwendigen Parameter zur Beschreibung der angenommenen Verteilung abzuschätzen (z. B. Weibull, Cauchy, Normal). Mit diesen Parametern kann ich einen Kolmogorov-Smirnov-Test durchführen, um abzuschätzen, ob meine Probendaten aus derselben Verteilung stammen wie meine angenommene Verteilung.

Wenn der p-Wert> 0,05 ist, kann ich davon ausgehen, dass die Probendaten aus derselben Verteilung stammen. Der p-Wert gibt aber keine Auskunft über die Anpassungsgöttlichkeit, oder?

Wenn der p-Wert meiner Beispieldaten für eine Normalverteilung und eine Weibullverteilung> 0,05 ist, wie kann ich dann herausfinden, welche Verteilung besser zu meinen Daten passt?

Dies ist im Grunde das, was ich getan habe:

> 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

Die p-Werte betragen 0,8669 für die Weibull-Verteilung und 0,5522 für die Normalverteilung. Somit kann ich davon ausgehen, dass meine Daten sowohl einer Weibull- als auch einer Normalverteilung folgen. Aber welche Verteilungsfunktion beschreibt meine Daten besser?


In Bezug auf elevendollar habe ich den folgenden Code gefunden, weiß aber nicht, wie ich die Ergebnisse interpretieren soll:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
Tobibo
quelle
5
Warum möchten Sie herausfinden, welche Distribution am besten zu Ihren Daten passt?
Roland
6
Weil ich Pseudozufallszahlen generieren möchte, die der gegebenen Verteilung folgen.
Tobibo
6
Sie können KS nicht verwenden, um zu überprüfen, ob eine Verteilung mit Parametern, die aus dem Dataset gefunden wurden, mit dem Dataset übereinstimmt. Siehe Nr. 2 auf dieser Seite sowie Alternativen (und andere Möglichkeiten, wie der KS-Test irreführend sein kann).
tpg2114
Eine weitere Diskussion hier mit Codebeispielen zur Anwendung des KS-Tests, wenn Parameter aus dem Beispiel geschätzt werden.
Aksakal
1
I used the fitdistr() function ..... Was ist fitdistrFunktion? Etwas aus Excel? Oder etwas, das Sie selbst in C geschrieben haben?
Wolfies

Antworten:

162

Hier sind zunächst einige kurze Kommentare:

  • Die Werte eines Kolmovorov-Smirnov-Tests (KS-Test) mit geschätzten Parametern sind völlig falsch. Sie können also leider nicht einfach eine Verteilung anpassen und dann die geschätzten Parameter in einem Kolmogorov-Smirnov-Test verwenden, um Ihre Probe zu testen.p
  • Ihre Stichprobe folgt niemals genau einer bestimmten Verteilung. Selbst wenn Ihre Werte aus dem KS-Test gültig wären und , würde dies 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. Aber die Antwort auf die Frage "Entsprechen meine Daten genau der Verteilung xy?" ist immer nein.p>0.05
  • Das Ziel kann hier nicht sein, mit Sicherheit zu bestimmen, welcher Verteilung Ihre Stichprobe folgt. Das Ziel ist das, was @whuber (in den Kommentaren) als sparsame ungefähre Beschreibung der Daten bezeichnet. Eine bestimmte parametrische Verteilung kann als Modell der Daten nützlich sein.

Aber lasst uns etwas erforschen. Ich werde das exzellente fitdistrplusPaket verwenden, das einige nette Funktionen für die Verteilungsanpassung bietet. Wir werden die Funktion nutzen descdist, um einige Ideen über mögliche Verteilungskandidaten zu erhalten.

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)

Jetzt können wir verwenden descdist:

descdist(x, discrete = FALSE)

Descdist

Die Kurtosis und die quadratische Schiefe Ihrer Probe werden als blauer Punkt mit der Bezeichnung "Beobachtung" aufgezeichnet. Es scheint, dass mögliche Verteilungen die Weibull-, Lognormal- und möglicherweise die Gamma-Verteilung umfassen.

Passen wir eine Weibull-Verteilung und eine Normalverteilung an:

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

Überprüfen Sie nun die Passform für den Normalfall:

plot(fit.norm)

Normale Passform

Und für die Weibull-Passform:

plot(fit.weibull)

Weibull fit

Beide sehen gut aus, aber gemessen am QQ-Plot sieht der Weibull vielleicht ein bisschen besser aus, besonders an den Schwänzen. Dementsprechend 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 die hier erläuterte Prozedur von @ Aksakal verwenden , 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()

Simulierte KS-Statistik

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

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 Probe mit einer Weibull-Verteilung kompatibel ist.

Wie hier erläutert , 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)

CI_Dichte

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

CI_CDF


Automatische Verteilungsanpassung mit GAMLSS

Das gamlssPaket für Rbietet die Möglichkeit, viele verschiedene Distributionen auszuprobieren und die "besten" gemäß GAIC (dem verallgemeinerten Akaike-Informationskriterium) auszuwählen. Die Hauptfunktion ist fitDist. Eine wichtige Option in dieser Funktion ist die Art der ausprobierten Verteilungen. Zum Beispiel werden bei der Einstellung type = "realline"alle implementierten Verteilungen ausprobiert, die für die gesamte reale Linie definiert sind, wohingegen type = "realsplus"nur Verteilungen ausprobiert werden, die für die reale positive Linie definiert sind. Eine weitere wichtige Option ist der Parameter , der die Strafe für den GAIC darstellt. Im folgenden Beispiel habe ich den Parameter was bedeutet, dass die "beste" Verteilung gemäß dem klassischen AIC ausgewählt wird. Sie können auf einen beliebigen einstellen , zkk=2klog(n) für den 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 WEI2passt die Weibull-Verteilung (genauer gesagt eine spezielle Parametrisierung davon) am besten zu den Daten. Die genaue Parametrisierung der Verteilung WEI2wird in detaillierten diesem Dokument auf Seite 279. Lassen Sie sich den Sitz kontrollieren , indem sie in einem an den Residuen sucht Wurm Grundstück (im Grunde ein de-tendierte QQ-Plot):

WormPlot

Wir erwarten, dass die Residuen nahe der mittleren horizontalen Linie liegen und zu 95% zwischen der oberen und der unteren gepunkteten Kurve liegen, die als 95% -Punkt-Konfidenzintervalle fungieren. In diesem Fall erscheint mir das Wurmdiagramm gut, was darauf hinweist, dass die Weibull-Verteilung angemessen ist.

COOLSerdash
quelle
1
+1 Schöne Analyse. Eine Frage allerdings. Können positive Schlussfolgerungen zur Kompatibilität mit einer bestimmten Hauptverteilung (in diesem Fall Weibull) das Vorhandensein einer Mischungsverteilung ausschließen? Oder müssen wir eine ordnungsgemäße Gemischanalyse durchführen und GoF prüfen, um diese Option auszuschließen?
Aleksandr Blekh
18
@AleksandrBlekh Es ist unmöglich, genügend Leistung zu haben, um ein Gemisch auszuschließen: Wenn das Gemisch zwei nahezu identische Verteilungen aufweist, kann es nicht nachgewiesen werden, und wenn alle Komponenten bis auf eine sehr kleine Anteile aufweisen, kann es auch nicht nachgewiesen werden. Typischerweise (in Ermangelung einer Theorie, die eine Verteilungsform nahe legt) passt man parametrische Verteilungen an, um sparsame ungefähre Beschreibungen von Daten zu erhalten. Gemische sind keine davon: Sie erfordern zu viele Parameter und sind zu flexibel für den Zweck.
whuber
4
@whuber: +1 Schätzen Sie Ihre ausgezeichnete Erklärung!
Aleksandr Blekh
1
@Lourenco Ich habe mir die Grafik von Cullen und Fey angesehen. Der blaue Punkt kennzeichnet unsere Stichprobe. Sie sehen, dass der Punkt in der Nähe der Linien von Weibull, Lognormal und Gamma liegt (zwischen Weibull und Gamma). Nachdem ich jede dieser Verteilungen angepasst hatte, verglich ich die Anpassungsgütestatistiken mit der Funktion gofstatund dem AIC. Es besteht kein Konsens darüber, wie die "beste" Verteilung am besten ermittelt werden kann. Ich mag grafische Methoden und den AIC.
COOLSerdash
1
@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 anschaue. In diesem Tutorial habe ich beschlossen, es nicht anzuzeigen, um den Beitrag kurz zu halten. Das Lognormal zeigt eine schlechtere Anpassung im Vergleich zur Weibull- und Normalverteilung. Der AIC liegt bei 537,59 und die Grafiken sehen auch nicht gut aus.
COOLSerdash
15

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 Parametern zu plotten, die Sie von fitdistr () erhalten haben.

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

Bildbeschreibung hier eingeben

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

Wenn Sie nur zwei konkurrierende Verteilungen haben (zum Beispiel diejenigen auswählen, die am besten in die Darstellung passen), können Sie mit einem Likelihood-Ratio-Test testen, welche Verteilungen besser passen .

elevendollar
quelle
20
Willkommen bei CrossValidated! Ihre Antwort könnte nützlicher sein, wenn Sie sie so bearbeiten könnten, dass sie (a) den Code enthält, den Sie zum Erstellen der Grafik verwendet haben, und (b) wie man die Grafik lesen würde.
Stephan Kolassa
2
Was wird dort geplottet? Ist das eine Art Exponentialitätsplot?
Glen_b
1
Aber wie entscheiden Sie, welche Distribution am besten zu Ihren Daten passt? Nur anhand der Grafik konnte ich Ihnen nicht sagen, ob logNormal oder weibull am besten zu Ihren Daten passt.
Tobibo
4
Wenn Sie einen Pseudozufallszahlengenerator erstellen möchten, verwenden Sie das empirische cdf. Möchten Sie Zahlen zeichnen, die über Ihre beobachtete Verteilung hinausgehen?
Elevendollar
6
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 sich wirklich um ein ecdf handelt.
Hong Ooi