Introducción

El conjunto de datos BATTERY incluido en el paquete PASWR2 contiene 100 observaciones de 2 variables correspondientes a la duración de dos tipos de baterías A y B (en horas). El conjunto de datos es un data.frame con las columnas lifetime y facility. Para realizar esta práctica, carga primero el conjunto de datos en tu espacio de trabajo, por ejemplo:

library(PASWR2)
datos <- BATTERY

Fíjate que tienes que tener instalado el paquete PASWR2 para poder acceder a este conjunto de datos. La variable de interés es lifetime, pero como sabemos que los datos se refieren a dos tipos distintos de baterías, posiblemente nos interese separarlos. En esta práctica vamos a realizar cálculo de probabilidades basados en este conjunto de datos para que se vea una aplicación, aunque tengamos que hacer uso de algún concepto de inferencia.

library(ggplot2)
library(dplyr)
library(magrittr)

Actividad 1

datos %>%
  select(lifetime) %>%
  ggplot(aes(x=lifetime)) +
  geom_histogram(aes(fill=..count..), binwidth=0.5, col='black') +
  xlim(170, 210) +
  xlab('lifetime (hours)')

Como se ve en la figura anterior, existe un hueco entre los valores 185 (aproximadamente) y 195 (aproximadamente). Este hueco forma dos subconjuntos de datos y esto quiere decir que será interesante analizar por separado ambos conjuntos. Podríamos decir que hay dos tipos de baterías en función de su tiempo de vida.

batteries_type_b <- datos %>%
                      filter(lifetime < 190)

batteries_type_a <- datos %>%
                      filter(lifetime > 190)

Aquellas baterías con tiempo de vida menor a 190 horas serán Baterías de tipo B y aquellas con un tiempo de vida mayor a 190 horas serán las Baterías de tipo A.

He decidido hacer esta diferenciación con el valor 190 porque no hay ninguna batería con tiempo de vida 190 horas y aquellos valores inferior a este pertenecerán a un tipo de batería y los mayores a otro tipo de batería.

batteries_type_a %>%
  select(lifetime) %>%
  ggplot(aes(x=lifetime)) +
  geom_histogram(aes(y = ..density.., fill=..count..), binwidth=0.5, col='black') +
  xlim(190, 210) +
  geom_density(alpha=.2, fill="#FF6666") +
  xlab('lifetime (hours)')

summary(batteries_type_a)
##     lifetime     facility
##  Min.   :194.1   A:50    
##  1st Qu.:198.8   B: 0    
##  Median :200.3           
##  Mean   :200.5           
##  3rd Qu.:202.9           
##  Max.   :206.6
batteries_type_b %>%
  select(lifetime) %>%
  ggplot(aes(x=lifetime)) +
  geom_histogram(aes(y = ..density.., fill=..count..), binwidth=0.25, col='black') +
  xlim(170, 185) + 
  geom_density(alpha=.2, fill="#FF6666") +
  xlab('lifetime (hours)')

summary(batteries_type_b)
##     lifetime     facility
##  Min.   :174.2   A: 0    
##  1st Qu.:178.5   B:50    
##  Median :179.6           
##  Mean   :179.7           
##  3rd Qu.:181.1           
##  Max.   :183.6

Es difícil dar un pronóstico ya que la anchura de las barras del histograma que he seleccionado no permite apreciar bien si se forma una gráfica típica de una distribución normal. Si dibujamos por encima del histograma el diagrama de densidad de probabilidad, aunque parece que se asemeja a una distribución normal ya que es más o menos simétrica y centrada en la media, el dibujo no es exactamente igual a una campana de Gauss.

Será necesario realizar test de normalidad o gráficos que nos ayuden a comprobar esto.

library(scales)

qqnorm(batteries_type_a$lifetime, pch = 20, col = alpha("red4", 0.5), las = 1)
grid()
qqline(batteries_type_a$lifetime, lwd = 2)

qqnorm(batteries_type_b$lifetime, pch = 20, col = alpha("red4", 0.5), las = 1)
grid()
qqline(batteries_type_b$lifetime, lwd = 2)

shapiro.test(batteries_type_a$lifetime)
## 
##  Shapiro-Wilk normality test
## 
## data:  batteries_type_a$lifetime
## W = 0.9848, p-value = 0.7632
shapiro.test(batteries_type_b$lifetime)
## 
##  Shapiro-Wilk normality test
## 
## data:  batteries_type_b$lifetime
## W = 0.98395, p-value = 0.7256

