Come determinare quale distribuzione si adatta meglio ai miei dati?

> 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 

Innanzitutto, ecco alcuni rapidi commenti:

• I $ p $ -valori di un Kolmovorov -Smirnov-Test (KS-Test) con parametri stimati sarà completamente sbagliato. Quindi, sfortunatamente, non puoi semplicemente adattare una distribuzione e quindi utilizzare i parametri stimati in un test di Kolmogorov-Smirnov per testare il tuo campione.
• Il tuo campione mai seguirà uno specifico distribuzione esatta. Quindi, anche se i tuoi valori $ p $ di KS-Test fossero validi e $ > 0,05 $ , significherebbe solo che non puoi “escludere che i tuoi dati seguano questa specifica distribuzione. Unaltra formulazione potrebbe essere che il tuo campione sia compatibile con una certa distribuzione. Ma la risposta alla domanda ” I miei dati seguono esattamente la distribuzione xy? ” è sempre no
• Lobiettivo qui non può essere quello di determinare con certezza quale distribuzione segue il tuo campione. Lobiettivo è ciò che @whuber (nei commenti) chiama descrizioni approssimative parsimoniose dei dati. Avere una distribuzione parametrica specifica può essere utile come modello dei dati.

---

Ma facciamo un po di esplorazione. Userò leccellente fitdistrplus che offre alcune belle funzioni per ladattamento della distribuzione. Useremo la funzione descdist per ottenere alcune idee sulle possibili distribuzioni 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) 

Ora usiamo descdist:

descdist(x, discrete = FALSE) 

La curtosi e lasimmetria al quadrato del campione sono rappresentate da un punto blu denominato ” Osservazione “. Sembra che le possibili distribuzioni includano la distribuzione Weibull, Lognormal e possibilmente Gamma.

Adattiamo una Distribuzione di Weibull e distribuzione normale:

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

Ora controlla ladattamento per il normale:

plot(fit.norm) 

E per il fit Weibull:

plot(fit.weibull) 

Entrambi sembrano buoni ma, a giudicare dal QQ-Plot, il Weibull sembra forse un po meglio, specialmente in coda. Di conseguenza, lAIC delladattamento di Weibull è inferiore rispetto alladattamento normale:

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

---

Simulazione del test di Kolmogorov-Smirnov

Userò la procedura di @Aksakal spiegata qui per simulare la statistica KS sotto il valore null.

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

LECDF delle statistiche KS simulate ha il seguente aspetto:

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

Infine, il nostro valore $ p $ che utilizza la distribuzione nulla simulata delle statistiche KS è:

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

Ciò conferma la nostra conclusione grafica che il campione è compatibile con una distribuzione di Weibull.

Come spiegato qui , possiamo utilizzare il bootstrap per aggiungere intervalli di confidenza puntuali al PDF o CDF Weibull stimato:

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

---

Adattamento automatico della distribuzione con GAMLSS

gamlss il pacchetto per R offre la possibilità di provare molte distribuzioni diverse e di selezionare best ” secondo il GAIC (il criterio di informazione Akaike generalizzato). La funzione principale è fitDist. Unopzione importante in questa funzione è il tipo di distribuzioni che vengono provate. Ad esempio, limpostazione di type = "realline" proverà tutte le distribuzioni implementate definite sullintera riga reale mentre type = "realsplus" proverà solo le distribuzioni definite sulla riga positiva reale . Unaltra opzione importante è il parametro $ k $ , che è la penalità per il GAIC. Nellesempio seguente, ho impostato il parametro $ k = 2 $ , il che significa che il ” best ” la distribuzione è selezionata secondo il classico AIC. Puoi impostare $ k $ su qualsiasi cosa tu voglia, come $ \ log (n) $ per 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 *** 

Secondo lAIC, la distribuzione di Weibull (più precisamente WEI2, una sua parametrizzazione speciale ) si adatta meglio ai dati. Lesatta parametrizzazione della distribuzione WEI2 è dettagliata in questo documento a pagina 279. Ispezioniamo ladattamento con guardando i residui in un grafico del worm (fondamentalmente un grafico QQ de-trending):

Ci aspettiamo che i residui siano vicini alla linea orizzontale centrale e il 95% di essi si trovi tra la parte superiore e curve punteggiate inferiori, che agiscono come intervalli di confidenza puntuali al 95%. In questo caso, il grafico del worm mi sembra a posto, indicando che la distribuzione di Weibull è adeguata.

