Ejemplo 1:

El peso que deben contener ciertas bolsas de detergente es de 750 g, con una tolerancia de ± 5 g. Se desea verificar si es razonable suponer que la distribución del peso es normal. Para ello, se toma una muestra aleatoria de 25 productos, se pesan y se obtienen los datos dados en el vector Peso:

Peso <- c(750.0, 749.3, 752.5, 748.9, 749.9, 748.6, 750.2, 748.4, 747.8, 749.3, 749.6, 749.0, 747.7, 748.3, 750.5, 750.6, 750.0, 750.4, 752.0, 750.2, 751.4, 750.9, 752.4, 751.7, 750.6)

El primer paso para determinar la distribución de los pesos es realizar un análisis exploratorio de éstos. Dado que la variable de interés es continua se pueden determinar medidas de resumen y realizar gráficos descriptivos como histogramas, y densidades. Así mismo, se puede realizar el grafico de cuantiles teóricos vs cuantiles muestrales para una distribución normal.

par(mfrow=c(1,3))
hist(Peso, xlab="Peso", ylab="Frecuencia", las=1, main="")
plot(density(Peso), xlab="Peso", ylab="Densidad", las=1, main="")
qqnorm(Peso, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")
qqline(Peso)

Los tres gráficos muestran que el comportamiento de los Pesos es simétrico y podrían seguir una distribución normal.

El segundo paso es estimar los parámetros de la distribución hipótetica a partir de la función fitdistr del paquete MASS. En este caso la distribución hipótetica es la normal. Esta función permite estimar parámetros de las siguientes distribuciones: “beta”, “cauchy”, “chi-squared”, “exponential”, “f”, “gamma”, “geometric”, “log-normal”, “lognormal”, “logistic”, “negative binomial”, “normal”, “Poisson”, “t” y “weibull”.

La función fitdistr tiene dos argumentos básicos, el primero es el nombre de los datos y el segundo es el nombre de la distribución hipótetica.

require(MASS)
ajuste <- fitdistr(Peso,"normal")
ajuste
##       mean           sd     
##   750.0080000     1.3233050 
##  (  0.2646610) (  0.1871436)

El resultado de la función es la estimación de los parámetros de la distribución hipótetica y su respectiva desviación estándar (valor entre parentesis). En este caso el peso promedio (media estimada) es de 750.008 y la desviación estándar del peso (desviación estándar estimada) es de 1.323305.

El tercer paso es usar una prueba de bondad de ajuste para probar si los pesos siguen o no una distribución normal. La primera prueba que se usará será la prueba de Kolmogorov-Smirnov. Esta prueba se obtiene a partir de la función ks.test, la cual tiene tres argumentos básicos: el nombre de la variable, el nombre de la distribución hipótetica antecedida por la letra p y los parámetros estimados.

Ks<- ks.test(Peso, "pnorm", mean =ajuste$estimate[1], sd= ajuste$estimate[2])
## Warning in ks.test(Peso, "pnorm", mean = ajuste$estimate[1], sd = ajuste
## $estimate[2]): ties should not be present for the Kolmogorov-Smirnov test
Ks
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  Peso
## D = 0.087306, p-value = 0.9911
## alternative hypothesis: two-sided

La segunda prueba que se usará será la prueba de Anderson Darling. La función asociada a esta prueba es ad.test del paquete goftest y tiene los mismos argumentos de la prueba ks.

require(goftest)
## Warning: package 'goftest' was built under R version 3.2.3
Ad<- ad.test(Peso, "pnorm", mean =ajuste$estimate[1], sd= ajuste$estimate[2])
Ad
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 750.008, sd = 1.32330495351601
## 
## data:  Peso
## An = 0.19712, p-value = 0.9913

De acuerdo al valor-p de la prueba de Kolmogorov-Smirnov (0.9911409) y de la prueba de Anderson Darling (0.9912932) se puede concluir con un nivel de significancia de 0.05 que los pesos siguen una distribución normal, puesto que son mayores a un nivel de significancia de 0.05. Cabe anotar que se usaron estas dos pruebas porque el peso de los detergentes es una variable continua. Si la variable de inetrés fuera discreta deberiamos usar la prueba Chi-cuadrado y no las dos anteriores.

La siguiente gráfica muestra la función de distribución acumulada empirica de los pesos y la teórica de la distribución normal y se observa que las dos siguen el mismo patron, lo cual corrobora la conclusión anterior.

xe <- seq(min(Peso), max(Peso), by=0.0001)
plot(xe, pnorm(xe, mean=ajuste$estimate[1], sd=ajuste$estimate[2]), type="l", col="red", xlab="x", ylab="pnorm(x, mean, sd)")
plot(ecdf(Peso), add=TRUE)

Ejemplo 2:

Se quiere determinar la distribución de las horas transcurridas hasta la falla de componentes de un computador. Para ello se toma una muestra de 20 componentes y se determinan los tiempos, los cuales están definidos en el vector Tiempos:

Tiempo <- c(3.70, 3.75, 12.18, 28.55, 29.37, 31.61, 36.78, 51.14, 108.71, 125.21, 125.35, 131.76, 158.61, 172.96, 177.12, 185.37, 212.98, 280.40, 351.28, 441.79)

El primer paso para determinar la distribución de los tiempos es realizar un histograma y la densidad. No se hará un gráfico qqnorm porque áun no se sabe cual es la distribución hipótetica. En el ejemplo 1 se hizo el gráfico qqnorm porque se sabia que la distribución hipótetica era la normal.

par(mfrow=c(1,2))
hist(Tiempo, xlab="Tiempo", ylab="Frecuencia", las=1, main="")
plot(density(Tiempo), xlab="Tiempo", ylab="Densidad", las=1, main="")

La forma de la distribución es asimetrica con sesgo positivo similar a una distribución exponencial.

El segundo paso es estimar los parámetros de la distribución hipótetica, exponencial.

require(MASS)
ajuste <- fitdistr(Tiempo,"exponential")
ajuste
##       rate    
##   0.007494510 
##  (0.001675823)

El resultado de la función es que la tasa promedio de falla de los componentes es de 0.0074945.

El tercer paso es usar las pruebas de Kolmogorov-Smirnov y de Anderson Darling, usando como argumentos el paramétro de la distribución exponencial

Ks<- ks.test(Tiempo, "pexp", rate =ajuste$estimate[1])
Ad<- ad.test(Tiempo, "pexp", rate =ajuste$estimate[1])

Ks
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  Tiempo
## D = 0.15874, p-value = 0.6381
## alternative hypothesis: two-sided
Ad
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.00749451027122633
## 
## data:  Tiempo
## An = 0.43001, p-value = 0.8167

De acuerdo al valor-p de la prueba de Kolmogorov-Smirnov (0.6380994) y de la prueba de Anderson Darling (0.8167356) se puede concluir con un nivel de significancia de 0.05 que los Tiempos siguen una distribución exponencial, puesto que son mayores a un nivel de significancia de 0.05.

La siguiente gráfica muestra la función de distribución acumulada empirica de los tiempos y la teórica de la distribución exponencial y se observa que las dos siguen el mismo patron, lo cual corrobora la conclusión anterior.

xe <- seq(min(Tiempo), max(Tiempo), by=0.0001)
plot( xe, pexp(xe, rate=ajuste$estimate[1]), type="l", col="red", xlab="x",ylab="pexp(x, rate)")
plot(ecdf(Tiempo), add=TRUE)