Con el test de normalidad Shapiro-Wilk para un valor \(\alpha\) de 0.05, como la hipótesis nula es que la población sigue una distribución normal y nuestros p-valor son mayores que \(\alpha\), no podemos rechazar que la distribución sea normal.

Gracias al gráfico Q-Q podemos ver que ambos conjuntos de datos se asemejan bastante a una distribución normal, especialmente batteries_type_b (los puntos rojos, que son las muestras de cuantiles, se aproximan mucho a los cuantiles teóricos, la recta negra).

Actividad 2

Ahora que sabemos que nuestros datos siguen aproximadamente una distribución normal, tendríamos que estimar sus parámetros µ y σ. A partir de ahí, podemos realizar cálculo de probabilidades de la normal.

library(fitdistrplus)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: survival
print('Estimaciones para las baterías de tipo A')
## [1] "Estimaciones para las baterías de tipo A"
fitdist(batteries_type_a$lifetime, "norm")
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##       estimate Std. Error
## mean 200.50866  0.3844088
## sd     2.71818  0.2718179
print('Estimaciones para las baterías de tipo B')
## [1] "Estimaciones para las baterías de tipo B"
fitdist(batteries_type_b$lifetime, "norm")
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##        estimate Std. Error
## mean 179.680524  0.2918968
## sd     2.064022  0.2064020

Para realizar las estimaciones de la media y de la desviación típica he hecho uso del paquete fitdistrplus. Con el parámetro ‘norm’ se indica que se quieren estimar los parámetros de una distribución normal (\(\mu\) y \(\sigma\)). El método empleado es el de máxima verosimilitud.

Para las baterías de tipo A, los valores estimados de los parámetros son \(\mu\) = 179.680524 y \(\sigma\) = 2.064022.

Para las baterías de tipo B, los valores estimados de los parámetros son \(\mu\) = 200.50866 y \(\sigma\) = 2.71818.

pnorm(q = 210, mean = 200.50866, sd = 2.71818, lower.tail = FALSE)
## [1] 0.000239889

La probabilidad de que haya una batería tomada al azar del tipo A y que dure más de 210 horas es del 0%.

pnorm(q = 175, mean = 179.680524, sd = 2.064022)
## [1] 0.01167462

La probabilidad de que haya una batería tomada al azar del tipo B y que dure menos de 175 horas es del 1.167462%.

qnorm(p = 0.03, mean = 179.680524, sd = 2.064022)
## [1] 175.7985

La duración máxima del 3% de las pilas del tipo B que duran menos es de 195.3963 horas.

Actividad 3

Vamos a centrarnos ahora en las baterías de tipo B. Supongamos que una duración por debajo de 175 horas no es aceptable para el usuario de la batería. En la actividad anterior hemos calculado la probabilidad p de que esto suceda. Entonces, si tomamos una batería del tipo B al azar y comprobamos si dura menos de 175 horas, estamos realizando un experimento de Bernoulli con probabilidad p.

La distribución binomial cuenta el número de éxitos en una muestra de tamaño n. En nuestro caso, el tamaño de la muestra es n = 10 y la probabilidad de éxito del suceso (que una batería tomada al azar del tipo B dure menos de 175 horas) es de p = 0.01167462.

pbinom(q = 0, size = 10, prob = 0.01167462)
## [1] 0.8892001

La probabilidad de que haya 0 baterías defectuosas en un lote de 10 es de 88.92001%.

La distribución geométrica es una distribución de probabilidad del número de experimentos necesarios hasta obtener el primer éxito en una serie de pruebas independientes de Bernoulli con probabilidad de éxito p. En este caso, p = 0.01167462, siendo el caso de éxito que una batería tomada al azar del tipo B dure menos de 175 horas.

dgeom(x = 4, prob = 0.01167462)
## [1] 0.01113891

La probabilidad de que la batería producida en quinto lugar sea la primera defectuosa es del 1.113891%.

La distribución hipergeométrica mide la probabilidad de obtener x (\(0 \leq x \leq d\)) elementos del grupo A, de una submuestra sin reemplazamiento de n elementos pertenecientes a una muestra de N elementos con d elementos de clase A.

En este caso tenemos que d = 3, n = 5, N = 20. En la función dhyper m = d, n = N - d y k = n.

phyper(q = 0, m = 3, n = 17, k = 5, lower.tail = FALSE)
## [1] 0.6008772

La probabilidad de que al tomar una muestra sin reposición de 5 baterías al menos una sea defectuosa es del 60.08772%.

Actividad 4