Commenti

• +1 Bella analisi. Una domanda, però. La conclusione positiva sulla compatibilità con una particolare distribuzione principale (Weibull, in questo caso) consente di escludere la possibilità di una distribuzione mista ‘ s presenza? Oppure dobbiamo eseguire una corretta analisi della miscela e controllare GoF per escludere questa opzione?
• @AleksandrBlekh È impossibile avere abbastanza potenza per escludere una miscela: quando la miscela è di due distribuzioni quasi identiche non può essere rilevata e quando tutti i componenti tranne uno hanno proporzioni molto piccole non può nemmeno essere rilevato. Tipicamente (in assenza di una teoria che possa suggerire una forma distributiva), si adattano le distribuzioni parametriche al fine di ottenere descrizioni approssimative parsimoniose dei dati. Le miscele non sono nessuna di queste: richiedono troppi parametri e sono troppo flessibili allo scopo.
• @whuber: +1 Apprezzo la tua eccellente spiegazione!
• @Lourenco Ho guardato il grafico di Cullen e Fey. Il punto blu indica il nostro campione. Vedi che il punto è vicino alle linee di Weibull, Lognormal e Gamma (che è tra Weibull e Gamma). Dopo aver adattato ciascuna di queste distribuzioni, ho confrontato le statistiche sulla bontà di adattamento utilizzando la funzione gofstat e lAIC. Non vi è ‘ consenso su quale sia il modo migliore per determinare la ” migliore ” distribuzione è. Mi piacciono i metodi grafici e lAIC.
• @Lourenco Intendi il lognormal? La distribuzione logistica (il segno ” + “) è un po lontana dai dati osservati. Il lognormale sarebbe anche un candidato che ‘ guardo normalmente. Per questo tutorial, ‘ ho scelto di non mostrarlo per mantenere il post breve. Il lognormale mostra un adattamento peggiore rispetto sia alla distribuzione di Weibull che a quella normale. LAIC è 537,59 e anche i grafici non ‘ sembrano troppo belli.

I grafici sono principalmente un buon modo per avere unidea migliore dellaspetto dei dati. Nel tuo caso, consiglierei di rappresentare la funzione di distribuzione cumulativa empirica (ecdf) rispetto ai cdf teorici con i parametri ottenuti da fitdistr ().

Lho fatto una volta per i miei dati e ho incluso anche gli intervalli di confidenza. Ecco limmagine che ho ottenuto usando ggplot2 ().

La linea nera è la funzione di distribuzione cumulativa empirica e le linee colorate sono cdf di diverse distribuzioni che utilizzano parametri che ho ottenuto utilizzando il metodo di massima verosimiglianza. Si può facilmente vedere che la distribuzione esponenziale e normale non si adattano bene ai dati, perché le linee hanno una forma diversa dallecdf e le linee sono abbastanza lontane dallecdf. Purtroppo le altre distribuzioni sono abbastanza vicine. Ma direi che la linea logNormal è la più vicina alla linea nera. Utilizzando una misura della distanza (ad esempio MSE) si potrebbe convalidare lipotesi.

Se si hanno solo due distribuzioni concorrenti (ad esempio scegliendo quelle che sembrano adattarsi meglio al grafico) si potrebbe utilizzare un Test del rapporto di verosimiglianza per verificare quali distribuzioni si adattano meglio.

Commenti

• Benvenuto in CrossValidated! La tua risposta potrebbe essere più utile se potessi modificarla per includere (a) il codice che hai usato per produrre il grafico e (b) come si leggerebbe il grafico.
• Cosa viene tracciato lì? È una specie di grafico dellesponenzialità?
• Ma come decidi quale distribuzione si adatta meglio ai tuoi dati? Solo secondo il grafico non potrei ‘ t dirti se logNormal o weibull si adattano meglio ai tuoi dati.
• Se vuoi creare un generatore di numeri pseudo-casuali perché non usi il cdf empirico? Vuoi disegnare numeri che vadano oltre la tua distribuzione osservata?
• Prendendo il tuo grafico al valore nominale, sembrerebbe che nessuna delle tue distribuzioni candidate si adatta bene ai dati. Inoltre, il tuo ecdf sembra avere un asintoto orizzontale inferiore a 0,03 che ‘ non ha senso, quindi ‘ non ne sono sicuro è davvero un ecdf in primo luogo.