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