Practica de probabilidad_modelos de distribucion (1).Rmd

---
title: "Practica final: Modelos de distribución"
author: 
- Mónica Alexa
- Raúl Salazar de Torres
- Jorge Fernández Hernández
- Carlos Sánchez Vega
date: "19 de noviembre de 2017"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### 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.


```{r, warning=FALSE,message=FALSE}
library(PASWR2)
library(dplyr)
library(ggplot2)
library(scales)
library(gridExtra)
library(evd)

datos = BATTERY
```


### Actividad 1

##### 1.1 Realiza un histograma de todas las filas de la variable `lifetime` y comprueba que efectivamente nos interesa separar los datos.

```{r,fig.align='center'}
ggplot(BATTERY, aes(lifetime)) + 
  geom_histogram(aes(fill = ..count..), binwidth = 0.7)
```



**Los datos se distribuyen en dos subconjuntos ya que aparecen separados en el histograma. Por estas razones conviene separar los datos.**


##### 1.2 Crea dos conjuntos de datos diferentes para los dos tipos de baterías, por ejemplo datosA y datosB

```{r}
datosA = BATTERY %>%
  filter(facility == 'A')

datosB = BATTERY %>%
  filter(facility == 'B')
```


##### 1.3 Realiza ahora un histograma de cada uno de los tipos y comenta si te parece que los datos siguen una distribucion normal

```{r,fig.align='center'}
datosA_plot = ggplot(datosA, aes(lifetime)) + 
  geom_histogram(aes(fill = ..count..),binwidth = 0.5) + 
  ggtitle('Histograma datos A') + 
  theme(plot.title = element_text(hjust = 0.5))

datosB_plot = ggplot(datosB, aes(lifetime)) + 
  geom_histogram(aes(fill = ..count..),binwidth = 0.5) + 
  ggtitle('Histograma datos B') + 
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(datosA_plot, datosB_plot, ncol = 2, nrow = 1)
```

*A primera vista es dificil determinar si los datos siguen una distribución normal. Vamos a añadir una curva de densidad a los datos:*

```{r,fig.align='center'}
datosA_plot = ggplot(datosA, aes(lifetime)) + 
  geom_histogram(aes(fill = ..count..),binwidth = 0.5, fill='grey') + 
  geom_density(aes(y=0.5 * ..count..), alpha = 0.3, fill = '#FF6666') + 
  ggtitle('Histograma datos A') + 
  theme(plot.title = element_text(hjust = 0.5))

datosB_plot = ggplot(datosB, aes(lifetime)) + 
  geom_histogram(aes(fill = ..count..),binwidth = 0.5, fill='grey') + 
  geom_density(aes(y=0.5 * ..count..), alpha = 0.3, fill = '#FF6666') + 
  ggtitle('Histograma datos B') + 
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(datosA_plot, datosB_plot, ncol = 2, nrow = 1)
```

**Aún con la curva de densidad, parece que tiene forma acampanada y parece simetrica no obstante no se aprecia claramente la forma de campana de Gauss. Por lo tanto sería adecuado realizar alguno de los test de normalidad.**



##### 1.4 Confirma tus conclusiones con alguna/s de las herramientas vistas en clase (test de normalidad, gráfico Quantil-Quantil... )

*Según el gráfico Quantil-Quantil si se trata de una distribución normal, los puntos deberían coincidir aproximadamente con la línea. Aplicamos el gráfico sobre los dos conjuntos de datos:*

**Gráfico Quantil-Quantil DatosA**
```{r,fig.align='center'}
qqnorm(datosA$lifetime, pch = 20, col = alpha("red4", 0.5), las = 1) 
grid() 
qqline(datosA$lifetime, lwd = 2)

```

**Gráfico Quantil-Quantil DatosB**

```{r,fig.align='center'}
qqnorm(datosB$lifetime, pch = 20, col = alpha("red4", 0.5), las = 1)
grid()
qqline(datosB$lifetime, lwd = 2)

```



*Observamos en el caso de los Datos A que los puntos se alejan ligeramente de la línea en varios segmentos. Por lo tanto vamos a realizar un contraste de normalidad y aplicamos el test Shapiro-Wilk. Establecemos con hipotesis nula que los datos son normales. Si el p-valor toma un valor pequeño, rechazamos la hipotesis nula.*