Seguimos con las baterías de tipo B, pero en vez de hacer experimentos de bernoulli queremos estudiar el número de baterías defectuosas fabricadas cada día. Supongamos que se fabrican 1000 baterías cada día. Entonces, cada día en promedio se estarán produciendo aproximadamente 1000 × p baterías, y el número de baterías defectuosas por día sigue una distribución de Poisson. Tomemos 12 como ese promedio de baterías defectuosas cada día.

Ahora estamos trabajando con la distribución de Poisson, muy útil para conocer la probabilidad con que ocurren eventos por unidad de tiempo. En este ejercicio λ = 12.

ppois(q = 20, lambda = 12, lower.tail = FALSE)
## [1] 0.01159774

La probabilidad de que en un día se produzcan más de 20 baterías defectuosas es del 1.159774%.

ppois(q = 0, lambda = 12)
## [1] 6.144212e-06

La probabilidad de que un día no salga ninguna batería defectuosa de la fábrica es del 0.0006144212%.

Siguiendo la suposición anterior, como cada día se fabrican 12 baterías defectuosas en promedio (\(\lambda = 12\)), el número promedio de baterías defectuosas a la semana será de 60 (\(\lambda = 60\)), ya que la fábrica funciona sólo de lunes a viernes, es decir, 5 días a la semana (12 baterías defectuosas por día * 5 días que funciona la fábrica por semana).

Esto se puede realizar porque la suma de variables aleatorias de Poisson da como resultado otra variable aleatoria de Poisson cuyo valor es la suma de las variables aleatorias de Poisson.

Se tienen en cuenta 5 días de la semana únicamente, porque hay 2 días que la fábrica no está en funcionamiento. Para cada día que la fábrica funciona, este proceso sigue una distribución de Poisson con parámetro \(\lambda = 12\).

Por la propiedad vista anteriormente, podemos sumar los \(\lambda\) de cada día, resultando: \(\lambda_{lunes} + \lambda_{martes} + \lambda_{miércoles} + \lambda_{jueves} + \lambda_{viernes} = 12 + 12 + 12 + 12 + 12 = 60\) baterías defectuosas por día.

Actividad 5

El departamento de I+D de la empresa que fabrica las baterías tipo B está investigando nuevos materiales y métodos para mejorar la vida útil de las baterías. En particular, quieren llegar a diseñar una batería cuya duración siga una distribución de Weibull con parámetros a = 100 y b = 185.

En esta ocasión emplearemos la distribución de Weibull porque es interesante para modelar tiempos de vida de mecanismos o elementos.

weibull_simulation <- data.frame(rweibull(n = (5 * 1000), shape = 100, scale = 185))
colnames(weibull_simulation) <- c('values')

weibull_simulation %>%
  ggplot(aes(x=values)) +
  geom_histogram(aes(fill=..count..), binwidth=0.25, col='black') +
  xlab('lifetime (hours)')

weibull_simulation %>%
  ggplot(aes(y=values, x='Baterías de tipo B')) +
  geom_boxplot() +
  xlab('') +
  ylab('lifetime (hours)')

Con la función rweibull generamos n valores aleatorios a partir de una función de distribución Weibull con parámetros shape (representa el parámetro a) y scale (representa el parámetro b).

summary(weibull_simulation)
##      values     
##  Min.   :168.9  
##  1st Qu.:182.7  
##  Median :184.4  
##  Mean   :184.0  
##  3rd Qu.:185.6  
##  Max.   :188.8

Este es el output del tercer apartado del primer ejercicio: lifetime facility Min. :174.2 A: 0
1st Qu.:178.5 B:50
Median :179.6
Mean :179.7
3rd Qu.:181.1
Max. :183.6

La duración media de las baterías mejora desde las 179.7 horas hasta las 184.0 horas: supone una mejora de 4.3 horas de media, un 2.3929% más de duración.

pweibull(q = 175, shape = 100, scale = 185)
## [1] 0.003852956

La probabilidad de que haya baterías defectuosas es del 0.3852956%. Es sin duda un porcentaje notablemente inferior al anterior, que era del 1.167462%. Supone una reducción de baterías defectuosas del 0.7821664% y una disminución del 66.99716% respecto a la probabilidad anterior.

Conclusiones y comentarios

Me ha resultado una práctica amena y útil para refrescar la teoría sobre funciones de distribución de probabilidad. Esta práctica se ha realizado usando la sintáxis RMarkdown, además de las librerías magrittr, dplyr y knittr y el IDE RStudio. Como en cada práctica, he utilizado la herramienta de control de versiones Git y el repositorio GitHub de esta práctica se encuentra en https://github.com/jvicentem/probability-distributions-assignment.