Van egy adatkészletem, és szeretném kideríteni, melyik eloszlás felel meg legjobban az adataimnak.
A fitdistr() függvénnyel becsültem meg a feltételezett eloszlás leírásához szükséges paramétereket (azaz Weibull, Cauchy, Normal). Ezen paraméterek felhasználásával Kolmogorov-Smirnov tesztet tudok végezni annak becslésére, hogy a mintaadataim ugyanarról az eloszlásról származnak-e, mint a feltételezett eloszlásom. ugyanabból az eloszlásból merített. De a p-érték nem ad információt az illeszkedés isteniségéről, nem igaz?
Tehát, ha a mintaadataim p-értéke> 0,05 egy normál eloszláshoz, valamint egy weibull eloszláshoz, honnan tudhatom, melyik eloszlás illik jobban az adataimhoz?
Alapvetően ezt tettem:
> 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 A p-értékek 0,8669 a Weibull eloszlásnál és 0,5522 a normális eloszlás. Így feltételezhetem, hogy adataim Weibull mellett normális eloszlást is követnek. De melyik terjesztési függvény írja le jobban az adataimat?
Az elevendollar ra hivatkozva a következő kódot találtam, de nem tudom, hogyan kell értelmezni az eredményeket:
fits <- list(no = fitdistr(mydata, "normal"), we = fitdistr(mydata, "weibull")) sapply(fits, function(i) i$loglik) no we -259.6540 -257.9268 Megjegyzések
Válasz
Először íme néhány gyors megjegyzés:
- A Kolmovorov $ p $ -értékei -Smirnov-Test (KS-Test) a becsült paraméterekkel elég hibás lesz. Tehát sajnos nem lehet csak elosztást illeszteni, majd a Kolmogorov-Smirnov-teszt becsült paramétereit felhasználni a minta teszteléséhez.
- A mintája soha nem követ egy adott pontosan akkor is, ha a KS-tesztből származó $ p $ -értékei érvényesek lennének, és $ > 0,05 $ , ez csak azt jelentené, hogy nem zárhatja ki , hogy az adatai ezt a konkrét terjesztést kövessék. Egy másik megfogalmazás az lenne, hogy a mintád kompatibilis egy bizonyos eloszlással. De a válasz a ” kérdésre pontosan követi-e az xy elosztást? ” mindig nem .
- A cél itt nem lehet az, hogy pontosan meghatározzuk, milyen eloszlást követ a mintája. A cél az, amit @whuber (a megjegyzésekben) az adatok párhuzamos hozzávetőleges leírása nak nevez. Egy adott paraméteres eloszlás hasznos lehet az adatok modelljeként.
De végezzünk egy kis feltárást. A kiváló fitdistrplus csomag, amely néhány szép funkciót kínál a disztribúció illesztéséhez. A descdist függvényt nyerjük néhány ötlet a lehetséges jelöltek terjesztéséről.
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) Most használhatjuk a descdist alkalmazást:
descdist(x, discrete = FALSE) 
A minta kurtosisa és négyzetbeli ferdesége kék pontként ábrázolva, a div neve: div div = “2faea39abb”>
Megfigyelés “. Úgy tűnik, hogy a lehetséges disztribúciók magukban foglalják a Weibull, a Lognormal és esetleg a Gamma eloszlást.
Legyen “illő” Weibull-eloszlás és normál eloszlás:
fit.weibull <- fitdist(x, "weibull") fit.norm <- fitdist(x, "norm") Most ellenőrizze a normálist:
plot(fit.norm) 
És a Weibull illeszkedéshez:
plot(fit.weibull) 
Mindkettő jól néz ki, de a QQ-Plot alapján megítélve a Weibull talán egy kicsit jobban néz ki, főleg a faroknál. Ennek megfelelően a Weibull-illesztés AIC-értéke alacsonyabb a normál illeszkedéshez képest:
fit.weibull$aic [1] 519.8537 fit.norm$aic [1] 523.3079 Kolmogorov-Smirnov tesztszimuláció
Az @Aksakal itt leírt eljárást a KS-statisztika szimulálására használom a null alatt.
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 ) }) A szimulált KS-statisztika ECDF-je a következőképpen néz ki:
plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7) grid() 
Végül a $ p $ -értékünk a szimulált null-elosztás felhasználásával A KS-statisztika a következő:
fit <- logspline(stats) 1 - plogspline(ks.test(x , "pweibull" , shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"])$statistic , fit ) [1] 0.4889511 Ez megerősíti azt a grafikus következtetést, hogy a minta kompatibilis egy Weibull-terjesztéssel.
Amint azt elmagyaráztuk itt , a bootstrapeléssel pontszerű megbízhatósági intervallumokat adhatunk hozzá a becsült Weibull PDF vagy CDF fájlhoz:
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") 
Automatikus terjesztés illesztése a GAMLSS-szel
A gamlss a R csomagja lehetőséget kínál arra, hogy sokféle disztribúciót kipróbálhasson, és kiválassza a best ” a GAIC szerint (az általánosított Akaike információs kritérium). A fő funkció a fitDist. Fontos lehetőség ebben a függvényben a kipróbált disztribúciók típusa. Például az type = "realline" beállítás az összes megvalósított disztribúciót kipróbálja az egész valós vonalon, míg a type = "realsplus" csak a valódi pozitív vonalon definiált disztribúciókat próbálja meg. . Egy másik fontos lehetőség a $ k $ paraméter, amely a GAIC büntetését jelenti. Az alábbi példában beállítottam a $ k = 2 $ paramétert, ami azt jelenti, hogy a ” legjobb terjesztést a klasszikus AIC szerint választják meg. A $ k $ értéket bármire beállíthatja, például $ \ log (n) $ 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 *** Az AIC szerint a Weibull eloszlás (pontosabban WEI2, ennek egy speciális paraméterezése ) az adatoknak felel meg a legjobban. Az WEI2 eloszlás pontos paraméterezését ebben a dokumentumban a 279. oldalon részletezzük. Ellenőrizzük az illesztést a maradékokat megnézni egy féregtáblázatban (alapvetően egy trend nélküli QQ-diagramban):
Arra számítunk, hogy a maradványok közel vannak a középső vízszintes vonalhoz, és 95% -uk a felső felső közé esik és az alacsonyabb pontozott görbék, amelyek 95% -os pontosságú megbízhatósági intervallumként működnek. Ebben az esetben a féregtáblázat jól néz ki számomra, jelezve, hogy a Weibull-eloszlás megfelelő illeszkedés.
Megjegyzések
- +1 Szép elemzés. Azonban egy kérdés. A pozitív következtetés a kompatibilitásról egy adott fő eloszlással (ebben az esetben Weibull) lehetővé teszi-e a keverékeloszlás lehetőségének kizárását ‘ jelenléte? Vagy megfelelő keverékelemzést kell végeznünk, és ellenőriznünk kell a GoF kizárja ezt a lehetőséget?
- @AleksandrBlekh Lehetetlen elegendő erővel rendelkezni egy keverék kizárásához: ha a keverék két, szinte azonos eloszlású, akkor nem lehet kimutatni, és ha az összes kivételével az összes komponensnek nagyon kicsi az aránya azt sem lehet kimutatni. Jellemzően (disztribúciós formára utaló elmélet hiányában) a paraméteres eloszlások illeszkednek az adatok párhuzamos hozzávetőleges leírása eléréséhez. A keverékek nem tartoznak ezek közé: túl sok paramétert igényelnek, és túl rugalmasak a cél érdekében.
- @whuber: +1 Értékelje kitűnő magyarázatát!
- @Lourenco Megnéztem a Cullen és Fey grafikont. A kék pont a mintánkat jelöli. Látja, hogy a pont közel van a Weibull, a Lognormal és a Gamma vonalához (amely a Weibull és a Gamma között van). Miután mindegyik disztribúciót illesztettem, összehasonlítottam az illeszkedés statisztikáját a
gofstatfüggvény és az AIC segítségével. Nincs ‘ ta konszenzus arról, hogy mi a legjobb módszer a ” legjobb ” eloszlás meghatározására van. Szeretem a grafikus módszereket és az AIC-t. - @Lourenco A lognormálisra gondolsz? A logisztikai eloszlás (a ” + ” jel) meglehetősen távol áll a megfigyelt adatoktól. A lognormális egy olyan jelölt is, akit általában nézek. Ehhez az oktatóanyaghoz ‘ úgy döntöttem, hogy nem jelenik meg annak érdekében, hogy rövid legyen a bejegyzés. A lognormal rosszabb illeszkedést mutat mind a Weibull, mind a Normal eloszláshoz képest. Az AIC értéke 537,59, és a grafikonok szintén nem néznek ki túlságosan jónak ‘.
Válasz
A diagramok többnyire jó módszerek arra, hogy jobban megértsék adatait. Az Ön esetében javasolnám a empirikus kumulatív elosztási függvény (ecdf) ábrázolását az fitdistr-től kapott paraméterekkel ellátott elméleti cdf-ekkel szemben ().
Ezt egyszer megcsináltam az adataimhoz, és a bizalmi intervallumokat is belefoglaltam. Itt van a kép, amelyet a ggplot2 () használatával kaptam.
A fekete vonal az empirikus kumulatív elosztási függvény és a színes vonalak különböző eloszlásokból származó cdf-k, a Maximum Likelihood módszerrel kapott paraméterek felhasználásával. Könnyen belátható, hogy az exponenciális és a normál eloszlás nem illik jól az adatokhoz, mert a vonalaknak más formájuk van, mint az ecdf-nek, és a vonalak meglehetősen messze vannak az ecdf-től. Sajnos a többi disztribúció meglehetősen közel áll egymáshoz. De azt mondanám, hogy a logNormal vonal áll a legközelebb a fekete vonalhoz. A távolság mérésével (például MSE) érvényesíteni lehet a feltételezést.
Ha csak két versengő eloszlásod van (például kiválasztod azokat, amelyek a legjobban illeszkednek a diagramhoz), használhatsz egy Likelihood-Ratio-Test annak tesztelésére, hogy melyik eloszlás illik jobban.
Megjegyzések
- Üdvözöljük a CrossValidated oldalán! Válasza hasznosabb lehet, ha szerkesztheti (a) a grafika előállításához használt kódot és (b) hogyan olvassa el a grafikát.
- Mi van ott ábrázolva? Ez valamiféle exponenciális ábrázolás?
- De hogyan döntheti el, hogy melyik terjesztés felel meg az adatainak a legjobban? Csak a grafika szerint nem tudtam ‘ megmondani, hogy a logNormal vagy a weibull illik-e a legjobban az adataidhoz.
- Ha ál-véletlenszerű számgenerátort akarsz létrehozni nem használja az empirikus cdf-t? Szeretne olyan számokat rajzolni, amelyek meghaladják a megfigyelt eloszlást?
- Ha grafikonját névértékre vesszük, akkor úgy tűnik, hogy a jelölt disztribúciók közül egyik egyáltalán nem felel meg jól az adatoknak. Ezenkívül úgy tűnik, hogy az ecdf vízszintes aszimptotája kevesebb, mint 0,03, aminek ‘ nincs értelme, ezért nem vagyok biztos abban, hogy ‘ ez valójában egy ecdf.
I used the fitdistr() function… ..Milyen ‘ sfitdistrfunkció? Valami az Excelből? Vagy valami, amit magad írtál C-be?