21_DDisc.Rmd

# Distribuciones discretas {#discretas}

En este capítulo se mostrarán las funciones de R para distribuciones discretas.

## Funciones disponibles para distribuciones discretas
Para cada distribución discreta se tienen 4 funciones, a continuación el listado de funciones y su utilidad.

```{r, eval=FALSE}
dxxx(x, ...)  # Función de masa de probabilidad, f(x)
pxxx(q, ...)  # Función de distribución acumulada hasta q, F(x)
qxxx(p, ...)  # Cuantil para el cual P(X <= q) = p
rxxx(n, ...)  # Generador de números aleatorios.
```

En el lugar de las letras `xxx` se de debe colocar el nombre de la distribución en R, a continuación el listado de nombres disponibles para las 6 distribuciones discretas básicas.

```{r, eval=FALSE}
binom     # Binomial
geo       # Geométrica
nbinom    # Binomial negativa
hyper     # Hipergeométrica
pois      # Poisson
multinom  # Multinomial
```

Combinando las funciones y los nombres se tiene un total de 24 funciones, por ejemplo, para obtener la función de masa de probabilidad $f(x)$ de una binomial se usa la función `dbinom( )` y para obtener la función acumulada $F(x)$ de una Poisson se usa la función `ppois( )`.

### Ejemplo binomial {-}
Suponga que un grupo de agentes de tránsito sale a una vía principal para revisar el estado de los buses de transporte intermunicipal. De datos históricos se sabe que un 10\% de los buses generan una mayor cantidad de humo de la permitida. En cada jornada los agentes revisan siempre 18 buses, asuma que el estado de un bus es independiente del estado de los otros buses.

1) Calcular la probabilidad de que se encuentren exactamente 2 buses que generan una mayor cantidad de humo de la permitida.

Aquí se tiene una distribucion $Binomial(n=18, p=0.1)$ y se desea calcular $P(X=2)$. Para obtener esta probabilidad se usa la siguiente instrucción.

```{r}
dbinom(x=2, size=18, prob=0.10)
```

Así $P(X=2)=0.2835$.

2) Calcular la probabilidad de que el número de buses que sobrepasan el límite de generación de gases sea al menos 4.

En este caso interesa calcular $P(X \geq 4)$, para obtener esta probabilidad se usa la siguiente instrucción.

```{r}
sum(dbinom(x=4:18, size=18, prob=0.10))
```

Así $P(X \geq 4)=0.0982$

3) Calcular la probabilidad de que tres o menos buses emitan gases por encima de lo permitido en la norma.

En este caso interesa $P(X\leq3)$ lo cual es $F(x=3)$, por lo tanto, la instrucción para obtener esta probabilidad es

```{r}
pbinom(q=3, size=18, prob=0.10)
```

Así $P(X\leq3)=F(x=3)=0.9018$

4) Dibujar la función de masa de probabilidad.

Para dibujar la función de masa de probabilidad para una $Binomial(n=18, p=0.1)$ se usa el siguiente código.

```{r binom1, fig.cap='Función de masa de probabilidad para una $Binomial(n=18, p=0.1)$.', fig.asp=0.7, fig.width=8}
x <- 0:18  # Soporte (dominio) de la variable
Probabilidad <- dbinom(x=x, size=18, prob=0.1)
plot(x=x, y=Probabilidad, 
     type='h', las=1, lwd=6)
```

En la Figura \@ref(fig:binom1) se muestra la función de masa de probabilidad para la $Binomial(n=18, p=0.1)$, de esta figura se observa claramente que la mayor parte de la probabilidad está concentrada para valores pequeños de $X$ debido a que la probabilidad de éxito individual es $p=0.10$. Valores de $X \geq 7$ tienen una probabilidad muy pequeña y es por eso que las longitudes de sus barras son muy cortas.

5) Generar con 100 de una distribución $Binomial(n=18, p=0.1)$ y luego calcular las frecuencias muestrales y compararlas con las probabilidades teóricas.

La muestra aleatoria se obtiene con la función `rbinom` y los resultados se almacenan en el objeto `m`, por último se construye la tabla de frecuencias relativas, a continuación el código usado.

```{r}
m <- rbinom(n=100, size=18, prob=0.1)
m  # Para ver lo que hay dentro de m
prop.table(table(m))  # Tabla de frecuencia relativa
```
A pesar de ser una muestra aleatoria de sólo 100 observaciones, se observa que las frecuencias relativas obtenidas son muy cercanas a las mostradas en la Figura \@ref(fig:binom1).

### Ejemplo geométrica {-}
En una línea de producción de bombillos se sabe que sólo el 1\% de los bombillos son defectuosos. Una máquina automática toma un bombillo y lo prueba, si el bombillo enciende, se siguen probando los bombillos hasta que se encuentre __un__ bombillo defectuoso, ahí se para la línea de producción y se toman los correctivos necesarios para mejorar el proceso.

