Taller Pruebas de bondad y ajuste

1.Ejercicio de normalidad

En una fábrica de harina para arepas, se piensa que el peso de las bolsas de harina sigue una distribución normal. Se toma una muestra de pesos en gramos de 30 bolsas de manera aleatoria. Los pesos se consignan a continuación:

pesos<-c(25.74, 24.99, 26.75, 26.37, 24.44, 26.43, 25.48, 25.60, 24.75, 24.69,
         24.03, 24.88, 24.56, 23.52, 26.78, 25.79, 24.67, 25.84, 25.86, 26.10,
         24.84, 26.25, 23.88, 26.89, 26.28, 24.51, 25.43, 23.95, 25.81, 25.36)

Partiendo de la situación mencionada, se plantea la prueba de hipótesis. Sea \(X\) el peso en gramos de las bolsas de harina

\(H_0: X \sim Normal\)

\(H_1: X \nsim Normal\)

1.1 Análisis exploratorio: Pruebas gráficas

Son útiles para hacer una aproximación gráfica al comportamiento de los datos

require(car)
par(mfrow=c(1,3))
hist(pesos, xlab = "pesos de bolsas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(pesos), xlab = "pesos de bolsas", ylab = "Densidad", las=1, main = "")
qqPlot(pesos, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1] 14 24

Se puede observar un comportamiento que no parece ser del todo normal en el histograma y el diagrama de densidad. Sin embargo, los valores están dentro de las bandas de confianza del qqplot.

1.2 Estimación de los parámetros de la distribución hipotética

Los parámetros de la distribución hipotética, que en este caso es la normal, pueden estimarse a partir de la función fitdistr perteneciente al paquete MASS, que también estima parámetros de muchas otras distribuciones.

La función fitdistr tiene como argumentos: el nombre de la variable y el nombre de la distribución hipotética. Veamos

require(MASS)
Ajuste_N<-fitdistr(pesos, "normal")
Ajuste_N
##       mean          sd    
##   25.3490000    0.9165637 
##  ( 0.1673409) ( 0.1183279)

Se puede observar que la distribución hipotética implica una media estimada del peso de bolsas de harina \(\bar X= 25.349\) gramos y una desviación estándar estimada de \(s=0.916\) gramos

1.3 Prueba de bondad y ajuste

Nos permitirá concluir sobre el comportamiento de la variable analizada pesos.

Si la variable analizada es continua, se utiliza habitualmente las pruebas de Kolmogorov Smirnov y Anderson Darling.

Si se cree en principio que se tiene una distribución normal, puede usarse la prueba Shapiro Wilk.

Si la variable analizada es discreta, se utiliza la prueba Chi-cuadrado.

Como en este caso estamos tratando con una variable continua, que se cree normal, trabajaremos con las pruebas de Kolmogorov Smirnov, Anderson Darling y Shapiro Wilk.

Kolmogorov Smirnov

Se obtiene a través de la instrucción ks.test. Esta prueba tiene tres argumentos: nombre de la variable, distribución hipotética (pnorm para el caso normal) y los parámetros estimados con fitdistr.

KSN<- ks.test(pesos, "pnorm", mean =Ajuste_N$estimate[1], sd= Ajuste_N$estimate[2])
KSN
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pesos
## D = 0.098496, p-value = 0.9056
## alternative hypothesis: two-sided

En el resultado, puede observarse el valor-p\(=0.9056\) y el estadístico de prueba \(D=0.0984\).

Anderson Darling

Se obtiene a partir de la función ad.test del paquete goftest y sus argumentos son iguales a los mencionados en la prueba Kolmogorov Smirnov.

require(goftest)
ADN<- ad.test(pesos, "pnorm", mean =Ajuste_N$estimate[1], sd= Ajuste_N$estimate[2])
ADN
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 25.349, sd = 0.916563691185724
## 
## data:  pesos
## An = 0.3405, p-value = 0.9042

Shapiro Wilk

Se obtiene a partir de la instrucción shapiro.test. Solo tiene un argumento, que es la variable estudiada.

SWN<-shapiro.test(pesos)
SWN
## 
##  Shapiro-Wilk normality test
## 
## data:  pesos
## W = 0.96742, p-value = 0.4712

Análisis

Una vez obtenidos los valores-p para cada una de las pruebas, procedemos a compararlos con el valor de \(\alpha=0.05\).

Para la prueba Kolmogorov Smirnov, obtuvimos valor-p=\(0.9056 > \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\).

Para la prueba Anderson Darling, obtuvimos valor-p=\(0.9042 > \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\).

Para la prueba Shapiro Wilk, obtuvimos valor-p=\(0.4712 > \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\).

Por tanto, se concluye con un nivel de significancia del 5%, que los datos de los pesos de bolsas de harina siguen una distribución normal.

El test de Shapiro Wilk se considera uno de los más potentes para el contraste de normalidad, especialmente para muestras pequeñas.

Adicionalmente, es posible comparar la distribución acumulada empírica de las bolsas de harina con la distribución acumulada teórica normal.

par(mfrow=c(1,1))
xen <- seq(min(pesos), max(pesos), by=0.0001)
plot(xen, pnorm(xen, mean=Ajuste_N$estimate[1], sd=Ajuste_N$estimate[2]), type="l", col="red", xlab="Pesos de bolsas", ylab="pnorm(pesos, mean, sd)")
plot(ecdf(pesos), add=TRUE)

Note que el comportamiento es muy similar

2. Normalidad de proporción

Se toma 30 muestras de fusibles, cada una de tamaño 8, para compromar su calidad. Se desea analizar la proporción de fusibles defectuosos para determinar si sigue una distribución normal. Para ello, se define un nivel de significancia de \(\alpha=0.01\) .Los datos de las muestras pueden verse a continuación:

defectuosos<-c(1, 2, 1, 1, 5, 3, 2, 2, 1, 1,
               1, 4, 3, 0, 1, 3, 3, 4, 1, 0,
               2, 3, 1, 2, 2, 4, 2, 2, 3, 3)
proporcion<-defectuosos/8
proporcion
##  [1] 0.125 0.250 0.125 0.125 0.625 0.375 0.250 0.250 0.125 0.125 0.125
## [12] 0.500 0.375 0.000 0.125 0.375 0.375 0.500 0.125 0.000 0.250 0.375
## [23] 0.125 0.250 0.250 0.500 0.250 0.250 0.375 0.375

Partiendo de la situación mencionada, se plantea la prueba de hipótesis. Sea \(p\) la proporción de defectuosos

\(H_0: p \sim Normal\)

\(H_1: p \nsim Normal\)

2.1 Análisis exploratorio

require(car)
par(mfrow=c(1,3))
hist(proporcion, xlab = "proporcion de defectuosos", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(proporcion), xlab = "proporcion de defectuosos", ylab = "Densidad", las=1, main = "")
qqPlot(proporcion, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")

## [1]  5 14

2.2 Estimación de los parámetros de la distribución hipotética

require(MASS)
Ajuste_N2<-fitdistr(proporcion, "normal")
Ajuste_N2
##       mean          sd    
##   0.26250000   0.15258195 
##  (0.02785752) (0.01969824)

Observe que el comportamiento dista bastante de una distribución normal.

2.3 Prueba de bondad y ajuste

Ksn2<- ks.test(proporcion, "pnorm", mean =Ajuste_N2$estimate[1], sd= Ajuste_N2$estimate[2])
Adn2<- ad.test(proporcion, "pnorm", mean=Ajuste_N2$estimate[1], sd=Ajuste_N2$estimate[2])
Swn2<-shapiro.test(proporcion)
Ksn2
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  proporcion
## D = 0.18291, p-value = 0.268
## alternative hypothesis: two-sided
Adn2
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Normal distribution
##  with parameters mean = 0.2625, sd = 0.152581945196671
## 
## data:  proporcion
## An = 0.99404, p-value = 0.3597
Swn2
## 
##  Shapiro-Wilk normality test
## 
## data:  proporcion
## W = 0.92986, p-value = 0.04869

Análisis

Para la prueba Kolmogorov Smirnov, obtuvimos un valor-p\(=0.268> \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\)

Para la prueba Anderson Darling, obtuvimos valor-p=\(0.3597 > \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\).

Para la prueba Shapiro Wilk, obtuvimos valor-p=\(0.04869 < \alpha=0.05\), lo que sugiere que no se rechaza \(H_0\).

Las 3 pruebas sugieren que no se rechaza \(H_0\). Por tanto, la conclusión sería no rechazar la hipótesis nula. No obstante, la consideración dependerá del criterio del observador e incluye el análisis gráfico.

Adicionalmente, es posible comparar la distribución acumulada empírica de la proporción con la distribución acumulada teórica normal.

par(mfrow=c(1,1))
xen <- seq(min(proporcion), max(proporcion), by=0.0001)
plot(xen, pnorm(xen, mean=Ajuste_N2$estimate[1], sd=Ajuste_N2$estimate[2]), type="l", col="red", xlab="proporcion defectuosos", ylab="pnorm(proporcion, mean, sd)")
plot(ecdf(proporcion), add=TRUE)

Note que el comportamiento parece ser similar

3. Eventos en el tiempo

En el supermercado doña Celia, se está estudiando el comportamiento del número de personas que llegan cada hora. Se analizaron 20 horas, cuyos datos se consignan a continuación:

personas<-c(13, 14, 14, 19, 17, 14, 13,  9, 16, 16,
            13, 13, 15, 13,  7, 14, 14, 13, 20, 15)

3.1 Análisis exploratorio

Se analiza los gráficos para determinar una distribución hipotética.

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

Se observa que la forma de la distribución es asimétrica. Dado que la variable de estudio es discreta, la distribución hipotética a elegir debe corresponder a las discretas.

En este caso, la variable de interés registra un número de eventos por unidad de tiempo, por lo que se sugiere analizar el ajuste a una distribución poisson. Se muestra la respectiva prueba de hipótesis. Sea \(X\) el número de clientes que visitan Celia Express.

\(H_0:X \sim Poisson\)

\(H_1: X \nsim Poisson\)

Se realizará el análisis con un \(\alpha = 0.01\).

3.2 Estimación de los parámetros de la distribución hipotética

Para estimar los parámetros de una distribución Poisson, se requiere la función goodfit del paquete vcd. Esta función también realiza la prueba de bondad y ajuste y sus argumentos son: variable de interés, tipo de distribución y método.

Se usará el test de Chi-cuadradado a través del argumento “MinChisq”

require(vcd)
gf<-goodfit(personas, type = "poisson", method = "MinChisq")
gf$par
## $lambda
## [1] 13.60833
summary(gf)
## Warning in summary.goodfit(gf): Chi-squared approximation may be incorrect
## 
##   Goodness-of-fit test for poisson distribution
## 
##              X^2 df  P(> X^2)
## Pearson 19.30042 19 0.4377217

Análisis

En el resultado, puede observarse los grados de libertad del estadístico (df), el estadístico \(\chi ^2\) y el valor-p asociado a la prueba.

Según el resultado arrojado por la función, la tasa de llegadas es de \(13.6083\) personas por hora. De acuerdo al valor-p de la prueba Chi-cuadrado (\(0.4377\)), como este valor es superior a \(\alpha=0.01\), no se rechaza \(H_0\) y el número de personas que ingresan al banco por hora sigue una distribución Poisson.

xempp <- seq(min(personas), max(personas), by=0.0001)
plot(xempp, ppois(xempp, lambda=13.60833), type="l", col="red", xlab="Número de personas",ylab="ppois(personas, lambda)")
plot(ecdf(personas), add=TRUE)

En la gráfica se compara la función de distribución acumulada empírica del número de personas que ingresan al Celia express por hora y la teórica de la distribución poisson.

4. Tiempo entre eventos

Celia quiere ahora medir el tiempo de atención a los usuarios. Se seleccionaron 20 personas y los tiempos de atención en minutos.

tiempo<-c(3.69, 39.50,  4.43,  2.70,  9.11, 10.21, 10.44,  2.57,  5.68,  0.80,
          12.63,  2.35, 25.47, 8.07,  0.96,  0.21, 12.06, 10.79,  6.58, 13.06)

4.1 Análisis explotatorio

Procedemos a analizar el comportamiento a través de gráficos

par(mfrow=c(1,3))
hist(tiempo, xlab = "Tiempo", ylab = "Frecuencia", las=1, main = "", col = "gray")
boxplot(tiempo, col = "gray", ylab="Tiempo")
plot(density(tiempo), xlab = "Tiempo", ylab = "Densidad", las=1, main = "")

Dadas las condiciones del problema, se procede a revisar el ajuste con respecto a una distribución exponencial con un \(\alpha=0.05\). Sea \(X\) el tiempo entre llegadas a Celia Express.

\(H_0:X \sim exp\)

\(H_1: X \nsim exp\)

4.2 Estimación de parámetros

Se estima los parámetros de la distribución hipotética

Ajustex <- fitdistr(tiempo, "exponential")
Ajustex
##       rate   
##   0.11030831 
##  (0.02466569)

Observe que la tasa promedio de demora es de \(0.1103\)

4.3 Pruebas de bondad y ajuste

Usando Kolmogorov Smirnov y Anderson Darling, se espera concluir sobre la variable analizada

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

Ks
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  tiempo
## D = 0.13678, p-value = 0.8006
## alternative hypothesis: two-sided
Ad
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: exponential distribution
##  with parameter rate = 0.110308311731289
## 
## data:  tiempo
## An = 0.26651, p-value = 0.9607

Según las pruebas realizadas, no se rechaza la hipótesis nula y por tanto, se asume la distribución exponencial.

Dada la naturaleza del problema, también es posible una distribución Weibull. La correspondiente prueba de hipótesis sería:

Por ello, se procede a realizar el ajuste correspondiente

Ajustew <- fitdistr(tiempo, "weibull")
Ajustew
##      shape       scale  
##   1.0279041   9.1691599 
##  (0.1771424) (2.1014327)

Los ajustes se analizan con las pruebas habituales.

Ksw<- ks.test(tiempo, "pweibull",shape=Ajustew$estimate[1],scale=Ajustew$estimate[2])
Adw<- ad.test(tiempo, "pweibull",shape=Ajustew$estimate[1],scale=Ajustew$estimate[2])
Ksw
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  tiempo
## D = 0.13729, p-value = 0.797
## alternative hypothesis: two-sided
Adw
## 
##  Anderson-Darling test of goodness-of-fit
##  Null hypothesis: Weibull distribution
##  with parameters shape = 1.02790407066519, scale = 9.169159935727
## 
## data:  tiempo
## An = 0.26144, p-value = 0.9638

Observe que las pruebas tampoco brindan evidencia para rechazar \(H_0\), por lo que es posible asumir la distribución Weibull

Finalmente, comparamos el comportamiento de la distribución acumulada

par(mfrow=c(1,2))
xemp <- seq(min(tiempo), max(tiempo), by=0.0001)
plot(xemp, pexp(xemp, rate=Ajustex$estimate[1]), type="l", col="red", xlab="Tiempo",ylab="pexp(tiempo, rate)", main = "Distribución exponencial")
plot(ecdf(tiempo), add=TRUE)

xempw <- seq(min(tiempo), max(tiempo), by=0.0001)
plot(xempw, pweibull(xempw, shape=Ajustew$estimate[1],scale = Ajustew$estimate[2]), type="l", col="red", xlab="Tiempo",ylab="pweibull(tiempo, shape, scale)", main ="Distribución Weibull")
plot(ecdf(tiempo), add=TRUE)