Jak ustalić, która dystrybucja najlepiej pasuje do moich danych?

133

Mam zestaw danych i chciałbym dowiedzieć się, która dystrybucja najlepiej pasuje do moich danych.

Użyłem tej fitdistr()funkcji do oszacowania niezbędnych parametrów do opisania założonego rozkładu (tj. Weibull, Cauchy, Normal). Korzystając z tych parametrów, mogę przeprowadzić test Kołmogorowa-Smirnowa, aby oszacować, czy moje przykładowe dane pochodzą z tego samego rozkładu, co założony rozkład.

Jeśli wartość p wynosi> 0,05, mogę założyć, że przykładowe dane pochodzą z tego samego rozkładu. Ale wartość p nie dostarcza żadnych informacji o boskości dopasowania, prawda?

Jeśli więc wartość p moich przykładowych danych wynosi> 0,05 dla rozkładu normalnego, a także rozkładu Weibulla, to skąd mam wiedzieć, który rozkład lepiej pasuje do moich danych?

Zasadniczo to zrobiłem:

> 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

Wartości p wynoszą 0,8669 dla rozkładu Weibulla i 0,5522 dla rozkładu normalnego. Dlatego mogę założyć, że moje dane są zgodne z Weibullem, a także z normalnym rozkładem. Ale która funkcja dystrybucji lepiej opisuje moje dane?


Odnosząc się do elevendollar znalazłem następujący kod, ale nie wiem, jak interpretować wyniki:

fits <- list(no = fitdistr(mydata, "normal"),
             we = fitdistr(mydata, "weibull"))
sapply(fits, function(i) i$loglik)
       no        we 
-259.6540 -257.9268 
tobibo
źródło
5
Dlaczego chcesz dowiedzieć się, która dystrybucja najlepiej pasuje do twoich danych?
Roland
6
Ponieważ chcę wygenerować liczby pseudolosowe po podanym rozkładzie.
tobibo
6
Nie można użyć KS do sprawdzenia, czy rozkład z parametrami znalezionymi w zestawie danych jest zgodny z zestawem danych. Zobacz na przykład # 2 na tej stronie , a także alternatywy (i inne sposoby, w których test KS może wprowadzać w błąd).
tpg2114
Kolejna dyskusja tutaj z przykładowymi kodami, jak zastosować test KS, gdy parametry są szacowane z próbki.
Aksakal
1
I used the fitdistr() function ..... Jaka jest fitdistrfunkcja? Coś z Excela? A może coś napisałeś w C?
wilki

Odpowiedzi:

162

Po pierwsze, oto kilka szybkich komentarzy:

  • Wartości testu Kołmoworowa-Smirnowa (test KS) z oszacowanymi parametrami będą zupełnie błędne. Niestety nie można po prostu dopasować rozkładu, a następnie użyć oszacowanych parametrów w teście Kołmogorowa-Smirnowa, aby przetestować próbkę.p
  • Twoja próbka nigdy nie będzie dokładnie odpowiadać określonej dystrybucji. Więc nawet jeśli twoje wartości z testu KS byłyby prawidłowe i , oznaczałoby to po prostu, że nie możesz wykluczyć, że twoje dane są zgodne z tym konkretnym rozkładem. Innym sformułowaniem byłoby, że twoja próbka jest kompatybilna z pewnym rozkładem. Ale odpowiedź na pytanie „Czy moje dane dokładnie odpowiadają rozkładowi xy?” zawsze jest nie.p>0.05
  • Celem tutaj nie może być ustalenie, z jakim rozkładem następuje próbka. Celem jest to, co @whuber (w komentarzach) nazywa oszczędnym przybliżonym opisem danych. Posiadanie określonego rozkładu parametrycznego może być przydatny jako model danych.

Ale zajmijmy się eksploracją. Wykorzystam doskonały fitdistrpluspakiet, który oferuje kilka fajnych funkcji do dopasowania dystrybucji. Użyjemy tej funkcji, descdistaby uzyskać pomysły 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 pozwala użyć descdist:

descdist(x, discrete = FALSE)

Descdist

Kurtoza i kwadratowa skośność próbki jest wykreślona jako niebieski punkt o nazwie „Obserwacja”. Wydaje się, że możliwe rozkłady obejmują Weibulla, Lognormala i być może 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 normalnego:

plot(fit.norm)

Normalne dopasowanie

A dla dopasowania Weibull:

plot(fit.weibull)

Dopasowanie Weibull

Oba wyglądają dobrze, ale sądząc po QQ-Plot, Weibull może wygląda nieco lepiej, szczególnie na ogonach. Odpowiednio, AIC dopasowania Weibull 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

Wykorzystam wyjaśnioną tutaj procedurę @ Aksakala, 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()

Symulowane statystyki KS

Wreszcie, nasza wartość wykorzystująca symulowany rozkład zerowy statystyki KS wynosi: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

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

Jak wyjaśniono tutaj , możemy użyć ładowania początkowego, aby dodać punktowe przedziały ufności do szacowanego pliku Weibull PDF lub 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)

CI_Density

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


Automatyczne dopasowanie rozdzielacza z GAMLSS

