Kuinka selvittää, mikä jakelu sopii parhaiten tietooni?

> 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 

Ensinnäkin tässä on joitain pikakommentteja:

• Kolmovorovin $ p $ -arvot -Smirnov-testi (KS-testi), jossa on arvioitu parametri, on melko väärä. Joten valitettavasti et voi vain sovittaa jakaumaa ja käyttää sitten arvioituja parametreja Kolmogorov-Smirnov-testissä testataksesi näytteesi.
• Näytteesi ei koskaan seuraa tiettyä jakauma täsmälleen. Joten vaikka $ p $ -arvosi KS-testistä olisivat kelvollisia ja $ > 0,05 $ , se tarkoittaisi vain sitä, että et voi sulkea pois sitä, että tietosi seuraavat tätä tiettyä jakelua. Toinen muotoilu olisi, että näyte on yhteensopiva tietyn jakauman kanssa. Mutta vastaus kysymykseen ” Seuraavatko tietoni tarkkaan jakaumaa xy? ” on aina ei .
• Tavoite ei voi tässä olla määrittää varmasti, mitä jakaumaa näytteesi seuraa. Tavoitteena on, mitä @whuber (kommenteissa) kutsuu tietojen yksinkertaisiksi likimääräisiksi kuvauksiksi . Tietyn parametrisen jakauman käyttäminen voi olla hyödyllistä tietojen mallina.

---

Mutta tehkäämme etsintää. Käytän erinomaista fitdistrplus paketti, joka tarjoaa hienoja toimintoja jakeluasennukseen. Käytämme funktiota descdist saadaksesi joitain ideoita mahdollisista ehdokasjakoista.

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) 

Nyt voidaan käyttää descdist:

descdist(x, discrete = FALSE) 

Näytteen kurtoosi ja neliön vinous on piirretty sinisenä pisteenä nimeltä ” Tarkkailu ”. Vaikuttaa siltä, että mahdolliset jakelut sisältävät Weibull-, Lognormal- ja mahdollisesti Gamma-jakelun.

Anna ”sovittaa” Weibull-jakauma ja normaalijakauma:

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

Tarkista nyt normaalin sopivuus:

plot(fit.norm) 

Ja Weibull-sovitus:

plot(fit.weibull) 

Molemmat näyttävät hyvältä, mutta QQ-Plotin perusteella Weibull näyttää ehkä hieman paremmalta, etenkin hännän kohdalta. Vastaavasti Weibull-sovituksen AIC on pienempi kuin normaali sovitus:

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

---

Kolmogorov-Smirnov-testisimulaatio

Käytän @Aksakalin menettelyä, joka on selitetty täällä simuloimaan KS-tilastoa nollan alla.

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

Simuloidun KS-tilaston ECDF näyttää tältä:

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

Lopuksi $ p $ -arvo käyttämällä simuloitua nollajakaumaa KS-tilastoista on:

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

Tämä vahvistaa graafisen johtopäätöksemme siitä, että näyte on yhteensopiva Weibull-jakelun kanssa.

Kuten selitettiin täällä , voimme käynnistyshihnan avulla lisätä pistemäiset luottamusvälit arvioituun Weibull-PDF-tiedostoon tai CDF: ään:

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

---

Automaattinen jakelu GAMLSS: n kanssa

gamlss paketti kohteelle R tarjoaa mahdollisuuden kokeilla monia erilaisia jakeluja ja valita paras ” GAIC: n (yleinen Akaike-tietokriteeri) mukaan. Päätoiminto on fitDist. Tärkeä vaihtoehto tässä toiminnossa on kokeiltujen jakelujen tyyppi. Esimerkiksi asettamalla type = "realline" yritetään kaikkia toteutettuja jakeluja, jotka on määritelty koko reaalirivillä, kun taas type = "realsplus", vain todellisella positiivisella linjalla määriteltyjä jakeluja. . Toinen tärkeä vaihtoehto on parametri $ k $ , joka on rangaistus GAIC: lle. Seuraavassa esimerkissä asetin parametrin $ k = 2 $ , mikä tarkoittaa, että ” paras ” jakelu valitaan klassisen AIC: n mukaan. Voit asettaa $ k $ mihin tahansa haluamaasi, kuten $ \ 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 *** 

AIC: n mukaan Weibull-jakauma (tarkemmin sanottuna WEI2, sen erityinen parametrisointi ) sopii parhaiten tietoihin. Jakelun WEI2 tarkka parametrointi on selostettu tässä asiakirjassa sivulla 279. Tarkastakaamme sovitus tarkastelemalla jäännöksiä matomatkalla (pohjimmiltaan suuntaamaton QQ-juoni):

Odotamme, että jäännökset ovat lähellä keskimmäistä vaakasuoraa viivaa ja 95% niistä on ylemmän välissä ja alemmat pisteviivat käyrät, jotka toimivat 95%: n pisteiden luottamusväleinä. Tässä tapauksessa mato-juoni näyttää hyvältä osoittaen, että Weibull-jakauma on sopiva.