1) Calcular la probabilidad de que se necesiten probar 125 bombillos para encontrar el primer bombillo defectuoso.

En la distribución geométrica, la variable $X$ representa el número de fracasos antes de encontrar el único éxito, por lo tanto, en este caso el interés es calcular $P(X=124)$. La instrucción para obtener esta probabiliad es la siguiente.

```{r}
dgeom(x=124, prob=0.01)
```

2) Calcular $P(X \leq 8)$.

En este caso interesa $P(X \leq 50)$ lo que equivale a $F(8)$, la instrucción para obtener la probabilidad es la siguiente.

```{r}
pgeom(q=50, prob=0.01)
```

3) Encontrar el cuantil $q$ tal que $P(X \leq q) = 0.40$.

En este caso interesa encontrar el cuantil $q$ que cumpla la condición de que hasta $q$ esté el 40\% de las observaciones, por esa razón se usa la función `qgeom` como se muestra a continuación.

```{r}
qgeom(p=0.4, prob=0.01)
```

```{block2, type='rmdnote'}
Note que las funciones `pxxx` y `qxxx` están relacionadas, `pxxx` entrega la probabilidad hasta el cuantil $q$ mientras `qxxx` entrega el cuantil en el que se acumula $p$ probabilidad.
```

### Ejemplo binomial negativa {-}
Una familia desea tener hijos hasta conseguir  __2 niñas__, la probabilidad individual de obtener una niña es 0.5 y se supone que todos los nacimientos son individuales, es decir, un sólo bebé.

1) Calcular la probabilidad de que se necesiten 4 hijos, es decir, 4 nacimientos para consguir las dos niñas.

En este problema se tiene una distribución binomial negativa con $r=2$ niñas, los éxitos deseados por la familia. La variable $X$ representa los fracasos, es decir los niños, hasta que se obtienen los éxitos $r=2$ deseados.

En este caso lo que interesa es $P(\text{familia tenga 4})$, en otras palabras interesa $P(X=2)$, la instrucción para calcular la probabilidad es la siguiente.

```{r}
dnbinom(x=2, size=2, prob=0.5)
```

2) Calcular $P(\text{familia tenga al menos 4 hijos})$.

Aquí interesa calcular $P(X \geq 2)=P(X=2)+P(X=3)+\ldots$, como esta probabilidad va hasta infinito, se debe usar el complemento así:

$$P(X \geq 2) = 1 - [P(X=0)+P(X=1)]$$
y para obtener la probabilidad solicitada se puede usar la función `dnbinom` de la siguiente manera.

```{r}
1 - sum(dnbinom(x=0:1, size=2, prob=0.5))
```

Otra forma para obtener la probabilidad solicitada es por medio de la función `pnbinom` de la siguiente manera.

```{r}
1 - pnbinom(q=1, size=2, prob=0.5)
```

### Ejemplo hipergeométrica {-}
Un lote de partes para ensamblar en una empresa está formado por 100 elementos del proveedor A y 200 elementos del proveedor B. Se selecciona una muestra de 4 partes al azar sin reemplazo de las 300 para una revisión de calidad.

1) Calcular la probabilidad de que todas las 4 partes de la muestra sean del proveedor A.

Aquí se tiene una situación que se puede modelar por medio de una distribución hipergeométrica con $m=100$ éxitos en la población, $n=200$ fracasos en la población y $k=4$ el tamaño de la muestra. El objetivo es calcular $P(X=4)$, para obtener esta probabilidad se usa la siguiente instrucción.

```{r}
dhyper(x=4, m=100, k=4, n=200)
```

2) Calcular la probabilidad de que dos o más de las partes sean del proveedor A.

Aquí interesa $P(X \geq 2)$, la instrucción para obtener esta probabilidad es.

```{r}
sum(dhyper(x=2:4, m=100, k=4, n=200))
```

### Ejemplo Poisson {-}
En una editorial se asume que todo libro de 250 páginas tiene en promedio 50 errores.

1) Encuentre la probabilidad de que en una página cualquiera no se encuentren errores.

Este es un problema de distribución Poisson con tasa promedio de éxitos dada por:

$$\lambda=\frac{50 \quad errores}{libro}=\frac{0.2 \quad errores}{pagina}$$
El objetivo es calcular $P(X=0)$, para obtener esta probabilidad de usa la siguiente instrucción.

```{r}
dpois(x=0, lambda=0.2)
```
Así $P(X=0)=0.8187$.

## Distribuciones discretas generales
En la práctica nos podemos encontramos con variables aleatorias discretas que no se ajustan a una de las distribuciones mostradas anteriormente, en esos casos, es posible manejar ese tipo de variables por medio de unas funciones básicas de R como se muestra en el siguiente ejemplo.