gamlssPakiet Rdaje możliwość spróbować wielu różnych rozkładów i wybierz „najlepsze” według GAIC (uogólniona Akaike kryterium informacji). Główną funkcją jest fitDist. Ważną opcją w tej funkcji jest typ wypróbowanych dystrybucji. Na przykład ustawienie type = "realline"spróbuje wypróbować wszystkie zaimplementowane rozkłady zdefiniowane na całej linii rzeczywistej, podczas gdy type = "realsplus"wypróbuje tylko rozkłady zdefiniowane na prawdziwej linii dodatniej. Inną ważną opcją jest parametr , który jest karą dla GAIC. W poniższym przykładzie ustawiłem parametr co oznacza, że ​​„najlepszy” rozkład jest wybierany zgodnie z klasycznym AIC. Możesz ustawić na dowolne, na przykładkk=2klog(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 ***

Według AIC rozkład Weibulla (a konkretniej WEI2jego specjalna parametryzacja) najlepiej pasuje do danych. Dokładna parametryzacja rozkładu WEI2jest szczegółowo opisana w tym dokumencie na stronie 279. Sprawdźmy dopasowanie, patrząc na resztki na wykresie robaka (w zasadzie zniekształconego wykresu QQ):

WormPlot

Oczekujemy, że reszty będą zbliżone do środkowej linii poziomej, a 95% z nich będzie leżeć między górną i dolną krzywą kropkowaną, które działają jak 95% przedziały ufności punktowej. W tym przypadku wykres robaka wygląda dla mnie dobrze, wskazując, że rozkład Weibulla jest odpowiedni.

COOLSerdash
źródło
1
+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ść obecności rozkładu mieszaniny? Czy musimy przeprowadzić odpowiednią analizę mieszaniny i sprawdzić GoF, aby wykluczyć tę opcję?
Aleksandr Blekh
18
@AleksandrBlekh Nie można mieć wystarczającej mocy, aby wykluczyć mieszaninę: gdy mieszanina ma dwa prawie identyczne rozkłady, nie można jej wykryć, a gdy wszystkie składniki oprócz jednego mają bardzo małe proporcje, nie można jej również wykryć. Zazwyczaj (przy braku teorii, która mogłaby sugerować formę dystrybucyjną), dopasowuje się rozkłady parametryczne, aby uzyskać skąpe przybliżone opisy danych. Mieszaniny nie są żadnymi z nich: wymagają zbyt wielu parametrów i są zbyt elastyczne do tego celu.
whuber
4
@whuber: +1 Doceń swoje doskonałe wyjaśnienie!
Aleksandr Blekh
1
@Lourenco Spojrzałem na wykres Cullen i Fey. Niebieski punkt oznacza naszą próbkę. Widzisz, że punkt jest zbliżony do linii Weibulla, Lognormala i Gammy (która jest pomiędzy Weibullem a Gammą). Po dopasowaniu każdego z tych rozkładów porównałem statystyki dobroci dopasowania za pomocą funkcji gofstati AIC. Nie ma zgody co do tego, jaki jest najlepszy sposób określenia „najlepszego” rozkładu. Lubię metody graficzne i AIC.
COOLSerdash
1
@Lourenco Czy masz na myśli lognormal? Rozkład logistyczny (znak „+”) jest dość oddalony od obserwowanych danych. Lognormal byłby również kandydatem, na którego normalnie patrzę. W tym samouczku postanowiłem nie pokazywać go, aby post nie 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 wyglądają zbyt dobrze.
COOLSerdash
15

Wykresy to przede wszystkim dobry sposób, aby uzyskać lepszy obraz tego, jak wyglądają Twoje dane. W twoim przypadku poleciłbym wykreślenie empirycznej funkcji kumulatywnej dystrybucji (ecdf) względem teoretycznych plików cdf z parametrami uzyskanymi z fitdistr ().

Zrobiłem to raz dla moich danych, a także uwzględniłem przedziały ufności. Oto zdjęcie, które otrzymałem za pomocą ggplot2 ().

wprowadź opis zdjęcia tutaj

Czarna linia jest empiryczną funkcją kumulatywnego rozkładu, a kolorowe linie to pliki cdf z różnych rozkładów przy użyciu parametrów, które otrzymałem przy użyciu metody Maximum Likelihood. Łatwo zauważyć, że rozkład wykładniczy i normalny nie jest dobrze dopasowany do danych, ponieważ linie mają inną formę niż ecdf, a linie są dość daleko od ecdf. Niestety inne dystrybucje są dość bliskie. Powiedziałbym jednak, że linia logNormal jest najbliższa czarnej linii. Za pomocą miary odległości (na przykład MSE) można zweryfikować to założenie.

Jeśli masz tylko dwie konkurencyjne dystrybucje (na przykład wybierając te, które wydają się najlepiej pasować do wykresu), możesz użyć testu ilorazu wiarygodności, aby sprawdzić, które dystrybucje pasują lepiej.

jedenastoletni
źródło
20
Witamy w CrossValidated! Twoja odpowiedź może być bardziej przydatna, jeśli możesz ją edytować w celu włączenia (a) kodu użytego do wytworzenia grafiki oraz (b) sposobu, w jaki odczytujesz grafikę.
Stephan Kolassa
2
Co się tam dzieje? Czy to jakiś wykres wykładniczy?
Glen_b
1
Ale w jaki sposób decydujesz, która dystrybucja najlepiej pasuje do twoich danych? Tylko zgodnie z grafiką nie mogłem powiedzieć, czy logNormal czy weibull najlepiej pasuje do twoich danych.
tobibo
4
Jeśli chcesz stworzyć generator liczb pseudolosowych, skorzystaj z empirycznego cdf? Czy chcesz narysować liczby, które wykraczają poza obserwowany rozkład?
elevendollar
6
Biorąc wykres za wartość nominalną, wydaje się, że żadna z twoich kandydujących dystrybucji w ogóle nie pasuje do danych. Ponadto wydaje się, że twój ecdf ma poziomą asymptotę mniejszą niż 0,03, co nie ma sensu, więc nie jestem pewien, czy to naprawdę jest ecdf.
Hong Ooi