Jak określić, która dystrybucja najlepiej pasuje do moich danych?

> 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 

Najpierw kilka krótkich komentarzy:

• $ p $ -wartości Kolmovorova -Test-Smirnova (test KS) z oszacowanymi parametrami będzie całkiem błędny. Więc niestety nie możesz po prostu dopasować rozkładu, a następnie użyć oszacowanych parametrów w teście Kołmogorowa-Smirnowa, aby przetestować swoją próbkę.
• Twoja próbka nigdy nie będzie podążać za określonym dystrybucji dokładnie. Więc nawet jeśli Twoje $ p $ -wartości z KS-Test byłyby prawidłowe i $ 0,05 $ , oznaczałoby to po prostu, że nie można wykluczyć , że dane są zgodne z tą konkretną dystrybucją. Innym sformułowaniem byłoby to, że twoja próbka jest zgodna z pewnym rozkładem. Ale odpowiedź na pytanie ” Czy moje dane są dokładnie zgodne z rozkładem xy? ” zawsze brzmi nie.
• Celem nie może być tutaj ustalenie z całą pewnością, jaki rozkład podąża Twoja próbka. Celem jest to, co @whuber (w komentarzach) nazywa oszczędnymi, przybliżonymi opisami danych. Posiadanie określonego rozkładu parametrycznego może być przydatne jako model danych.

---

Ale zróbmy trochę eksploracji. Użyję doskonałego fitdistrplus , który oferuje kilka przydatnych funkcji dopasowywania dystrybucji. Użyjemy funkcji descdist, aby uzyskać kilka pomysłów na temat możliwych dystrybucji kandydatów.

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) 

Teraz użyjmy descdist:

descdist(x, discrete = FALSE) 

Kurtoza i kwadratowa skośność twojej próbki jest wykreślana jako niebieski punkt o nazwie ” Obserwacja „. Wydaje się, że możliwe dystrybucje obejmują rozkład Weibulla, Lognormal i prawdopodobnie rozkład Gamma.

Dopasujmy Rozkład Weibulla i rozkład normalny:

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

Teraz sprawdź dopasowanie do normy:

plot(fit.norm) 

A dla dopasowania Weibull:

plot(fit.weibull) 

Oba wyglądają dobrze, ale sądząc po QQ-Plot, Weibull może wygląda trochę lepiej, zwłaszcza na ogonach. Odpowiednio, AIC dopasowania Weibulla jest niższe w porównaniu do normalnego dopasowania:

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

---

Symulacja testu Kołmogorowa-Smirnowa

Użyję procedury @Aksakal opisanej tutaj , aby zasymulować statystykę KS pod wartością zerową.

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 symulowanej statystyki KS wygląda następująco:

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

Na koniec nasza $ p $ -wartość wykorzystująca symulowaną dystrybucję zerową statystyk KS to:

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

Potwierdza to nasz graficzny wniosek, że próbka jest zgodna z rozkładem Weibulla.

Jak wyjaśniono tutaj , możemy użyć metody ładowania początkowego, aby dodać punktowe przedziały ufności do szacowanego pliku PDF lub CDF Weibulla:

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

---

Automatyczne dopasowanie dystrybucji za pomocą GAMLSS

gamlss pakiet dla R oferuje możliwość wypróbowania wielu różnych dystrybucji i wybrania best ” zgodnie z GAIC (uogólnione kryterium informacyjne Akaike). Główną funkcją jest fitDist. Ważną opcją w tej funkcji jest typ wypróbowanych dystrybucji. Na przykład ustawienie type = "realline" spowoduje wypróbowanie wszystkich zaimplementowanych dystrybucji zdefiniowanych w całej rzeczywistej linii, podczas gdy type = "realsplus" spowoduje wypróbowanie tylko dystrybucji zdefiniowanych na rzeczywistej linii dodatniej . Inną ważną opcją jest parametr $ k $ , który jest karą dla GAIC. W poniższym przykładzie ustawiłem parametr $ k = 2 $ , co oznacza, że ” best jest wybierana zgodnie z klasycznym AIC. Możesz ustawić $ k $ na cokolwiek chcesz, na przykład $ \ log (n) $ dla 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 *** 

Zgodnie z AIC, rozkład Weibulla (a dokładniej WEI2, jego specjalna parametryzacja ) najlepiej pasuje do danych. Dokładna parametryzacja rozkładu WEI2 jest szczegółowo opisana w tym dokumencie na stronie 279. Sprawdźmy dopasowanie, wykonując patrząc na pozostałości na wykresie robaka (w zasadzie wykresie QQ pozbawionym trendu):

Oczekujemy, że reszty będą blisko środkowej poziomej linii, a 95% z nich będzie leżeć między górną i dolne krzywe kropkowane, które stanowią 95% punktowych przedziałów ufności. W tym przypadku wykres robaka wygląda dobrze, co wskazuje, że rozkład Weibulla jest odpowiednio dopasowany.

