Se necesita verificar si es correcto suponer que el volumen de llenado (en onzas) de una máquina dispensadora de jugos sigue una distribución normal, por lo que se toman 25 botellas de forma aleatoria. Los datos del volumen de llenado obtenidos de la muestra se encuentran almacenados en el vector volumen.
Puede plantearse una prueba hipótesis de la siguiente manera:
H0: el volumen de llenado (en onzas) sigue una distribución normal.
H1: el volumen de llenado (en onzas) no sigue una distribución normal.
Nivel de significancia: 0.05 (Hipotético).
volumen <-c(8.39,12.14, 11.80,12.04,7.34,12.62,11.51,12.47,11.08,14.32,11.33,11.56,12.79,11.72,12.84,11.73,12.14,11.88,11.95,10.84,11.79,13.21,12.56,12.55,12.80)
Para verificar que el volumen de llenado sigue una distribución normal, primero debe realizarse un análisis descriptivo de los datos. Para ello, pueden realizarse gráficos como histogramas y densidades, además de analizar el gráfico de cuantiles teóricos vs cuantiles muestrales que ayudarán a dar una idea acerca de la normalidad de los datos.
require(car)
par(mfrow=c(1,3))
hist(volumen, xlab = "Volumen de llenado", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(volumen), xlab = "Volumen de llenado", ylab = "Densidad", las=1, main = "")
qqPlot(volumen, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")
## [1] 5 1
Tanto el histograma como el gráfico de densidad parecen mostrar una distribución no tan normal, mientras que en el gráfico de cuantiles teóricos vs cuantiles muestrales se observan algunos puntos por fuera de las bandas de confianza.
Una vez hecho el análisis exploratorio, deben estimarse los parámetros de la distribución hipotética, que en este caso es la normal, a partir de la función fitdistr perteneciente al paquete MASS, que también estima parámetros de muchas otras distribuciones.
Debe tenerse en cuenta que la función fitdistr tiene como argumentos el nombre de la variable donde se encuentran los datos y el nombre de la distribución hipotética.
Lo que arrojará la función, será la estimación de los parámetros de la distribución hipotética y su desviación estándar que será el valor entre paréntesis debajo de cada parámetro estimado. Por ejemplo, para el caso en cuestión, el volumen de llenado promedio estimado es de 11.816 onzas y la desviación estándar estimada es de 1.375 onzas.
require(MASS)
Ajusten<-fitdistr(volumen, "normal")
Ajusten
## mean sd
## 11.8160000 1.3755959
## ( 0.2751192) ( 0.1945386)
Una vez fueron estimados los parámetros, se usa una prueba de bondad de ajuste para finalmente concluir si la variable analizada sigue o no la distribución hipotética planteada. Debe tenerse en cuenta que si la variable en cuestión es continua, pueden utilizarse las pruebas Kolmogorov Smirnov y Anderson Darling, y si se cree en un principio que la variable sigue una distribución normal, también puede utilizarse la prueba Shapiro Wilk. Sin embargo, cuando se quiere determinar la distribución de una variable discreta, debe usarse la prueba Chi-cuadrado.
En este caso, se utilizarán las pruebas Kolmogorov Smirnov, Anderson Darling y Shapiro Wilk, debido a que la variable de interés (volumen de llenado) es continua.
La prueba Kolmogorov Smirnov se obtiene a través de la función ks.test, que consta de tres argumentos: el nombre de la variable, el nombre de la distribución hipotética con una p al inicio y los parámetros que fueron estimados anteriormente a partir de la función fitdistr. La función arrojará el valor-p asociado a la prueba y el estadístico de la misma (estadístico de Kolmogorov-Smirnov) representado por la letra D.
Ksn<- ks.test(volumen, "pnorm", mean =Ajusten$estimate[1], sd= Ajusten$estimate[2])
Ksn
##
## One-sample Kolmogorov-Smirnov test
##
## data: volumen
## D = 0.21198, p-value = 0.2112
## alternative hypothesis: two-sided
Por su parte, la prueba Anderson Darling se calcula a partir de la función ad.testdel paquete goftest y sus argumentos son iguales a los mencionados en la prueba Kolmogorov Smirnov. La función arrojará el valor-p asociado a la prueba y el estadístico de la misma (estadístico de Anderson Darling) representado por las letras An.
require(goftest)
Adn<- ad.test(volumen, "pnorm", mean =Ajusten$estimate[1], sd= Ajusten$estimate[2])
Adn
##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: Normal distribution
## with parameters mean = 11.816, sd = 1.3755958708865
##
## data: volumen
## An = 1.6222, p-value = 0.1501
Finalmente, la prueba shapiro-Wilk, simplemente tiene como argumento la variable de interés. La función arrojará el valor-p asociado a la prueba y el estadístico de la misma que es W.
Swn<-shapiro.test(volumen)
Swn
##
## Shapiro-Wilk normality test
##
## data: volumen
## W = 0.81611, p-value = 0.0004272
De acuerdo al valor-p de las pruebas Kolmogorov Smirnov (0.2112) y Anderson Darling (0.1501), el volumen de llenado sigue una distribución normal, puesto que dichos valores son mayores al nivel de significancia utilizado en esta prueba (0.05). Sin embargo, al analizar la prueba Shapiro-Wilk, se observa que el valor-p (0.00043) es inferior al nivel de significancia, por lo tanto se rechazaría la hipótesis nula y se concluiría, con un nivel de confianza del 95%, que los datos no siguen una distribución normal y le daría sentido a lo observado en el análisis exploratorio. Debe mirarse entonces, qué tan precisas son las diferentes pruebas de bondad de ajuste utilizadas y cuál tiene mayor peso. El test de Shapiro Wilk se considera uno de los más potentes para el contraste de normalidad, especialmente para muestras pequeñas.
Para ilustrar un poco más la conclusión, la siguiente gráfica muestra la distribución acumulada empírica del volumen de llenado y la teórica de una distribución normal. Note que no es muy claro que sigan el mismo patrón.
par(mfrow=c(1,1))
xen <- seq(min(volumen), max(volumen), by=0.0001)
plot(xen, pnorm(xen, mean=Ajusten$estimate[1], sd=Ajusten$estimate[2]), type="l", col="red", xlab="Volumen de llenado", ylab="pnorm(volumen, mean, sd)")
plot(ecdf(volumen), add=TRUE)
Por otro lado, para tomar acciones sobre el proceso actual de llenado, el departamento de calidad toma 15 muestras cada una de tamaño 8 y registra el número de botellas que no cumplen con las especificaciones de volumen estipuladas por el departamento, dichos datos están dados en el vector rechazadas. Sin embargo, el departamento encuentra más útil analizar la proporción de botellas que no cumplen con las especificaciones y desean verificar que dicha variable efectivamente sigue una distribución normal utilizando un nivel de significancia de 0.01.
Puede plantearse una prueba hipótesis de la siguiente manera:
H0: la proporción de botellas rechazadas sigue una distribución normal.
H1: la proporción de botellas rechazadas no sigue una distribución normal.
Nivel de significancia: 0.01.
rechazadas<-c(5,1,3,3,1,4,2,2,6,4,2,3,4,3,5)
proporcion<-rechazadas/8
proporcion
## [1] 0.625 0.125 0.375 0.375 0.125 0.500 0.250 0.250 0.750 0.500 0.250
## [12] 0.375 0.500 0.375 0.625
Como se había descrito antes, primero debe realizarse un análisis exploratorio de los datos, esto se hace con gráficos como el histograma, densidad y el gráfico de cuantiles teóricos vs cuantiles muestrales.
par(mfrow=c(1,3))
hist(proporcion, xlab = "Proporción de botellas rechazadas", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(proporcion), xlab = "Proporción de botellas rechazadas", ylab = "Densidad", las=1, main = "")
qqPlot(proporcion, xlab="Cuantiles teóricos", ylab="Cuantiles muestrales", las=1,main="")
## [1] 9 2
Los tres gráficos muestran que la proporción de botellas rechazadas tiene un comportamiento simétrico y podrían seguir una distribución normal.
Ahora, deben estimarse los parámetros de la distribución hipotética, en este caso, de nuevo es la normal a partir de la función fitdistr.
Ajusten2<-fitdistr(proporcion, "normal")
Ajusten2
## mean sd
## 0.40000000 0.17795130
## (0.04594683) (0.03248931)
Luego, deben usarse las pruebas de bondad de ajuste. Para este caso se utilizarán las mismas del numeral anterior.
Ksn2<- ks.test(proporcion, "pnorm", mean =Ajusten2$estimate[1], sd= Ajusten2$estimate[2])
Adn2<- ad.test(proporcion, "pnorm", mean=Ajusten2$estimate[1], sd=Ajusten2$estimate[2])
Swn2<-shapiro.test(proporcion)
Ksn2
##
## One-sample Kolmogorov-Smirnov test
##
## data: proporcion
## D = 0.15586, p-value = 0.8594
## alternative hypothesis: two-sided
Adn2
##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: Normal distribution
## with parameters mean = 0.4, sd = 0.177951304200522
##
## data: proporcion
## An = 0.34829, p-value = 0.8965
Swn2
##
## Shapiro-Wilk normality test
##
## data: proporcion
## W = 0.95224, p-value = 0.5604
De acuerdo al valor-p de las pruebas Kolgomorov Smirnov (0.8594), Anderson Darling (0.8965) y Shapiro-Wilk (0.5604), puede concluirse con un nivel de confianza del 99% que la proporción de botellas que no cumplen las especificaciones siguen una distribución normal, dado que los valores-p mencionados son mayores al nivel de significancia (0.01), corroborando lo observado en el análisis exploratorio.
En la siguiente gráfica de la función de distribución acumulada empírica de la proporción de botellas rechazadas y la teórica de la distribución normal, se observa el patrón muy similar que siguen ambas funciones.
par(mfrow=c(1,1))
xen2 <- seq(min(proporcion), max(proporcion), by=0.0001)
plot(xen2, pnorm(xen2, mean=Ajusten2$estimate[1], sd=Ajusten2$estimate[2]), type="l", col="red", xlab="Proporción de botellas rechazadas", ylab="pnorm(proporción, mean, sd)")
plot(ecdf(proporcion), add=TRUE)
El personal de servicio al cliente del Banco de Los Pobres, ha venido haciendo una medición del número de personas que ingresan al banco cada hora para evaluar la capacidad y eficiencia del banco en cuanto a atención. Se analizaron 20 horas y los datos asociados al número de personas que ingresan al banco por hora se encuentran definidos en el vector usuarios.
usuarios<-c(50,42,60,39,44,54,48,44,43,50,66,62,43,50,45,43,47,46,52,55)
En primer lugar, se hará un análisis exploratorio de los datos a través de gráficos como el histograma y el de densidad para observar su comportamiento y definir una distribución hipotética.
par(mfrow=c(1,2))
hist(usuarios, xlab = "Usuarios", ylab = "Frecuencia", las=1, main = "", col = "gray")
plot(density(usuarios), xlab = "Usuarios", ylab = "Densidad", las=1, main = "")
Se observa que la forma de la distribución es asimétrica con sesgo positivo. Debe tenerse en cuenta que la variable de interés en este caso es discreta, por lo tanto, la distribución hipotética a elegir debe ser de la misma clasificación. En este caso, la variable de interés registra un número de eventos por unidad de tiempo, por lo que podría ajustarse a una distribución poisson.
Puede plantearse una prueba hipótesis de la siguiente manera:
H0: el número de personas que ingresan al banco por hora sigue una distribución poisson.
H1: el número de personas que ingresan al banco por hora no sigue una distribución poisson.
Nivel de significancia: 0.1.
Posteriormente, deben estimarse los parámetros de la distribución hipotética, en este caso, la poisson. Para ello, se utilizará la función goodfitdel paquete vcd. Además de estimar el parámetro lambda de la distribución Poisson, esta función también realizará la prueba de bondad de ajuste y tiene como argumentos la variable de interés, el tipo de distribución y el método de ajuste que en este caso será “MinChisq” que es la opción para que realice el test Chi-cuadrado. El resultado serán los grados de libertad del estadístico, el estadístico X2 y el valor-p asociado a la prueba. Las pruebas de bondad de ajuste utilizadas anteriormente (Kolmogorov Smirnov y Anderson Darling) no se utilizarán para este caso debido a que la variable de interés es discreta.
require(vcd)
gf<-goodfit(usuarios, type = "poisson", method = "MinChisq")
gf$par
## $lambda
## [1] 49.63808
summary(gf)
##
## Goodness-of-fit test for poisson distribution
##
## X^2 df P(> X^2)
## Pearson 25.89737 65 0.9999964
Según el resultado arrojado por la función, la tasa de llegadas es de 49.63808 personas por hora. De acuerdo al valor-p de la prueba Chi-cuadrado (0.9999964), la hipótesis nula no se rechaza debido a que dichos valores son mayores que el nivel de significancia (0.05), por lo tanto, el número de personas que ingresan al banco por hora sigue una distribución Poisson.
La siguiente gráfica de la función de distribución acumulada empírica del número de personas que ingresan al banco por hora y la teórica de la distribución poisson siguen el mismo patrón apoyando la conclusión dada anteriomente.
xempp <- seq(min(usuarios), max(usuarios), by=0.0001)
plot(xempp, ppois(xempp, lambda=49.63808), type="l", col="red", xlab="Número de usuarios",ylab="ppois(usuarios, lambda)")
plot(ecdf(usuarios), add=TRUE)
Para mejorar la atención a los usuarios, el personal de servicio al cliente del banco también midió el tiempo que demora un usuario desde que ingresa al banco hasta que sale, esto con el objetivo de mejorar la eficiencia en los tiempos de atención y disminuir las filas que se forman en el banco. Para ello, aleatoriamente se seleccionaron 20 personas y los tiempos (en minutos) asociados a cada una se encuentran en el vector tiempo.
tiempo<-c(9.2,10.5,10.8,12.3,14.8,15.6,16.1,19.5,22.4,24.7,25.3,26.5,29.9,30.3,30.7,35.6,46.5,50.4,68.2,90.6)
Como se ha mencionado con anterioridad, para determinar la distribución hipotética a evaluar, se utilizarán gráficos como histogramas, boxplot y densidad.
par(mfrow=c(1,3))
hist(tiempo, xlab = "Tiempo", ylab = "Frecuencia", las=1, main = "", ylim = c(0,8), col = "gray")
boxplot(tiempo, col = "gray", ylab="Tiempo")
plot(density(tiempo), xlab = "Tiempo", ylab = "Densidad", las=1, main = "")
La distribución tiene forma asimétrica con sesgo positivo muy similar a la distribución exponencial.
Puede plantearse una prueba hipótesis de la siguiente manera:
H0: el tiempo que demora un usuario en el banco sigue una distribución exponencial.
H1: el tiempo que demora un usuario en el banco sigue una distribución exponencial.
Nivel de significancia: 0.05.
Ahora deben estimarse los parámetros de la distribución hipotética, en este caso, la exponencial.
Ajustex <- fitdistr(tiempo, "exponential")
Ajustex
## rate
## 0.033904052
## (0.007581176)
Según la función fitdistr, la tasa promedio de demora de los usuarios es de 0.0339.
A partir de las pruebas Kolmogorov Smirnov y Anderson Darling, finalmente puede concluirse si los datos siguen una distribución exponencial.
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.26796, p-value = 0.09301
## alternative hypothesis: two-sided
Ad
##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: exponential distribution
## with parameter rate = 0.0339040515341583
##
## data: tiempo
## An = 1.6093, p-value = 0.1527
De acuerdo al valor-p de las pruebas Kolmogorov Smirnov (0.09301) y Anderson Darling (0.1527), la hipótesis nula no se rechaza debido a que dichos valores son mayores que el nivel de significancia (0.05), por lo tanto, los tiempos de demora siguen una distribución exponencial.
Sin embargo, también pudo haberse ajustado la variable de interés a una distribución Weibull, ya que los datos también tienen forma similar a dicha distribución. Por otro lado, hay que recordar ue la distribución exponencial es un caso particular de la Weibull.
Puede plantearse una prueba hipótesis de la siguiente manera:
H0: el tiempo que demora un usuario en el banco sigue una distribución Weibull.
H1: el tiempo que demora un usuario en el banco sigue una distribución Weibull.
Nivel de significancia: 0.05.
Ahora deben estimarse los parámetros de la distribución hipotética, en este caso, la weibull.
Ajustew <- fitdistr(tiempo, "weibull")
Ajustew
## shape scale
## 1.6042319 33.2093572
## ( 0.2600694) ( 4.9167994)
A partir de las pruebas Kolmogorov Smirnov y Anderson Darling, finalmente puede concluirse si los datos siguen una distribución Weibull.
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.16413, p-value = 0.5972
## alternative hypothesis: two-sided
Adw
##
## Anderson-Darling test of goodness-of-fit
## Null hypothesis: Weibull distribution
## with parameters shape = 1.60423192886976, scale =
## 33.2093571772492
##
## data: tiempo
## An = 0.51385, p-value = 0.7306
De acuerdo al valor-p de las pruebas Kolmogorov Smirnov (0.5972) y Anderson Darling (0.7306), la hipótesis nula no se rechaza debido a que dichos valores son mayores que el nivel de significancia (0.05), por lo tanto, los tiempos de demora siguen una distribución weibull.
Obsérvese que el valor-p en ambas pruebas, para este ejrecicio, es mayor que en el caso de la distribución exponencial, por lo que los datos podrían tener un mejor ajuste bajo la distribución Weibull.
Las siguientes gráficas de la distribución acumulada empírica de los tiempos y la teórica de las distribuciones exponencial y Weibull dan una mejor idea acerca de cuál distribución se ajusta mejor a los datos.
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)
Efectivamente, la distribución Weibull se ajusta mejor a los datos de los tiempos al seguir un patrón mucho más parecido a la función de distribución acumulada teórica.