**Test Shapiro-Wilk DatosA**

```{r}
shapiro.test(datosA$lifetime)
```

**Test Shapiro-Wilk DatosB**
```{r}
shapiro.test(datosB$lifetime)
```
**El test Shapiro-Wilk nos muestra que para una confianza de 95%, el p-valor debería ser menor de 0.05 para aceptar la hipotesis nula. No obstante, en los dos casos no podemos rechazar que las variables sigan una distribución normal.**    



### 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.  

##### 2.1 Realiza una estimación puntual de la media y la desviación típica de la población de cada tipo de baterías (media y la desviación tipica)


**Media y desviación típica de datosA**
```{r}
mean(datosA$lifetime)
sd(datosA$lifetime)
```

**Representación gráfica de la média para los datos de A**
```{r}
datosA_plot + geom_vline(aes(xintercept=mean(datosA$lifetime)), color="blue", linetype="dashed", size=1) 
```

**Media y desviación típica de datosB**
```{r}
mean(datosB$lifetime)
sd(datosB$lifetime)
```

**Representación gráfica de la média para los datos de B**
```{r}
datosB_plot + geom_vline(aes(xintercept=mean(datosB$lifetime)), color="blue", linetype="dashed", size=1) 
```

##### 2.2 Calcula la probabilidad de que una batería tomada al azar del tipo A dure más de 210 horas.

```{r}
pnorm(q = 210, mean = mean(datosA$lifetime), sd = sd(datosA$lifetime), lower.tail = FALSE)

# También podemos calcularlo así:
1-pnorm(q = 210,mean = mean(datosA$lifetime), sd = sd(datosA$lifetime))
```



**La probabilidad de que una batería tomada al azar del tipo A dure más de 210 horas es practicamente nula.**



##### 2.3 Calcula la probabilidad de que una batería tomada al azar del tipo B dure menos de 175 horas.


```{r}
dataBlt175 = pnorm(q = 175, mean = mean(datosB$lifetime), sd = sd(datosB$lifetime))
dataBlt175
```


**La probabilidad de que una batería tomada al azar del tipo B dure menos de 175 horas es de 1.23%**


##### 2.4 Encuentra cuál es la duración máxima del 3% de las pilas del tipo B que duran menos (ayuda: esto es equivalente a encontrar el cuantil 0.03 de la distribución).


```{r}
qnorm(p = 0.03, mean = mean(datosB$lifetime), sd = sd(datosB$lifetime))
```

**La duración máxima del 3% de las pilas del tipo B que duran menos  es 175.7591 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.


##### 3.1 Calcula la probabilidad de que en un lote de 10 baterías, no haya ninguna defectuosa (ayuda: distribución binomial).

```{r}
pbinom(q = 0, size = 10, prob = dataBlt175)

# También se puede calcular asi:

dbinom(0, size = 10, prob = dataBlt175)
```

**La probabilidad de que en un lote de 10 baterías, no haya ninguna defectuosa es de 88.28%**


##### 3.2 Imagina que las baterías se fabrican en serie e independientemente. ¿Cuál es la probabilidad de que la batería producida en quinto lugar sea la primera defectuosa? (ayuda: distribución geométrica.)

```{r}
dgeom(x = 4, prob = dataBlt175)
```

**La probabilidad de que la batería producida en quinto lugar sea la primera defectuosa es de 1.17%**


##### 3.3 Supongamos que en una caja de 20 baterías van 3 defectuosas. ¿Cuál es la probabilidad de que al tomar una muestra sin reposición de 5 baterías al menos una sea defectuosa? (ayuda: distribución hipergeométrica)


```{r}
phyper(q =0, m = 3, k = 5, n = 20-3, lower.tail = FALSE)
```


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


### 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. (ayuda: repasa qué modelo de distribución modeliza estos recuentos de eventos raros con una tasa media por unidad de tiempo)


##### 4.1 ¿Cuál es la probabilidad de que un día se produzcan más de 20 baterías defectuosas?


