Jak zjistit, která distribuce nejlépe vyhovuje mým datům?

> 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 

Nejprve uvádíme několik rychlých komentářů:

• Hodnoty $ p $ – hodnoty Kolmovorov -Smirnovův test (KS-test) s odhadovanými parametry bude docela špatný. Bohužel tedy nemůžete přizpůsobit distribuci a poté použít odhadované parametry v Kolmogorovově-Smirnovově testu k otestování svého vzorku.
• Váš vzorek nikdy nebude následovat konkrétní distribuce přesně. Takže i kdyby vaše $ p $ hodnoty z KS testu byly platné a $ > 0,05 $ , znamenalo by to jen to, že nemůžete vyloučit , že vaše data sledují tuto konkrétní distribuci. Další formulace by byla, že váš vzorek je kompatibilní s určitou distribucí. Ale odpověď na otázku “ Sledují moje data přesně distribuci xy? “ je vždy ne.
• Cílem zde nemůže být s jistotou určit, jaké rozdělení následuje váš vzorek. Cílem je to, co @ whuber (v komentářích) nazývá šetrné přibližné popisy dat. Mít specifickou parametrickou distribuci může být užitečné jako model dat.

---

Ale pojďme udělat nějaký průzkum. Použiji vynikající fitdistrplus balíček, který nabízí některé pěkné funkce pro distribuční přizpůsobení. Funkci descdist použijeme k získání několik představ o možných distribucích kandidátů.

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) 

Nyní si dovolíme použít descdist:

descdist(x, discrete = FALSE) 

Křivost a čtvercová šikmost vašeho vzorku se vykreslí jako modrý bod s názvem “ Pozorování „. Zdá se, že k možným distribucím patří Weibullova, Lognormální a případně i gama distribuce.

Pojďme Weibullova distribuce a normální distribuce:

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

Nyní zkontrolujte vhodnost pro normální:

plot(fit.norm) 

A pro Weibull fit:

plot(fit.weibull) 

Oba vypadají dobře, ale soudě podle QQ-Plot, Weibull možná vypadá o něco lépe, zejména na ocasy. Odpovídajícím způsobem je AIC Weibullova přizpůsobení nižší ve srovnání s normálním přizpůsobením:

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

---

Kolmogorov-Smirnovova simulace testu

Použiji @Aksakalův postup vysvětlený zde , abych simuloval statistiku KS pod nulou.

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 simulovaných statistik KS vypadá takto:

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

Nakonec naše $ p $ hodnota pomocí simulované nulové distribuce statistik KS je:

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

To potvrzuje náš grafický závěr, že vzorek je kompatibilní s Weibullovou distribucí.

Jak je vysvětleno zde , můžeme pomocí bootstrappingu přidat bodové intervaly spolehlivosti do odhadovaného Weibullova PDF nebo 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") 

---

Automatická distribuční armatura s GAMLSS

gamlss balíček pro R nabízí možnost vyzkoušet mnoho různých distribucí a vybrat nejlepší “ podle GAIC (zobecněné informační kritérium Akaike). Hlavní funkcí je fitDist. Důležitou možností v této funkci je typ distribucí, které jsou vyzkoušeny. Například nastavení type = "realline" vyzkouší všechny implementované distribuce definované na celé reálné linii, zatímco type = "realsplus" zkusí pouze distribuce definované na skutečné kladné linii . Další důležitou možností je parametr $ k $ , což je pokuta za GAIC. V níže uvedeném příkladu nastavím parametr $ k = 2 $ , což znamená, že “ nejlepší se vybírá podle klasického AIC. $ k $ můžete nastavit na cokoli, co se vám líbí, například $ \ log (n) $ pro 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 *** 

Podle AIC je Weibullova distribuce (konkrétněji WEI2, její speciální parametrizace) ) nejlépe vyhovuje datům. Přesná parametrizace distribuce WEI2 je podrobně popsána v tomto dokumentu na straně 279. Nechte jej zkontrolovat pomocí při pohledu na rezidua v červovém spiknutí (v podstatě de-trendovaný QQ-spiknutí):

Očekáváme, že zbytky budou blízko střední vodorovné čáry a 95% z nich bude ležet mezi horní a spodní tečkované křivky, které fungují jako 95% bodové intervaly spolehlivosti. V tomto případě mi červ vypadá dobře, což naznačuje, že Weibullova distribuce je adekvátní fit.

