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

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 

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)

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)

A dla dopasowania Weibull:

plot(fit.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()

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)

#-----------------------------------------------------------------------------
# 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 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ład$k$$k = 2$$k$$\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 ***

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

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.

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 ().

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.