```{r}
1 - ppois(q = 20, lambda = 12)
```


**La probabilidad de que un día se produzcan más de 20 baterías defectuosas es de 1.159%**


##### 4.2 ¿Cuál es la probabilidad de que un día no salga ninguna batería defectuosa de la fábrica?


```{r}
dpois(x = 0, lambda = 12)
```


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


##### 4.3 La fábrica funciona de lunes a viernes. ¿Qué distribución sigue el número de baterías defectuosas por semana? Justifica qué propiedad se aplica.

*Una de las caracteristicas de una función de Poisson es la **aditividad**. Según el teorema de la aditividad, la variable suma de dos o más variables independientes que tengan una distribución de Poisson de distintos parámetros λ (de distintas medias) se distribuirá también con una distribución de Poisson con parámetro λ, equivalente a la suma de los parámetros λ (con media, la suma de las medias).*

$$ Y = \sum_{j=1}^{m} X{j}, X{j} \equiv Poiss(\lambda{j}) \to Y \equiv Poiss \left(\sum_{j=1}^{m}\lambda{j}\right)$$
*Por lo tanto, para 5 días (lunes - viernes) el número de baterias defectuosas sigue una distribución de Poisson con la media* $\lambda_{lunes} + \lambda_{martes} + \lambda_{miércoles} + \lambda_{jueves} + \lambda_{viernes} = 12 + 12 + 12 + 12 + 12 = 60$

**Si el valor de λ hubiera sido lo suficientemente grande, tendiendo a infinito, se podría considerar una aproximación utilizando la distribucción normal, aunque inicialmente se ha considerado como una poisson **

### 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.

##### 5.1 Realiza una simulación de la producción semanal de baterías (recuerda: 5 días de produccción, a 1000 baterías por día). Guarda los datos en un vector.

```{r}
sim = rweibull(n=(5*1000), shape = 100, scale = 185)

sim=data.frame(rweibull(1000*5, shape=100, scale=185))
colnames(sim) = c("values")

ggplot(sim, aes(values)) + 
  geom_histogram(aes(fill = ..count..),binwidth = 0.5, fill='grey') + 
  geom_density(aes(y=0.5 * ..count..), alpha = 0.3, fill = '#FF6666') + 
  ggtitle('Histograma datos simulados') + 
  theme(plot.title = element_text(hjust = 0.5))

```

##### 5.2 Con este nuevo proceso, ¿se mejora realmente la duración media de las baterías? (ayuda: puedes usar los datos simulados o la expresión de la esperanza de una Weibull)

**Duración media de las baterias usando los datos iniciales**
```{r}
# Usando los datos simulados
mean(datosB$lifetime)
```

**Duración media de las baterias usando la simulación**
```{r}
# Usando los datos simulados
mean(sim$values)

# Usando la expresión de la esperanza de una Weibull

rw = 185*gamma(1+(1/100))
mean(rw)

```

**Si comparamos la media de duración de las baterias tipo B con un valor de 179.6805 y las baterias provenientes de datos simulados con un valor de 183.9878 vemos que la segunda opción es mejor. **

```{r}
mean(sim$values)/mean(datosB$lifetime)-1
```
**Porcentaje de mejora: 2.37%**

##### 5.3 Los ingenieros no lo tienen muy claro (parece que la diferencia no es tanta en promedio y los nuevos materiales son costosos). Para demostrarles que merece la pena, calcula la proporción de baterías defectuosas que producirá probablemente el nuevo proceso y compárala con el anterior (la p que calculamos en la actividad 2)

*Calculamos la probabilidad p de que las baterias simuladas dure menos de 175 horas*

```{r}
dataBWeiblt175 =  pweibull(q=175, shape = 100, scale = 185 )
dataBWeiblt175
```

**Porcentaje de mejora:**
```{r}
dataBWeiblt175/dataBlt175-1
```

**La probabilidad de encontrar una bateria defectuosa usando los datos simulados es de 0.385% mientras que para la bateria tipo B, usando los datos iniciales era de 1.238%, por lo tanto podemos recomendar a la empresa el nuevo proceso de producción ya que hay una bajada de 68.898% en la proporción de baterias defectuosas.**