Komentáře

• +1 Pěkná analýza. Jedna otázka však. Má pozitivní závěr o kompatibilitě s konkrétní hlavní distribucí (v tomto případě Weibull), umožňuje vyloučit možnost distribuce směsi ‚ s přítomností? Nebo musíme provést správnou analýzu směsi a zkontrolovat GoF na vyloučit tuto možnost?
• @AleksandrBlekh Je nemožné mít dostatek energie k vyloučení směsi: pokud je směs dvou téměř identických distribucí, nelze ji detekovat a když mají všechny složky kromě jedné velmi malé rozměry ani to nelze zjistit. Typicky (při absenci teorie, která by mohla naznačovat distribuční formu) se hodí parametrické distribuce, aby se dosáhlo šetrných přibližných popisů dat. Směsi nejsou ničím z nich: vyžadují příliš mnoho parametrů a jsou příliš pro tento účel flexibilní.
• @whuber: +1 Oceňujte své vynikající vysvětlení!
• @Lourenco Podíval jsem se na graf Cullen a Fey. Modrý bod označuje náš vzorek. Vidíte, že bod je blízko čar Weibullova, Lognormálního a Gama (což je mezi Weibullovým a Gama). Po přizpůsobení každé z těchto distribucí jsem porovnal statistiku dobré shody pomocí funkce gofstat a AIC. Není ‚ shoda ohledně toho, jaký je nejlepší způsob určení “ nejlepší “ distribuce je. Mám rád grafické metody a AIC.
• @Lourenco Myslíš lognormální? Logistická distribuce (znak “ + „) je dost daleko od pozorovaných údajů. Lognormal by byl také kandidátem, na který se normálně podívám ‚. V tomto výukovém programu jsem se ‚ rozhodl jej nezobrazovat, aby byl příspěvek krátký. Lognormal vykazuje horší přizpůsobení ve srovnání s Weibullovým i normálním rozdělením. AIC je 537,59 a grafy také nevypadají příliš dobře.

Grafy jsou většinou dobrým způsobem, jak získat lepší představu o tom, jak vaše data vypadají. Ve vašem případě bych doporučil vykreslit empirickou kumulativní distribuční funkci (ecdf) proti teoretickým cdfs s parametry, které jste dostali od fitdistr ().

Udělal jsem to jednou pro svá data a zahrnoval také intervaly spolehlivosti. Tady je obrázek, který mám pomocí ggplot2 ().

Černá čára je empirická kumulativní distribuční funkce a barevné čáry jsou CDF z různých distribucí pomocí parametrů, které jsem dostal pomocí metody Maximum Likelihood. Lze snadno vidět, že exponenciální a normální rozdělení se nehodí k datům, protože řádky mají jinou formu než ecdf a řádky jsou od ecdf dost daleko. Ostatní distribuce jsou bohužel docela blízko. Ale řekl bych, že logNormal linka je nejblíže černé linii. Použitím míry vzdálenosti (například MSE) by bylo možné ověřit předpoklad.

Pokud máte pouze dvě konkurenční distribuce (například výběr těch, které se nejvíc hodí do grafu), můžete použít Likelihood-Ratio-Test k testování, které distribuce se hodí lépe.

Komentáře

• Vítejte v CrossValidated! Vaše odpověď by mohla být užitečnější, kdybyste ji mohli upravit tak, aby obsahovala (a) kód, který jste použili k vytvoření grafiky, a (b) to, jak by se dalo číst grafiku.
• Co se tam vykresluje? Je to nějaký druh exponenciálního grafu?
• Ale jak se rozhodnete, která distribuce nejlépe vyhovuje vašim datům? Pouze podle obrázku vám nemohu ‚ říci, zda se k vašim datům nejlépe hodí logNormal nebo weibull.
• Pokud chcete vytvořit generátor pseudonáhodných čísel, proč nepoužívat empirický CDF? Chcete nakreslit čísla, která přesahují vaši pozorovanou distribuci?
• Vezmeme-li graf v nominální hodnotě, zdá se, že žádná z vašich kandidátských distribucí se k datům vůbec nehodí. Zdá se také, že váš soubor ecdf má horizontální asymptotu na méně než 0,03, což však ‚ nedává smysl, takže si ‚ nejsem jistý, že na prvním místě je to skutečně ekdf.