### Ejemplo {-}
El cangrejo de herradura hembra se caracteriza porque su caparazón se adhieren los machos de la misma especie, en la Figura \@ref(fig:crab) se muestra una fotografía de este cangrejo. Los investigadores están interesado en determinar cual es el patrón de variación del número de machos sobre cada hembra, para esto, se recolectó una muestra de hembras a las cuales se les observó el color, la condición de la espina, el peso en kilogramos, el ancho del caparazón en centímetros y el número de satélites o machos sobre el caparazón, la base de datos está disponible en el siguiente [enlace](https://raw.githubusercontent.com/fhernanb/datos/master/crab). \label{crabs}

```{r crab, echo=F, fig.cap='Fotografía del cangrejo de herradura, tomada de http://sccoastalresources.com/home/2016/6/21/a-night-of-horseshoe-crab-tagging', dpi=250, fig.align='center'}
knitr::include_graphics("images/crab.png")
```

1) Encontrar la distribución de probabilidad para la variable `Sa` que corresponde al número de machos sobre el caparazón de cada hembra.

Primero se debe leer la base de datos usando la url suministrada y luego se construye la tabla de frecuencia relativa y se almacena en el objeto `t1`.

```{r, echo=T}
url <- 'https://raw.githubusercontent.com/fhernanb/datos/master/crab'
crab <- read.table(file=url, header=T)

t1 <- prop.table(table(crab$Sa))
t1
```

La anterior tabla de frecuencias relativas se puede representar gráficamente usando el siguiente código.

```{r pmfcrab, fig.cap='Función de masa de probabilidad para el número de satélites por hembra.', fig.asp=0.7, fig.width=8}
plot(t1, las=1, lwd=5, xlab='Número de satélites',
     ylab='Proporción')
```

2) Sea $X$ la variable número de satélites por hembra, construir la función $F(x)$.

Para construir $F(x)$ se utiliza la función `ecdf` o _empirical cumulative density function_, a esta función le debe ingresar el vector con la información de la variable cuantitativa, a continuación del código usado. En la Figura \@ref(fig:Fcrab) se muestra la función de distribución acumulada para para el número de satélites por hembra.

```{r Fcrab, fig.cap='Función de distribución acumulada para el número de satélites por hembra.', fig.asp=0.7, fig.width=8}
F <- ecdf(crab$Sa)
plot(F, las=1, main='')
```

3) Calcular $P(X \leq 9)$.

Para obtener esta probabilidad se usa el objeto `F` que es en realidad una función, a continuación la instrucción usada.

```{r}
F(9)
```

Así $P(X \leq 9)=0.9595$.

4) Calcular $P(X > 4)$.

Para obtener esta probabilidad se usa el hecho de que $P(X > 4) = 1 - P(X \leq 4)$, así la instrucción a usar es.

```{r}
1 - F(4)
```

Por lo tanto $P(X > 4)=0.2775$.

5) Suponga que el grupo 1 está formado por las hembras cuyo ancho de caparazón es menor o igual al ancho mediano, el grupo 2 está formado por las demás hembras. ¿Será $F(x)$ diferente para los dos grupos?

Para realizar esto vamos a particionar el vector `Sa` en los dos grupos de acuerdo a la nueva variable `grupo` creada como se muestra a continuacion.

```{r}
grupo <- ifelse(crab$Wt <= median(crab$Wt), 'Grupo 1', 'Grupo 2')
x <- split(x=crab$Sa, f=grupo)
```

El objeto `x` es una lista y para acceder a los vectores allí almacenados usamos dos corchetes `[[]]`, uno dentro del otro. Luego para calcular $F(x)$ para los dos grupos se procede así:

```{r}
F1 <- ecdf(x[[1]])
F2 <- ecdf(x[[2]])
```

Para obtener las dos $F(x)$ en la misma figura se usa el código siguiente.

```{r Fcrabs, fig.cap='Función de distribución acumulada para el número de satélites por hembra diferenciando por grupo.', fig.asp=0.7, fig.width=8}
plot(F1, col='blue', main='', las=1)
plot(F2, col='red', add=T)
legend('bottomright', legend=c('Grupo 1', 'Grupo 2'),
       col=c('blue', 'red'), lwd=1)
```

En la Figura \@ref(fig:Fcrabs) se muestran las dos $F(x)$, en color azul para el grupo 1 y en color rojo para el grupo 2. Se observa claramente que las curvas son diferentes antes de $x=9$. El hecho de que la curva azul esté por encima de la roja para valores menores de 9, es decir, $F_1(x) \geq F_2(x)$, indica que las hembras del grupo 1 tienden a tener menos satélites que las del grupo 2, esto es coherente ya que las del grupo 2 son más grandes en su caparazón.