Cum să stabilesc ce distribuție se potrivește cel mai bine cu datele mele?

> 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 

Mai întâi, iată câteva comentarii rapide:

• $ p $ -valorile unui Kolmovorov -Smirnov-Test (KS-Test) cu parametri estimate va fi destul de greșit. Așadar, din păcate, nu puteți doar să încadrați o distribuție și apoi să utilizați parametrii estimați într-un test Kolmogorov-Smirnov pentru a testa eșantionul.
• Eșantionul dvs. niciodată nu va urma distribuție exact. Deci, chiar dacă $ p $ -valorile din testul KS ar fi valide și $ 0.05 $ , ar însemna doar că nu puteți exclude ca datele dvs. să respecte această distribuție specifică. O altă formulare ar fi că eșantionul dvs. este compatibil cu o anumită distribuție. Dar răspunsul la întrebarea ” Datele mele respectă exact distribuția xy? ” este întotdeauna nu.
• Scopul de aici nu poate fi de a determina cu certitudine ce distribuție urmează eșantionul dvs. Scopul este ceea ce @whuber (în comentarii) numește descrieri aproximative parsimoniose ale datelor. Având o distribuție parametrică specifică poate fi util ca model al datelor.

---

Dar să facem o explorare. Voi folosi excelentul fitdistrplus pachet care oferă câteva funcții frumoase pentru montarea distribuției. Vom folosi funcția descdist pentru a câștiga câteva idei despre distribuții posibile candidate.

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) 

Acum permiteți utilizarea descdist:

descdist(x, discrete = FALSE) 

Curtoza și asimetria pătrată a probei dvs. sunt reprezentate ca punct albastru numit ” Observație „. Se pare că distribuțiile posibile includ distribuția Weibull, Lognormal și posibil distribuția Gamma.

Distribuție Weibull și o distribuție normală:

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

Acum inspectați potrivirea pentru normal:

plot(fit.norm) 

Și pentru potrivirea Weibull:

plot(fit.weibull) 

Ambele arată bine, dar judecate după QQ-Plot, Weibull poate arăta puțin mai bine, mai ales la cozi. În mod corespunzător, AIC al ajustării Weibull este mai mic comparativ cu ajustarea normală:

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

---

Simularea testului Kolmogorov-Smirnov

Voi folosi procedura lui @Aksakal explicată aici pentru a simula statistica KS sub nul.

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 al statisticilor KS simulate arată după cum urmează:

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

În cele din urmă, valoarea noastră $ p $ utilizând distribuția nulă simulată din statisticile KS este:

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

Acest lucru confirmă concluzia noastră grafică că eșantionul este compatibil cu o distribuție Weibull.

După cum sa explicat aici , putem utiliza bootstrapping-ul pentru a adăuga intervale de încredere punctuale la PDF-ul Weibull sau CDF estimat:

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

---

Montaj automat de distribuție cu GAMLSS

gamlss pachetul pentru R oferă posibilitatea de a încerca mai multe distribuții diferite și de a selecta cel mai bun ” conform GAIC (criteriul de informare generalizat Akaike). Funcția principală este fitDist. O opțiune importantă în această funcție este tipul distribuțiilor care sunt încercate. De exemplu, setarea type = "realline" va încerca toate distribuțiile implementate definite pe întreaga linie reală în timp ce type = "realsplus" va încerca doar distribuțiile definite pe linia reală pozitivă . O altă opțiune importantă este parametrul $ k $ , care reprezintă penalizarea pentru GAIC. În exemplul de mai jos, am setat parametrul $ k = 2 $ ceea ce înseamnă că ” cel mai bun ” distribuția este selectată în conformitate cu clasicul AIC. Puteți seta $ k $ la orice doriți, cum ar fi $ \ log (n) $ pentru 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 *** 

Conform AIC, distribuția Weibull (mai precis WEI2, o parametrizare specială a acesteia ) se potrivește cel mai bine cu datele. Parametrizarea exactă a distribuției WEI2 este detaliată în acest document la pagina 279. Să inspectăm potrivirea prin privind reziduurile dintr-un grafic de vierme (practic un grafic QQ de-trended):

Ne așteptăm ca reziduurile să fie aproape de linia orizontală mijlocie și 95% dintre ele să se afle între partea superioară și curbe punctate inferioare, care acționează ca intervale de încredere punctuale de 95%. În acest caz, graficul viermelui mi se pare bine, indicând că distribuția Weibull este o potrivire adecvată.

Comentarii

+1 Analiză frumoasă. O întrebare, totuși. Oare concluzia pozitivă a compatibilității cu o anumită distribuție majoră (Weibull, în acest caz) permite excluderea posibilității unei distribuții de amestec 74282493e5 „>

Graficele sunt în mare parte o modalitate bună de a vă face o idee mai bună despre aspectul datelor dvs. În cazul dvs., aș recomanda reprezentarea grafică a funcție empirică de distribuție cumulativă (ecdf) împotriva cdf-urilor teoretice cu parametrii pe care i-ați primit de la fitdistr ().

Am făcut asta o dată pentru datele mele și am inclus și intervalele de încredere. Iată imaginea pe care am obținut-o folosind ggplot2 ().

Linia neagră este funcția de distribuție cumulativă empirică iar liniile colorate sunt cdf-uri din diferite distribuții folosind parametrii pe care i-am obținut folosind metoda Maximum Likelihood. Se poate vedea cu ușurință că distribuția exponențială și normală nu se potrivește bine cu datele, deoarece liniile au o formă diferită de ecdf și liniile sunt destul de departe de ecdf. Din păcate celelalte distribuții sunt destul de apropiate. Dar aș spune că linia logNormal este cea mai apropiată de linia neagră. Folosind o măsură a distanței (de exemplu MSE) s-ar putea valida presupunerea.

Dacă aveți doar două distribuții concurente (de exemplu alegerea celor care par să se potrivească cel mai bine în complot) puteți utiliza un Likelihood-Ratio-Test pentru a testa distribuțiile care se potrivesc mai bine.

Comentarii

• Bine ați venit la CrossValidated! Răspunsul dvs. ar putea fi mai util dacă l-ați putea edita pentru a include (a) codul pe care l-ați folosit pentru a produce imaginea și (b) modul în care s-ar citi graficul.
• Ce se trasează acolo? Este acesta un fel de complot de exponențialitate?
• Dar cum decideți ce distribuție se potrivește cel mai bine cu datele dvs.? Numai conform graficului nu ți-am putut ‘ să-ți spun dacă logNormal sau weibull se potrivesc cel mai bine datelor tale.
• Dacă vrei să creezi un generator de numere pseudo-aleatorii de ce nu folosiți CD-ul empiric? Doriți să atrageți numere care depășesc distribuția observată?
• Luând graficul la valoarea nominală, se pare că niciunul din distribuțiile candidate nu se potrivește deloc bine cu datele. De asemenea, ecdf-ul dvs. pare să aibă o asimptotă orizontală la mai puțin de 0,03 ceea ce nu are ‘ t sens, așa că nu ‘ nu sunt sigur că este într-adevăr un ecdf în primul rând.