Komentarze

• +1 Niezła analiza. Jedno pytanie. Czy pozytywny wniosek dotyczący zgodności z konkretnym głównym rozkładem (w tym przypadku Weibull) pozwala wykluczyć możliwość rozkładu mieszaniny ' s? Albo musimy przeprowadzić odpowiednią analizę mieszaniny i sprawdzić GoF, aby wykluczyć tę opcję?
• @AleksandrBlekh Niemożliwe jest posiadanie wystarczającej mocy, aby wykluczyć mieszaninę: gdy mieszanina ma dwa prawie identyczne rozkłady, nie można jej wykryć i gdy wszystkie składniki oprócz jednego mają bardzo małe proporcje nie można go również wykryć. Zazwyczaj (przy braku teorii, która mogłaby sugerować formę dystrybucji), dopasowuje się rozkłady parametryczne w celu uzyskania oszczędnych przybliżonych opisów danych. Mieszaniny nie należą do tych: wymagają zbyt wielu parametrów i są zbyt elastyczne do tego celu.
• @whuber: +1 Doceń swoje doskonałe wyjaśnienie!
• @Lourenco Spojrzałem na wykres Cullena i Feya. Niebieski punkt oznacza naszą próbkę. Widzisz, że punkt znajduje się blisko linii Weibull, Lognormal i Gamma (która znajduje się między Weibull i Gamma). Po dopasowaniu każdego z tych rozkładów porównałem statystyki dobroci dopasowania za pomocą funkcji gofstat i AIC. Nie ma ' konsensusu co do najlepszego sposobu określenia ” najlepszej ” dystrybucji jest. Lubię metody graficzne i AIC.
• @Lourenco Czy masz na myśli lognormal? Rozkład logistyczny (znak ” + „) jest nieco odbiegający od obserwowanych danych. Lognormal byłby również kandydatem, na który ' normalnie patrzę. W tym samouczku ' zdecydowałem się go nie pokazywać, aby post był krótki. Lognormal wykazuje gorsze dopasowanie w porównaniu z rozkładem Weibulla i normalnym. AIC to 537,59, a wykresy również nie ' nie wyglądają zbyt dobrze.

Wykresy to przede wszystkim dobry sposób, aby lepiej zrozumieć, jak wyglądają Twoje dane. W twoim przypadku zalecałbym wykreślenie empirycznej dystrybucyjnej funkcji dystrybucyjnej (ecdf) z teoretycznymi plikami cdf z parametrami uzyskanymi z fitdistr ().

Zrobiłem to raz dla moich danych i uwzględniłem również przedziały ufności. Oto obraz, który otrzymałem za pomocą ggplot2 ().

Czarna linia to empiryczna funkcja dystrybucji kumulatywnej a kolorowe linie to pliki cdf z różnych dystrybucji przy użyciu parametrów, które otrzymałem przy użyciu metody największej wiarygodności. Można łatwo zauważyć, że rozkład wykładniczy i normalny nie jest dobrze dopasowany do danych, ponieważ linie mają inną postać niż ecdf, a linie są dość daleko od ecdf. Niestety inne dystrybucje są dość zbliżone. Ale powiedziałbym, że linia logNormal jest najbliższa czarnej linii. Używając miary odległości (na przykład MSE), można by potwierdzić założenie.

Jeśli masz tylko dwa konkurujące rozkłady (na przykład wybierając te, które wydają się najlepiej pasować do wykresu), możesz użyć Test współczynnika prawdopodobieństwa , aby sprawdzić, które dystrybucje pasują lepiej.

Komentarze

• Witamy w CrossValidated! Twoja odpowiedź mogłaby być bardziej przydatna, gdybyś mógł ją edytować, aby zawierała (a) kod użyty do stworzenia grafiki i (b) jak można ją odczytać.
• Co tam jest kreślone? Czy to jakiś wykres wykładniczy?
• Ale jak zdecydować, który rozkład najlepiej pasuje do danych? Tylko zgodnie z grafiką nie mogłem ' powiedzieć, czy logNormal czy weibull najlepiej pasuje do twoich danych.
• Jeśli chcesz utworzyć generator liczb pseudolosowych, dlaczego nie używać empirycznego cdf? Czy chcesz narysować liczby wykraczające poza obserwowany rozkład?
• Przyjmując wykres za wartość nominalną, wydaje się, że żaden z Twoich rozkładów kandydatów w ogóle nie pasuje do danych. Wygląda na to, że Twój plik ecdf ma asymptotę poziomą na poziomie mniejszym niż 0,03, co nie ' nie ma sensu, więc ' nie jestem pewien, czy to naprawdę jest ecdf na pierwszym miejscu.