Kommentit

• +1 Hieno analyysi. Yksi kysymys kuitenkin. Onko positiivinen johtopäätös yhteensopivuudesta tietyn pääjakauman kanssa (tässä tapauksessa Weibull), voidaan sulkea pois seosjakauman mahdollisuus ’ läsnäolo? Tai meidän on suoritettava asianmukainen seosanalyysi ja tarkistettava sulkea pois tämä vaihtoehto?
• @AleksandrBlekh On mahdotonta saada tarpeeksi voimaa sulkea pois seos: kun seoksella on kaksi melkein identtistä jakaumaa, sitä ei voida havaita ja kun kaikilla muilla kuin yhdellä komponentilla on hyvin pienet osuudet myöskään sitä ei voida havaita. Tyypillisesti (ilman teoriaa, joka voisi ehdottaa jakautumismuotoa), parametrijakaumat sopivat yksinkertaisiin likimääräisiin kuvauksiin dataan. Seokset eivät ole mikään niistä: ne vaativat liian monta parametria ja ovat liian joustavia tähän tarkoitukseen.
• @whuber: +1 Arvosta erinomainen selityksesi!
• @Lourenco Katsoin Cullenin ja Feyn kaaviota. Sininen piste tarkoittaa näytettä. Näet, että piste on lähellä Weibullin, Lognormalin ja Gamman viivoja (joka on Weibullin ja Gamman välillä). Kun olen sovittanut nämä jakaumat, verrasin sopivuustilastoja funktiolla gofstat ja AIC: llä. ’ ei ole yksimielisyyttä siitä, mikä on paras tapa määrittää ” paras ” jakelu On. Pidän graafisista menetelmistä ja AIC: stä.
• @Lourenco Tarkoitatko lognormaalia? Logistinen jakauma (” + ” -merkki) on melko kaukana havaituista tiedoista. Lognormaali olisi myös ehdokas, jota normaalisti tarkastelen ’. Tässä opetusohjelmassa olen ’ päättänyt olla näyttämättä sitä, jotta viesti olisi lyhyt. Lognormaali näyttää huonomman sovituksen verrattuna Weibull- ja Normal-jakaumiin. AIC on 537,59, ja kaaviot eivät myöskään näytä liian hyvältä.

Suunnitelmat ovat yleensä hyvä tapa saada parempi käsitys datasi ulkomuodosta. Sinun tapauksessasi suosittelisin, että piirretään empiirinen kumulatiivinen jakautumistoiminto (ecdf) teoreettisiin cdf-tiedostoihin parametreillä, jotka sait fitdistriltä ().

Tein sen kerran tietojeni vuoksi ja sisällytin myös luottamusvälit. Tässä on kuva, jonka sain käyttämällä ggplot2 ().

Musta viiva on empiirinen kumulatiivinen jakelutoiminto ja värilliset viivat ovat eri jakeluista peräisin olevia PDF-tiedostoja käyttäen parametreja, jotka sain käyttämällä Suurimman todennäköisyyden menetelmää. Voidaan helposti nähdä, että eksponentiaalinen ja normaali jakauma eivät sovi hyvin dataan, koska viivoilla on erilainen muoto kuin ecdf: llä ja linjat ovat melko kaukana ecdf: stä. Valitettavasti muut jakelut ovat melko lähellä. Mutta sanoisin, että logNormal-viiva on lähinnä mustaa viivaa. Etäisyysmittauksen (esimerkiksi MSE) avulla oletus voidaan vahvistaa.

Jos sinulla on vain kaksi kilpailevaa jakaumaa (esimerkiksi valitset ne, jotka näyttävät parhaiten sopivan kaavioon), voit käyttää Todennäköisyys-suhdetesti testataksesi, mitkä jakaumat sopivat paremmin.

Kommentit

• Tervetuloa CrossValidatediin! Vastauksestasi voi olla hyötyä, jos voisit muokata sitä sisällyttämällä (a) koodin, jota käytit grafiikan tuottamiseen, ja (b) miten luisit grafiikan.
• Mitä siellä piirretään? Onko tämä jonkinlainen eksponentiaalisuunnitelma?
• Mutta miten päätät, mikä jakelu sopii parhaiten tietoihisi? Vain graafisen kuvan perusteella en voinut ’ kertoa, sopiiko logNormal vai weibull parhaiten tietoihisi.
• Jos haluat luoda näennäissatunnaislukugeneraattorin, miksi käytä empiiristä cdf-tiedostoa? Haluatko piirtää lukuja, jotka ylittävät havaitun jakauman?
• Kun kaaviosi otetaan nimellisarvoon, näyttää siltä, että yksikään ehdokasjakelustasi ei sovi lainkaan dataan. Lisäksi ecdf: lläsi näyttää olevan vaakasuora asymptootti alle 0,03, jolla ei ole ’ järkeä, joten en ’ ole varma, että se on oikeastaan ensinnäkin ecdf.