Distribuciones de probabilidad discretas

Introducción a las distribuciones de probabilidad

Se puede describir de forma completa una variable aleatoria sin más que especificar la probabilidad asociada a cada uno de sus posibles valores. Esta especificación se conoce con el nombre de distribución de probabilidad. Sin embargo, la forma en que se puede especificar dicha distribución de probabilidad depende del tipo de variable aleatoria con la que se esté trabajando. En el caso de variables discretas basta con determinar la probabilidad de cada uno de los posibles resultados observables de la variable, pero no ocurre de la misma forma con las variables continuas donde es imposible conocer todas las probabilidades asociadas con cada uno de sus posibles valores.

Distribuciones de probabilidad discretas

Se les denomina discretas, debido a que sus eventos no asumen valores intermedios entre dos numeros enteros.

Simulación de distribuciones de probabilidad

Método de la transformación inversa

El método de la transformada inversa se utiliza para simular variables aleatorias continuas, mediante la función de distribución acumulada y la generación de números aleatorios con distribución uniforme en el intervalo (0,1).

Método de convolución

En algunas distribuciones de probabilidad la variable aleatoria a simular puede generarse mediante la suma de otras variables aleatorias, de manera tal que es más rápida la simulación que utilizando otros métodos; y puede expresarse como:

\[{Y} = {X}_{1} + {X}_{1} + \cdots + {X}_{n}\]

Método de transformación directa.

Este método se utiliza para generar variables aleatorias con distribución normal.

Metodo de composición.

A partir de este método se pueden generar variables aleatorias cuando la función de densidad se expresa como una combinación de \(m\) distribuciones, es decir, se encuentra definida por partes; aquí la función de densidad se puede expresar como:

\[ {f}\left({x}\right) = {\sum}_{i=1}^{n}{f}_{i}\left({x}\right){I}_{A}\left({x}\right) \]

Con

\[ {I}_{A}\left({x}\right) = \begin{cases} 1 & \text{if } x \in A\\ 0 & \text{if } x \notin A \end{cases} \]

Bootstrap paramétrico

El método bootstrap puede ser usado tanto para el cálculo de errores estándar como de intervalos de confianza en un modelo paramétrico.

• En bootstrap no paramétrico son obtenidas muestras \({X}_{1}^{*},{X}_{2}^{*},\ldots,{X}_{n}^{*}\) de la distribución empírica \({P}_{n}\).
• En bootstrap paramétrico las muestras se obtienen de \({P}(X,\widehat{\theta})\) donde \(\widehat{\theta}\) es una estimación de \({\theta}\).

La diferencia entre la versión no paramétrica y la paramétrica radida en cómo se contruye la distribución de la que van a ser seleccionadas las muestras.

Introducción a la simulación de distribuciones de probabilidad discretas

Los métodos de simulación de variables aleatorias continuas pueden adaptarse para variables aleatorias continuas, como en el caso particular del método de la transformada inversa.

En general, para utilizar este método para simular una variable aleatoria discreta se sigue el siguiente procedimiento:

• Calcular la probabilidad para la distribución de probabilidad, \(f(x)\), de la variable a modelar.

• Calcular las probabilidades acumuladas, \(F(x)\).

• Generar números psudoaleatorios de la distribución uniforme en \((0,1)\).

• Comparar el valor de \(F(x)\) y determinar los valores de \(x\) que correspondan a \(F(x)\).

Sea \(X\) una variable aleatoria discreta para la cual \({P}({X}={x}_{j})={p}_{j}\), y sea \({U}\) una variable aleatoria continua con distribución uniforme en el intervalo \((0,1)\) \({U}{\sim}U(0,1)\) y además:

\[ {X} = \begin{cases} {x}_{1}\text{ si }{U}{\leq}{p}_{1}\\ {x}_{2}\text{ si }{p}_{1}{<}{U}{\leq}{p}_{1}+{p}_{2}\\ {x}_{3}\text{ si }{p}_{1}+{p}_{2}{<}{U}{\leq}{p}_{1}+{p}_{2}+{p}_{3}\\ {\vdots}\\ {x}_{j}\text{ si }{p}_{1}+{p}_{2}+\cdots+{p}_{j-1}{<}{U}{\leq}{p}_{1}+{p}_{2}+\cdots+{p}_{j} \end{cases} \] con \({\sum}_{i=1}^{j}{p}_{j}=1\)

Simulación de una distribución empírica.

Ejemplos

Durante 50 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  3,  3/50,
  "1",  8,  8/50,
  "2", 10, 10/50,
  "3", 15, 15/50,
  "4",  9,  9/50,
  "5",  3,  3/50,
  "6",  2,  2/50
)

La simulación de esta distribución de probabilidad empírica se muestra a continuación.

library(car)
U <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.06=0;
       0.06:0.22=1;
       0.22:0.42=2;
       0.42:0.72=3;
       0.72:0.90=4;
       0.90:0.96=5;
       0.96:1=6"
X <- recode(U,CL)
tabla.absolutas <- table(x = X)
simula.empirica <- as.data.frame(x = tabla.absolutas / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.empirica, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.empirica <- as.data.frame(cbind(x = c(0,1,2,3,4,5,6), Freq = c(0.06,0.16, 0.2,0.3,0.18,0.06,0.04)))

library(ggplot2)
p2 <- ggplot(data = variable.empirica, aes(x = as.character(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

grid.arrange(p1, p2, ncol = 2, top = "Distribución empirica")

Ejercicios

• Durante 96 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  6,  6/96,
  "1",  16,  16/96,
  "2", 20, 20/96,
  "3", 30, 30/96,
  "4",  18,  18/96,
  "5",  6,  6/96,
)

• Durante 100 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  6,  6/100,
  "1",  16,  16/100,
  "2", 20, 20/100,
  "3", 30, 30/100,
  "4",  18,  18/100,
  "5",  6,  6/100,
  "6",  4,  4/100
)

• Durante 115 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  16,  16/116,
  "1",  17,  17/116,
  "2", 20, 20/116,
  "3", 35, 35/116,
  "4",  22,  22/116,
  "5",  6,  6/116,
)

• Durante 180 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  11,  11/180,
  "1",  33,  33/180,
  "2", 35, 35/180,
  "3", 23, 23/180,
  "4",  30,  30/180,
  "5",  25,  25/180,
  "6",  23,  23/180
)

• Durante 200 días se observó el número de unidades de un artículo que se demandaron en una tienda, obteniéndose los datos que se muestran a continuación:

library(tidyverse)
tribble(
  ~demanda, ~días, ~frecuencia,
  #--|--|----
  "0",  12,  12/200,
  "1",  35,  35/200,
  "2", 38, 38/200,
  "3", 27, 27/200,
  "4",  31,  31/200,
  "5",  29,  29/200,
  "6",  28,  28/200
)

Distribución de probabilidad uniforme \(X{\sim}Unif(min, max)\)

La función de probabilidad de una variable aleatoria discreta tiene la forma:

\[ P\left(X=x\right)=\frac{1}{max-min}\text{ con }min{\leq}{x}{\leq}max \]

Simulación de una distribución uniforme.

\[ P\left(X=x\right)=\frac{1}{10-0}\text{ con }0{\leq}{x}{\leq}10 \]

y entonces

\[ {X} = \begin{cases} {1}\text{ si }0{<}{U}{\leq}0.1\\ {2}\text{ si }0.1{<}{U}{\leq}0.2\\ {3}\text{ si }0.2{<}{U}{\leq}0.3\\ {4}\text{ si }0.3{<}{U}{\leq}0.4\\ {5}\text{ si }0.4{<}{U}{\leq}0.5\\ {6}\text{ si }0.5{<}{U}{\leq}0.6\\ {7}\text{ si }0.6{<}{U}{\leq}0.7\\ {8}\text{ si }0.7{<}{U}{\leq}0.8\\ {9}\text{ si }0.8{<}{U}{\leq}0.9\\ {10}\text{ si }0.9{<}{U}{\leq}1.0\\ \end{cases} \]

library(car)
U <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.1=1;
       0.1:0.2=2;
       0.2:0.3=3;
       0.3:0.4=4;
       0.4:0.5=5;
       0.5:0.6=6;
       0.6:0.7=7;
       0.7:0.8=8;
       0.8:0.9=9;
       0.9:1=10"
X <- recode(U,CL)
tabla.absolutas <- table(x = X)
simula.uniforme <- as.data.frame(x = tabla.absolutas / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.uniforme, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.uniforme <- as.data.frame(cbind(x = 1:10, Freq = c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)))

library(ggplot2)
p2 <- ggplot(data = variable.uniforme, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución uniforme")

Ejercicios

• Simule una distribución uniforme con la siguiente densidad

\[ P\left(X=x\right)=\frac{1}{2-0}\text{ con }0{\leq}{x}{\leq}2 \]

• Simule una distribución uniforme con la siguiente densidad

\[ P\left(X=x\right)=\frac{1}{7-0}\text{ con }0{\leq}{x}{\leq}7 \]

• Simule una distribución uniforme con la siguiente densidad

\[ P\left(X=x\right)=\frac{1}{13-0}\text{ con }0{\leq}{x}{\leq}13 \]

• Simule una distribución uniforme con la siguiente densidad

\[ P\left(X=x\right)=\frac{1}{18-0}\text{ con }0{\leq}{x}{\leq}18 \]

• Simule una distribución uniforme con la siguiente densidad

\[ P\left(X=x\right)=\frac{1}{26-0}\text{ con }0{\leq}{x}{\leq}26 \]

Distribución de probabilidad binomial \(X{\sim}Bin(size, prob)\)

Se le conoce tambien como distribucion de la probabilidad puntual o de Bernoulli, debido al suizo Jacques Bernoulli, quien por primera vez desarrollo el concepto de ensayos independientes.

Función de probabilidades \(f_{X}\)

\[ P\left(X=x\right)=\binom{size}{x}{prob}^{x}\left(1-{prob}\right)^{size-x}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,2,\ldots,{size} \]

densidad.binomial <- function(x, size, prob) choose(size,x)*prob**x*(1-prob)**(size-x)

Esperanza \(\mu_{x}\)

\[ \mu_{x}=\sum_{x=0}^{size}x{\cdot}\binom{size}{x}{prob}^{x}\left(1-{prob}\right)^{size-x}={size}{\cdot}{prob} \]

esperanza.binomial <- function(size, prob) sum((0:size)*choose(size,(0:size))*prob**(0:size)*(1-prob)**(size-(0:size)))

Varianza \(\sigma_{x}^{2}\)

\[ \sigma_{x}^{2}=\sum_{x=0}^{size}\left(x-\mu_{x}\right)^{2}{\cdot}\binom{size}{x}{prob}^{x}\left(1-{prob}\right)^{{size}-x}={size}{\cdot}{prob}{\cdot}\left(1-{prob}\right) \]

varianza.binomial <- function(size, prob) sum((0:size-esperanza.binomial(size, prob))**2*choose(size,(0:size))*prob**(0:size)*(1-prob)**(size-(0:size)))

Características:

1. Los ensayos de Bernoulli se llevan a cabo size (o n) veces,

2. Los ensayos son independientes,

3. La probabilidad de exito prob (o p) no cambia entre los ensayos,

4. Solo hay dos posibles resultados en cada ensayo.

Ejemplo

Una medida de aislamiento preventivo obligatorio para contener la expansión de un virus mortal, el COVID-19, ha exito hasta el punto de que 80% de los habitantes de un país se han quedado en sus casas. Para verificar el acatamiento de la medida un grupo de 400 familias ha sido seleccionada para hacer un seguimiento a dicha medida. Con base en lo anterior encuentre las siguientes probabilidades.

\[ P(\text{# de familias que acatan la norma}=x)=\binom{400}{x}{0.8}^{x}\left(1-{0.8}\right)^{400-x}\text{ con }0{\leq}{0.8}{\leq}1\text{ y }x=0,1,2,\ldots,{400} \]

1. ¿Cuál es la probabilidad de que en el grupo hayan cumplido estrictamente con la norma 200 familias?

\[ P(\text{# de familias que acatan la norma}=200)=\binom{400}{200}{0.8}^{200}\left(1-{0.8}\right)^{400-200} \]

choose(400,200)*0.8**200*(1-0.8)**(400-200)

## [1] 6.864888e-41

densidad.binomial(x = 200, size = 400, prob = 0.8)

## [1] 6.864888e-41

dbinom(x = 200, size = 400, prob = 0.8)

## [1] 6.864888e-41

2. ¿Cuál es la probabilidad de que en el grupo hayan cumplido estrictamente con la norma 320 familias?

\[ P(\text{# de familias que acatan la norma}=320)=\binom{400}{320}{0.8}^{320}\left(1-{0.8}\right)^{400-320} \]

choose(400,320)*0.8**320*(1-0.8)**(400-320)

## [1] 0.04981327

densidad.binomial(x = 320, size = 400, prob = 0.8)

## [1] 0.04981327

dbinom(x = 320, size = 400, prob=0.80)

## [1] 0.04981327

3. ¿Y cómo máximo 200 cumplan con la norma de aislamiento preventivo obligatorio?

\[ P(\text{# de familias que acatan la norma }{\leq}200)=\sum_{x=0}^{200}P(\text{# de familias que acatan la norma}=x) \] \[ \sum_{x=0}^{200}P(\text{# de familias que acatan la norma}=x)=\sum_{x=0}^{200}\binom{400}{x}{0.8}^{x}\left(1-{0.8}\right)^{400-x} \]

sum(densidad.binomial(0:200, 400, 0.8))

## [1] 9.128211e-41

sum(dbinom(x = 0:200, size = 400, prob = 0.80))

## [1] 9.128211e-41

pbinom(q = 200, size = 400, prob = 0.80, lower.tail = TRUE)

## [1] 9.128211e-41

4. ¿Y como mínimo 300 familias sean las que acatan la norma de aislamiento preventivo obligatorio?

\[ P(\text{# de familias que acatan la norma }{\geq}300)=1-P(\text{# de familias que acatan la norma}{\leq}299) \]

\[ 1-P(\text{# de familias que acatan la norma}{\leq}299)=1-\sum_{x=0}^{299}P(\text{# de familias que acatan la norma}=x) \]

\[ 1-\sum_{x=0}^{299}P(\text{# de familias que acatan la norma}=x)=1-\sum_{x=0}^{299}\binom{400}{x}{0.8}^{x}\left(1-{0.8}\right)^{400-x} \]

1 - sum(densidad.binomial(0:299, 400, 0.8))

## [1] 0.9938071

1 - sum(dbinom(x = 0:299, size = 400, prob = 0.80))

## [1] 0.9938071

1 - pbinom(q = 299, size = 400, prob = 0.80, lower.tail = TRUE)

## [1] 0.9938071

También se puede encontrar esta probabilidad haciendo el calculo siguiente:

\[ P(\text{# de familias que acatan la norma }{\geq}300)=\sum_{x=300}^{400}P(\text{# de familias que acatan la norma}=x) \]

\[ \sum_{x=300}^{400}P(\text{# de familias que acatan la norma}=x)=\sum_{x=300}^{400}\binom{400}{x}{0.8}^{x}\left(1-{0.8}\right)^{400-x} \]

sum(densidad.binomial(300:400, 400, 0.8))

## [1] 0.9938071

sum(dbinom(x = 300:400, size = 400, prob = 0.80))

## [1] 0.9938071

pbinom(q = 299, size = 400, prob = 0.80, lower.tail = FALSE)

## [1] 0.9938071

Nota

La suma de las probabilidades sobre todos los valores que puede tomar la variable aleatoria es igual a 1

\[ P(\text{# de familias que acatan la norma }{\leq}299)+P(\text{# de familias que acatan la norma }{\geq}300)=1 \]

pbinom(q = 299, size = 400, prob = 0.80, lower.tail = TRUE) + pbinom(q = 299, size = 400, prob = 0.80, lower.tail = FALSE)

## [1] 1

De lo anterior y despejando, se puede obtener la formula del complemento de una probabilidad, y es por esto que las probabilidades son iguales

all.equal(pbinom(q = 299, size = 400, prob = 0.80, lower.tail = FALSE), 1 - pbinom(q = 299, size = 400, prob = 0.80, lower.tail = TRUE))

## [1] TRUE

5. ¿Y qué sean como máximo las 400 las familias las que acatan la norma, es decir, cero familias o una familia o dos familias, etc?

\[ P(\text{# de familias que acatan la norma} {\leq} 400)=\sum_{x = 0}^{400}{P(\text{# de familias que acatan la norma}=x)} \]

\[ \sum_{x = 0}^{400}{P(\text{# de familias que acatan la norma}=x)}=\sum_{x=0}^{400}\binom{400}{x}{0.8}^{x}\left(1-{0.8}\right)^{400-x} \]

sum(densidad.binomial(0:400, 400, 0.8))

## [1] 1

sum(dbinom(x = 0:400, size = 400, prob = 0.80))

## [1] 1

pbinom(q = 400, size = 400, prob = 0.80, lower.tail = TRUE)

## [1] 1

6. Analizar el éxito o fracaso de la medida de asilamento preventivo en la muestra aleatoria de 400 famiilias.

ensayos <- 400
pobabilidad.de.exito <- 0.80
valores.posibles <- 0:ensayos
fx.binomial <- dbinom(x = valores.posibles, size = ensayos, prob = pobabilidad.de.exito)
fbb = pbinom(q = valores.posibles, size = ensayos, prob = pobabilidad.de.exito, lower.tail = TRUE)
fx.binomial = data.frame(x = valores.posibles, "f(X)" = fx.binomial, "F(X)" = fbb)
#tail(fx.binomial)

library(ggplot2)
p1 <- ggplot(fx.binomial, aes(x = x, y = f.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(ggplot2)
p2 <- ggplot(fx.binomial, aes(x = x, y = F.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Distribución F(x)")

library(gridExtra)
grid.arrange(p1, p2, nrow = 2, top = "Distribución binomial")

7. Calcular el valor esperado del número de familias que se espera acaten la norma al seleccionar una muestra aleatoria de 400

esperanza.binomial(400,0.8)

## [1] 320

400*0.8

## [1] 320

8. Calcular la varianza del número de familias que se espera acaten la norma al seleccionar una muestra aleatoria de 400

varianza.binomial(400,0.8)

## [1] 64

400*0.8*(1-0.8)

## [1] 64

Simulación de la distribución binomial

Una variable aleatoria discreta con distribución binomial con parámetros \(size\) y \(prob\) puede ser generada a partir de la suma de \(size\) variables aleatorias \(Bernoulli\) con parámetro \(prob\) (\({x}_{i}{\sim}Ber(prob)\)).

\[ {X}_{i}{\sim}Ber(1, prob){\implies}{X}=\sum_{x=1}^{size}{X}_{i}{\sim}Bin(size, prob) \]

Ejemplos

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=2\) y \(prob=0.7\)

library(car)
U1 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.3=0;0.3:1.0=1"
X1 <- car::recode(var = U1, recodes = CL)
U2 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.3=0;0.3:1.0=1"
X2 <- car::recode(var = U2, recodes = CL)
X <- X1 + X2
tabla.frecuencias <- table(x = X)
simula.binomial <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.binomial, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.binomial <- as.data.frame(cbind(x = 0:2, Freq = dbinom(x = 0:2, size = 2, prob = 0.7)))

library(ggplot2)
p2 <- ggplot(data = variable.binomial, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución bernoulli")

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=8\) y \(prob=0.5\)

library(car)
U1 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X1 <- car::recode(var = U1, recodes = CL)
U2 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X2 <- car::recode(var = U2, recodes = CL)
U3 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X3 <- car::recode(var = U3, recodes = CL)
U4 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X4 <- car::recode(var = U4, recodes = CL)
U5 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X5 <- car::recode(var = U5, recodes = CL)
U6 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X6 <- car::recode(var = U6, recodes = CL)
U7 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X7 <- car::recode(var = U7, recodes = CL)
U8 <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.5=0;0.5:1.0=1"
X8 <- car::recode(var = U8, recodes = CL)
X <- X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8
tabla.frecuencias <- table(x = X)
simula.binomial <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.binomial, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.binomial <- as.data.frame(cbind(x = 0:8, Freq = dbinom(x = 0:8, size = 8, prob = 0.5)))

library(ggplot2)
p2 <- ggplot(data = variable.binomial, aes(x = as.character(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución binomial")

Ejercicios

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=2\) y \(prob=0.5\)

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=2\) y \(prob=0.7\)

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=7\) y \(prob=1\)

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=12\) y \(prob=0.3\)

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=20\) y \(prob=0.3\)

Distribución de probabilidad geométrica \(X{\sim}Geom(1,prob)\)

Es cualquiera de las dos distribuciones de probabilidad discretas siguientes:

1. la distribución de probabilidad del número \(X\) del ensayo de Bernoulli necesaria para obtener un éxito, contenido en el conjunto \(\left\{ 1, 2, 3,\ldots,\right\}\) o
2. la distribución de probabilidad del número \(Y = X − 1\) de fallos antes del primer éxito, contenido en el conjunto \(\left\{ 0, 1, 2,\ldots\right\}\).

Función de probabilidades \(f_{X}\)

1. Aquí la idea es contar el número de ensayos antes de obtener el primer éxito con probabilidad prob

\[ P(X=x)={prob}\left(1-{prob}\right)^{x-1}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=1,2,\ldots \]

2. Aquí cuento el número de fracasos antes de obtener el primer éxito con probabilidad prob

\[ P(X=x)={prob}\left(1-{prob}\right)^{x}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,\ldots \]

densidad.geometrica <- function(x, prob) prob*(1-prob)**(x)

Esperanza \(\mu_{x}\)

\[ \mu_{x}=\sum_{x=0}^{\infty}x{\cdot}{prob}\left(1-{prob}\right)^{x}=\frac{1-prob}{prob} \]

esperanza.geometrica <- function(size, prob) sum((0:size)*prob*(1-prob)**(0:size))

Varianza \(\sigma_{x}^{2}\)

\[ \sigma_{x}^{2}=\sum_{x=0}^{\infty}\left(x-\mu_{x}\right)^{2}{\cdot}{prob}\left(1-{prob}\right)^{x}=\frac{1-{prob}}{{prob}^{2}} \]

varianza.geometrica <- function(size, prob) sum((0:size-esperanza.geometrica(size, prob))**2*prob*(1-prob)**(0:size))

Ejemplo

Si con relación a cierta pandemia que ha azotado la población de un país se sabe por una muestra aleatoria que un viajero ha contraido la enfermedad y que ha tenido contacto con cien personas, y dado que esta enfermedad es muy contagiosa teniendo un contacto cercano se transmite en el 99% de los casos.

\[ P(X=x)={0.99}\left(1-{0.99}\right)^{x}\text{ con }0{\leq}{0.99}{\leq}1\text{ y }x=0,1,\ldots \]

1. calcue la probabilidad que el viajero haya contagiado a la primera persona que contacto, es decir, no haya tenido que contactar a nadie previamente para transferir el mal al primer ser humano con quien estableció contacto.

\[ P(\text{# contactos previos al primer contagio}=0)={0.99}\left(1-0.99\right)^{0} \]

densidad.geometrica(0, 0.99)

## [1] 0.99

dgeom(x = 0, prob = 0.99)

## [1] 0.99

2. calcue la probabilidad de que tan solo la segunda persona en tener contacto con este viajero haya sido contagiada del virus, o en otras palabras, se haya dado tan solo un contacto previo del viajero antes de contaminar a alguien.

\[ P(\text{# contactos previos al primer contagio}=1)={0.99}\left(1-0.99\right)^{1} \]

densidad.geometrica(x = 1, prob = 0.99)

## [1] 0.0099

dgeom(x = 1, prob = 0.99)

## [1] 0.0099

3. calcue la probabilidad de que tan solo la tercera persona en tener contacto con este viajero haya sido contagiada del virus, o en otras palabras, el viajero haya tenido dos contactos previos antes de contaminar a un ser humano.

\[ P(\text{# contactos previos al primer contagio}=2)={0.99}\left(1-0.99\right)^{2} \]

densidad.geometrica(x = 2, prob = 0.99)

## [1] 9.9e-05

dgeom(x = 2, prob = 0.99)

## [1] 9.9e-05

4. halle la probabilidad de que el viajero no haya contagiado a los 99 primeros seres humanos en tener contacto con él pero si haya contagiado al último, es decir, se hayan dado 99 contactos previos a la contaminación del primer individuo.

\[ P(\text{# contactos previos al primer contagio}=99)={0.99}\left(1-0.99\right)^{99} \]

densidad.geometrica(x = 99, prob = 0.99)

## [1] 9.9e-199

dgeom(x = 99, prob = 0.99)

## [1] 9.9e-199

5. halle la probabilidad de que el viajero no haya contagiado ninguna de los primeros 100 seres humanos con los que tuvo contacto, es decir, se hayan dado 100 contactos previos y el primer contagiado sea el siguiente individuo en ser contactado.

\[ P(\text{# contactos previos al primer contagio}=99)={0.99}\left(1-0.99\right)^{100} \]

densidad.geometrica(x = 100, prob = 0.99)

## [1] 9.9e-201

dgeom(x = 100, prob = 0.99)

## [1] 9.9e-201

6. Si usted fuera el Ministro de salud de dicho país, estaría interesado en saber la probabilidad de que un viajero en partícular no contagie a otros con quien tuvo contacto; estudie dicha situación.

\[ P(\text{# contactos previos al primer contagio}>x)=1-P(\text{# contactos previos al primer contagio}{\leq}x) \]

\[ 1-P(\text{# contactos previos al primer contagio}{\leq}x)=1-\sum_{x=0}^{x}P(\text{# contactos previos al primer contagio}=x) \]

\[ 1-\sum_{x=0}^{x}P(\text{# contactos previos al primer contagio}=x)=1-\sum_{x=0}^{x}{0.99}\left(1-{0.99}\right)^{100-x} \]

ensayos <- 100
probabilidad.de.contagio <- 0.99
valores.posibles <- 0:ensayos
fx.geometrica <- dgeom(x = valores.posibles, prob = probabilidad.de.contagio)
fbg = pgeom(q = valores.posibles, prob = probabilidad.de.contagio, lower.tail = TRUE)
fx.geometrica = data.frame(x = valores.posibles, "f(X)" = fx.geometrica, "F(X)" = fbg)
#head(fx.geometrica)

library(ggplot2)
p3 <- ggplot(fx.geometrica, aes(x = x, y = f.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(ggplot2)
p4 <- ggplot(fx.geometrica, aes(x = x, y = F.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Distribución F(x)")

library(gridExtra)
grid.arrange(p3, p4, nrow = 2, top = "Distribución geométrica")

7. El número esperado de contactos previos al primer contagio es igual a:

esperanza.geometrica(size = 1000000, prob = 0.99)

## [1] 0.01010101

(1-0.99)/0.99

## [1] 0.01010101

8. La distancia promedio, elevada al cuadrado, al valor esperado del número de contactos previos al primer contagiado es igual a:

varianza.geometrica(size = 10000000, prob = 0.99)

## [1] 0.01020304

(1-0.99)/0.99**2

## [1] 0.01020304

Simulación de la distribución geométrica

La distribución geométrica tiene función de probabilidad dada por: \(P(X=x)={prob}\left(1-{prob}\right)^{x}\) con \(0{\leq}{prob}{\leq}1\) y \(x=0,1,\ldots\). Se realizan ensayos hasta que ocurre el primer éxito. Sea \(X\) la variable aleatoria discreta que indica el número de ensayos hasta que ocurre el primer éxito \(X:=\text{"# fallos hasta que ocurre el primer éxito"}\), y entonces:

\[ \begin{align} \sum_{i=1}^{size-1}P\left({X}={i}\right)&={1}-P\left({X}\geq{i}\right)\\ &={1}-\left({1}-{prob}\right)^{{size}-{1}} \end{align} \]

Luego, es posible simular una variable uniforme, \(U\), en el intervalo \(\left(0,1\right)\) y hacer \(X\) igual al \(j\)-ésimo ensayo para el cual se ha obtenido el primer éxtido, o en forma equivalente se puede definir una variable aleatoria \(1−U\) que tiene la misma distribución que \(U\), con lo que ahora es posible definir \(X\) de la siguiente manera:

\[ \begin{align} {X} &= \min_{j}{\left({1}-{prob}\right)^{{size}} < U}\\ &= \min_{j}{\left[\ln{\left({1}-{prob}\right)^{{size}}}\right] < \ln{U}}\\ &= \min_{j}{{\left[{size}\ln{\left({1}-{prob}\right)}\right]} < \ln{U}}\\ &= \min_{j}{{\left[{size}\right]} \geq \frac{\ln{U}}{\ln{\left({1}-{prob}\right)}}} \end{align} \]

Con lo que finalmente

\[ \begin{align} {X} &\geq {\Bigg\lfloor}{\frac{\ln{U}}{\ln{\left({1}-{prob}\right)}}}{\Bigg\rfloor} \end{align} \]

Ejemplos

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetros \(prob=0.8\)

prob = 0.8
U <- runif(n = 100000, min = 0, max = 1)
X <- floor(log(x = U, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
tabla.frecuencias <- table(x = X)
simula.geometrica <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.geometrica, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.geometrica <- as.data.frame(cbind(x = 0:max(x = X), Freq = dgeom(x = 0:max(x = X), prob = prob)))

library(ggplot2)
p2 <- ggplot(data = variable.geometrica, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución geométrica")

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetros \(prob=0.5\)

prob = 0.5
U <- runif(n = 100000, min = 0, max = 1)
X <- floor(log(x = U, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
tabla.frecuencias <- table(x = X)
simula.geometrica <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.geometrica, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.geometrica <- as.data.frame(cbind(x = 0:max(x = X), Freq = dgeom(x = 0:max(x = X), prob = prob)))

library(ggplot2)
p2 <- ggplot(data = variable.geometrica, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución geométrica")

Ejercicios

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetro \(prob=0.2\)

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetro \(prob=0.4\)

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetro \(prob=0.8\)

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetro \(prob=0.6\)

• Simular una variable aleatoria discreta que tiene distribución geompetrica con parámetro \(prob=0.3\)

Distribución de probabilidad binomial negativa \(X{\sim}BN(size,pob)\)

Es cualquiera de las dos distribuciones de probabilidad discretas siguientes:

1. la distribución de probabilidad del número \(X\) del ensayo de Bernoulli necesaria para obtener r éxitos, contenido en el conjunto \(\left\{ 1, 2, 3,\ldots,\right\}\) o
2. la distribución de probabilidad del número \(Y = X − 1\) de fallos antes de que se obtengan r éxitos, contenido en el conjunto \(\left\{ 0, 1, 2,\ldots\right\}\).

Función de probabilidades \(f_{X}\)

1. Aquí la idea es contar el número de ensayos para obtener size éxitos con probabilidad prob

\[ P(X=x)=\binom{x-1}{size-1}{prob}^{size}\left(1-{prob}\right)^{x-size}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,\ldots \]

o

\[ P(X=x)=\binom{x-1}{x-size}{prob}^{size}\left(1-{prob}\right)^{x-size}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,\ldots \]

2. Aquí cuento el número de fracasos para obtener el size-ésimo éxito con probabilidad prob

\[ P(X=x)=\binom{x+size-1}{size-1}{prob}^{size}\left(1-{prob}\right)^{x}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,\ldots \]

o

\[ P(X=x)=\binom{x+size-1}{x}{prob}^{size}\left(1-{prob}\right)^{x}\text{ con }0{\leq}{prob}{\leq}1\text{ y }x=0,1,\ldots \]

densidad.binomialnegativa <- function(x, size, prob) choose(x+size-1,x)*prob**size*(1-prob)**x

Esperanza \(\mu_{x}\)

\[ \mu_{x}=\sum_{x=0}^{\infty}x{\cdot}\binom{x+size-1}{x}{prob}^{size}\left(1-{prob}\right)^{x}=size\cdot\frac{1-prob}{prob} \]

esperanza.binomialnegativa <- function(size, prob) sum((0:size)*choose((0:size)+size-1,(0:size))*prob**size*(1-prob)**(0:size))

Varianza \(\sigma_{x}^{2}\)

\[ \sigma_{x}^{2}=\sum_{x=0}^{\infty}\left(x-\mu_{x}\right)^{2}{\cdot}{prob}\left(1-{prob}\right)^{x}=size\cdot\frac{1-{prob}}{{prob}^{2}} \]

varianza.binomialnegativa <- function(size, prob) sum((0:size-esperanza.binomialnegativa(size, prob))**2*choose((0:size)+size-1,(0:size))*prob**size*(1-prob)**(0:size))

Ejemplo

Si con relación a cierta pandemia que ha azotado la población de un país, se sabe por una muestra aleatoria que un viajero ha contraido la enfermedad y que ha contactado a cien personas. Dado que esta enfermedad es muy contagiosa, pues, teniendo un contacto cercano se transmite en el 99% de los casos. Si se sabe que este viajero ha contagiado a 25 personas.

\[ P(X=x)=\binom{x+25-1}{x}{0.99}^{25}\left(1-{0.99}\right)^{x}\text{ con }0{\leq}{0.99}{\leq}1\text{ y }x=0,1,\ldots \]

1. calcue la probabilidad que el viajero al contactar al vigésimo quinto ser humano haya contagiado a veintícinco personas, es decir, haya contactado y contagiado a cero personas al transferir el mal a veinticínco seres humanos con quien estableció contacto.

\[ P(\text{# no contagiados hasta que se dan veintinco contagios}=0)=\binom{0+25-1}{0}{0.99}^{25}\left(1-{0.99}\right)^{0} \]

densidad.binomialnegativa(x = 0, size = 25, prob = 0.99)

## [1] 0.7778214

dnbinom(x = 0, size = 25, prob = 0.99)

## [1] 0.7778214

2. calcue la probabilidad que el viajero al contactar al quincuagésimo ser humano haya contagiado a veintícinco personas cualquiera sea la configuración de las personas que han sido contagiadas, es decir, que halla contactado y no contagiado a veinticinco personas al contagiar al mismo número de personas.

\[ P(\text{# no contagiados hasta que se dan veintinco contagios}=25)=\binom{25+25-1}{25}{0.99}^{25}\left(1-{0.99}\right)^{25} \]

densidad.binomialnegativa(x = 25, size = 25, prob = 0.99)

## [1] 4.916243e-37

dnbinom(x = 25, size = 25, prob = 0.99)

## [1] 4.916243e-37

3. calcue la probabilidad de que el viajero haya tenido cincuenta contactos no contagiados al contagiar a veintícinco personas, o en otras palabras, la persona número setenta y cinco en ser contactada por el viajero haya sido la número veinticínco en ser contagiada.

\[ P(\text{# no contagiados hasta que se dan veintinco contagios}=50)=\binom{50+25-1}{50}{0.99}^{25}\left(1-{0.99}\right)^{50} \]

densidad.binomialnegativa(x = 50, size = 25, prob = 0.99)

## [1] 1.363483e-81

dnbinom(x = 50, size = 25, prob = 0.99)

## [1] 1.363483e-81

4. halle la probabilidad de que el viajero haya tenido setenta y cinco contactos no contagiados al contagiar a veintícinco personas, o en otras palabras, la persona número cien en ser contactada por el viajero haya sido la número veinticínco en ser contagiada.

\[ P(\text{# no contagiados hasta que se dan veintinco contagios}=75)=\binom{75+25-1}{75}{0.99}^{25}\left(1-{0.99}\right)^{75} \]

densidad.binomialnegativa(x = 75, size = 25, prob = 0.99)

## [1] 4.715917e-128

dnbinom(x = 75, size = 25, prob = 0.99)

## [1] 4.715917e-128

5. halle la probabilidad de que el viajero haya tenido cien contactos no contagiados al contagiar a veintícinco personas, o en otras palabras, la persona número ciento veinticinco en ser contactada por el viajero haya sido la número veinticínco en ser contagiada.

\[ P(\text{# no contagiados hasta que se dan veintinco contagios}=100)=\binom{100+25-1}{100}{0.99}^{25}\left(1-{0.99}\right)^{100} \]

densidad.binomialnegativa(x = 100, size = 25, prob = 0.99)

## [1] 2.023189e-175

dnbinom(x = 100, size = 25, prob = 0.99)

## [1] 2.023189e-175

6. Si usted fuera el Ministro de salud de dicho país, estaría interesado en saber la probabilidad de que un viajero en partícular no contagie a otros con quien tuvo contacto; estudie dicha situación.

\[ P(\text{# contactos previos al primer contagio}>x)=1-P(\text{# contactos previos al primer contagio}{\leq}x) \]

\[ 1-P(\text{# no contagiados hasta que se dan veintinco contagios}{\leq}x)=1-\sum_{x=0}^{x}P(\text{# no contagiados hasta que se dan veintinco contagios}=x) \]

\[ 1-\sum_{x=0}^{x}P(\text{# no contagiados hasta que se dan veintinco contagios}=x)=1-\sum_{x=0}^{x}\binom{x+25-1}{x}{0.99}^{25}\left(1-{0.99}\right)^{x} \]

ensayos <- 100
probabilidad.de.contagio <- 0.99
valores.posibles <- 0:ensayos
fx.binomialnegativa <- dnbinom(x = valores.posibles, size = 25, prob = 0.99)
fbg = pnbinom(q = valores.posibles, size = 25, prob = 0.99, lower.tail = TRUE)
fx.binomialnegativa = data.frame(x = valores.posibles, "f(X)" = fx.binomialnegativa, "F(X)" = fbg)
#head(fx.binomialnegativa)

library(ggplot2)
p3 <- ggplot(fx.binomialnegativa, aes(x = x, y = f.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(ggplot2)
p4 <- ggplot(fx.binomialnegativa, aes(x = x, y = F.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Distribución F(x)")

library(gridExtra)
grid.arrange(p3, p4, nrow = 2, top = "Distribución binomial negativa")

7. El número esperado de contactos no contagiados al haber contagiado a veinticinco personas es:

esperanza.binomialnegativa(size = 25, prob = 0.99)

## [1] 0.2525253

25*(1-0.99)/0.99

## [1] 0.2525253

8. La distancia promedio, elevada al cuadrado, al valor esperado del número de contactos no contagiados al haber contagiado a veinticinco seres humanos es igual a:

varianza.binomialnegativa(size = 25, prob = 0.99)

## [1] 0.255076

25*(1-0.99)/0.99**2

## [1] 0.255076

Simulación de la distribución binomial negativa

Sea \(X\) una variable aleatoria discreta para la cual \({P}({X}={x}_{j})={p}_{j}\), y sea \({U}\) una variable aleatoria continua con distribución uniforme en el intervalo \((0,1)\) \({U}{\sim}U(0,1)\) y además:

\[ {P(X=x)} = \begin{cases} \binom{x+size-1}{x}{prob}^{size}\left(1-{prob}\right)^{x}\text{ con }0{\leq}{prob}{\leq}1\text{ si }x=0,1,\ldots\\ 0\text{ en otro caso. } \end{cases} \] La distribución binomial negativa de parámetros k y p puede simularse a partir de la distribución geométrica.

\[ {X}_{i}{\sim}Geom(1, prob){\implies}{X}=\sum_{x=1}^{size}{X}_{i}{\sim}BN(size, prob) \]

Ejemplos

• Simular una variable aleatoria discreta que tiene distribución binomial con parámetros \(size=2\) y \(prob=0.7\)

prob = 0.7
U1 <- runif(n = 100000, min = 0, max = 1)
X1 <- floor(log(x = U1, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U2 <- runif(n = 100000, min = 0, max = 1)
X2 <- floor(log(x = U2, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
X = X1 + X2
tabla.frecuencias <- table(x = X)
simula.binomial.negativa <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.binomial.negativa, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.binomial.negativa <- as.data.frame(cbind(x = 0:max(x = X), Freq = dnbinom(x = 0:max(x = X), size = 2, prob = prob)))

library(ggplot2)
p2 <- ggplot(data = variable.binomial.negativa, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución binomial negativa")

• Simular una variable aleatoria discreta que tiene distribución geométrica con parámetros \(size=6\) y \(prob=0.4\)

prob = 0.4
U1 <- runif(n = 100000, min = 0, max = 1)
X1 <- floor(log(x = U1, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U2 <- runif(n = 100000, min = 0, max = 1)
X2 <- floor(log(x = U2, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U3 <- runif(n = 100000, min = 0, max = 1)
X3 <- floor(log(x = U3, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U4 <- runif(n = 100000, min = 0, max = 1)
X4 <- floor(log(x = U4, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U5 <- runif(n = 100000, min = 0, max = 1)
X5 <- floor(log(x = U5, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
U6 <- runif(n = 100000, min = 0, max = 1)
X6 <- floor(log(x = U6, base = exp(x = 1)) / log(x = 1 - prob, base = exp(x = 1)))
X <- X1 + X2 + X3 + X4 + X5 + X6
tabla.frecuencias <- table(x = X)
simula.binomial.negativa <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.binomial.negativa, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.binomial.negativa <- as.data.frame(cbind(x = 0:max(x = X), Freq = dnbinom(x = 0:max(x = X), size = 6, prob = prob)))

library(ggplot2)
p2 <- ggplot(data = variable.binomial.negativa, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución binomial")

Ejercicios

• Simular una variable aleatoria discreta que tiene distribución binomial negativa con parámetros \(size=1\) y \(prob=0.9\)

• Simular una variable aleatoria discreta que tiene distribución binomial negativa con parámetros \(size=2\) y \(prob=0.2\)

• Simular una variable aleatoria discreta que tiene distribución binomial negativa con parámetros \(size=6\) y \(prob=1\)

• Simular una variable aleatoria discreta que tiene distribución binomial negativa con parámetros \(size=12\) y \(prob=0.6\)

• Simular una variable aleatoria discreta que tiene distribución binomial negativa con parámetros \(size=17\) y \(prob=0.6\)

Distribución de probabilidad hipergeométrica \(X{\sim}H(m,m+n)\)

Especialmente útil en todos aquellos casos en los que se extraigan muestras en donde se extraen secuencialmente elementos sin devolución del elemento extraído o, en otras palabras, sin retornar a la situación experimental inicial como es el caso de la binomial.

Función de probabilidades \(f_{X}\)

\[ P(X=x)=\frac{\binom{m}{x}\binom{n}{k-x}}{\binom{m+n}{k}}\text{ con }{0}{\leq}k{\leq}{n+m}\text{ y }\max{\left\{{0},{k}-{n}\right\}}{\leq}x{\leq}\min{\left\{{k},{m}\right\}} \]

densidad.hipergeometrica <- function(x, m, n, k) choose(m,x)*choose(n,k-x)/choose(m+n,k)

Esperanza \(\mu_{x}\)

\[ \mu_{x}=\sum_{x=0}^{\infty}x{\cdot}\frac{\binom{m}{x}\binom{n}{k-x}}{\binom{m+n}{k}}=k{\cdot}\frac{m}{m+n} \]

esperanza.hipergeometrica <- function(m, n, k) sum((0:k)*choose(m,(0:k))*choose(n,k-(0:k))/choose(m+n,k))

Varianza \(\sigma_{x}^{2}\)

\[ \sigma_{x}^{2}=\sum_{x=0}^{\infty}\left(x-\mu_{x}\right)^{2}{\cdot}\frac{\binom{m}{x}\binom{n}{k-x}}{\binom{m+n}{k}}=\frac{m+n-k}{m+n-1}{\cdot}k{\cdot}\frac{m}{m+n}\left(1-\frac{m}{m+n}\right) \]

varianza.hipergeometrica <- function(m, n, k) sum((0:k-esperanza.hipergeometrica(m, n, k))**2*choose(m,(0:k))*choose(n,k-(0:k))/choose(m+n,k))

Ejemplo

Una clínica ha adquirido quinientas pruebas rápidas para COVID-19, se seleccionan aleatoriamente ochenta de ellas y se someten a una prueba para encontrar posibles defectos. Si el fabricante en la anterior ocasion envío dieciséis de otras 500 pruebas defectuosas probadas la anterior ocasión.

1. ¿Cuál es la probabilidad de que la muestra contenga a lo más cinco pruebas defectuosas?

\[ P\left(\text{# de pruebas defectuosas}{\leq}5\right) = P\left(\text{# de pruebas defectuosas}=0\right) + \cdots + P\left(\text{# de pruebas defectuosas}=5\right) \]

sum(densidad.hipergeometrica(x=0:5, m = 16, n = 500 - 16, k = 80))

## [1] 0.9707682

sum(dhyper(x=0:5, m = 16, n = 500 - 16, k = 80))

## [1] 0.9707682

\[ P\left(\text{# de pruebas defectuosas}{\leq}5\right)=\sum_{x=0}^{5}{P\left(\text{# de pruebas defectuosas}=x\right)} \]

phyper(q = 5, m = 16, n = 500 - 16, k = 80, lower.tail = TRUE)

## [1] 0.9707682

2. ¿Cuál es la probabilidad de que la muestra contenga mas de cinco pruebas defectuosos?

\[ P\left(\text{# de pruebas defectuosas}>5\right) = \sum_{x=6}^{80}{P\left(\text{# de pruebas defectuosas}{\leq}x\right)} \]

sum(densidad.hipergeometrica(x=6:16, m = 16, n = 500 - 16, k = 80))

## [1] 0.02923178

sum(dhyper(x=6:16, m = 16, n = 500 - 16, k = 80))

## [1] 0.02923178

phyper(q = 5, m = 16, n = 500 - 16, k = 80, lower.tail = FALSE)

## [1] 0.02923178

La distribución de probabilidades del número de pruebas defectuosas sería:

extracciones <- 80
pruebas.malas <- 16
pruebas.buenas <- 500 - 16
valores.posibles <- max(0, extracciones - pruebas.buenas): min(extracciones, pruebas.malas)
fx.hipergeometrica <- dhyper(x = valores.posibles, m = pruebas.malas, n = pruebas.buenas, k = extracciones)
fbh = phyper(q = valores.posibles, m = pruebas.malas, n = pruebas.buenas, k = extracciones, lower.tail = TRUE)
fx.hipergeometrica = data.frame(x = valores.posibles, "f(X)" = fx.hipergeometrica, "F(X)" = fbh)
#tail(fx.hipergeometrica)

library(ggplot2)
p5 <- ggplot(fx.hipergeometrica, aes(x = x, y = f.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(ggplot2)
p6 <- ggplot(fx.hipergeometrica, aes(x = x, y = F.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Distribución F(x)")

library(gridExtra)
grid.arrange(p5, p6, nrow = 2, top = "Distribución hipergeométrica")

3. El número esperado de pruebas de mala calidad es igual a:

esperanza.hipergeometrica(m = 16,n = 500-16, k = 80)

## [1] 2.56

80*16/(16+500-16)

## [1] 2.56

4. La distancia promedio al valor esperado, elevada al cuadrado, del número de pruebas de mala calidad es igual a:

varianza.hipergeometrica(m = 16,n = 500-16, k = 80)

## [1] 2.085759

(16+500-16-80)/(16+500-16-1)*80*16/(16+500-16)*(1-16/(16+500-16))

## [1] 2.085759

Simulación de la distribución hipergeométrica

Una variable aleatoria discreta con distribución poisson con parámetros \(m\): cantidad de elementos del cual se desea extraer un número determinado de ellos, \(n\): resto de elementos que sumado al grupo del cual se extraen los elementos conforman el total en la población o \(m+n\) y \(k\): el tamaño de la muestra aleatoria a ser seleccionada sin reemplazo puede ser generada a partir de la simulación de \(size\rightarrow\infty\) variables aleatorias \(U(0,1)\) mediante el método de la transformación inversa; a través de la función de distribución de probabilidad acumulada \(F(\cdot)\).

\[ {X}{\sim}H(m,m+n) \]

Ejemplos

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetros: grupo del cual se desean extraer \(X\) elementos (variable aleatoria) \(m = 10\), resto de elementos en la población \(n = 15\) de un total de \(m + n = 10 + 15\) elementos y \(k = 7\) elementos seleccionados aleatoriamente.

U <- runif(100000)

m = 10
n = 15
k = 7

F <- phyper(0:k, m = m, n = n, k = k)

X <- c()
for(i in 1:100000){
  X[i] <- min(which(U[i] < F))-1
}
tabla.frecuencias <- table(x = X)
simula.hipergeometrica <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.hipergeometrica, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.hipergeometrica <- as.data.frame(cbind(x = 0:max(X), Freq = dhyper(x = 0:max(X), m = m, n = n, k = k)))

library(ggplot2)
p2 <- ggplot(data = variable.hipergeometrica, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución hipergeométrica")

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetros: grupo del cual se desean extraer \(X\) elementos (variable aleatoria) \(m = 15\), resto de elementos en la población \(n = 10\) de un total de \(m + n = 15 + 10\) elementos y \(k = 7\) elementos seleccionados aleatoriamente.

U <- runif(100000)

m = 15
n = 10
k = 7

F <- phyper(0:k, m = m, n = n, k = k)

X <- c()
for(i in 1:100000){
  X[i] <- min(which(U[i] < F))-1
}
tabla.frecuencias <- table(x = X)
simula.hipergeometrica <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.hipergeometrica, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.hipergeometrica <- as.data.frame(cbind(x = 0:max(X), Freq = dhyper(x = 0:max(X), m = m, n = n, k = k)))

library(ggplot2)
p2 <- ggplot(data = variable.hipergeometrica, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución hipergeométrica")

Ejercicios

• Simular una variable aleatoria discreta que tiene distribución hipergeométrica con parámetro \(m=80\), \(n=53\) y \(k=88\)

• Simular una variable aleatoria discreta que tiene distribución hipergeométrica con parámetro \(m=95\), \(n=91\) y \(k=100\)

• Simular una variable aleatoria discreta que tiene distribución hipergeométrica con parámetro \(m=66\), \(n=69\) y \(k=51\)

• Simular una variable aleatoria discreta que tiene distribución hipergeométrica con parámetro \(m=71\), \(n=23\) y \(k=49\)

• Simular una variable aleatoria discreta que tiene distribución hipergeométrica con parámetro \(m=55\), \(n=41\) y \(k=41\)

Distribución de probabilidad poisson \(X{\sim}P(\lambda)\)

Es una distribucion discreta al igual que la distribucion binomial, y solo tiene un solo parametro de la media que se le denomina (lambda o \(\lambda\)).

Función de probabilidades \(f_{X}\)

\[ P(X=x)=\frac{e^{-lambda}lambda^{x}}{x!}\text{ con }lambda{=}\mu_{x}{=}\sigma_{x}^{2}\text{ y }x=0,1,2,\ldots \]

densidad.poisson <- function(x, lambda) exp(-lambda)*lambda**x/factorial(x)

Esperanza \(\mu_{x}\)

\[ \mu_{x}=\sum_{x=0}^{\infty}x{\cdot}\frac{e^{-lambda}lambda^{x}}{x!}=lambda \]

esperanza.poisson <- function(size, lambda) sum((0:size)*exp(-lambda)*lambda**(0:size)/factorial(0:size))

Varianza \(\sigma_{x}^{2}\)

\[ \sigma_{x}^{2}=\sum_{x=0}^{\infty}\left(x-\mu_{x}\right)^{2}{\cdot}\frac{e^{-lambda}lambda^{x}}{x!}=lambda \]

varianza.poisson <- function(size, lambda) sum((0:size-esperanza.poisson(size, lambda))**2*exp(-lambda)*lambda**(0:size)/factorial(0:size))

Ejemplo

La tasa o índice de letalidad se refiere al cociente de fallecimientos en relación a las peronas que se han contagiado de COVID-19 y cuyo resultado suele mostrase en porcentaje. Si en cierto país la tasa de letalidad para el primer més de aislamiento ha sido de 1830 muertes entre 36210, es decir, una tasa de letalidad de \(1830/36210\); lo cual es aproximadamente el 5% o, en terminos diarios, 61 muertos por cada 1207 enfermos. Con base en lo anteriormente expuesto calcule:

1. la probabilidad de que haya más de un sesenta y un muertes diarias entre las personas contagiadas.

\[ P\left(\text{# de muertos}>61\right)=1-P\left(\text{# de muertos}{\leq}61\right) \]

\[ 1-P\left(\text{# de muertos}{\leq}61\right)=1-\sum_{x = 0}^{61}{P\left(\text{# de muertos}=x\right)} \]

\[ 1-\sum_{x = 0}^{61}{P\left(\text{# de muertos}=x\right)}=1-\sum_{x = 0}^{61}\frac{e^{-61}61^{x}}{x!} \]

1 - sum(densidad.poisson(x = 0:61, lambda = 61))

## [1] 0.4660183

1 - sum(dpois(x = 0:61, lambda = 61))

## [1] 0.4660183

\[ P\left(\text{# de infectados}>61\right)=\sum_{x = 62}^{\infty}{P\left(\text{# de infectados}=x\right)} \]

ppois(q = 61, lambda = 61, lower.tail = FALSE)

## [1] 0.4660183

2. ¿cuál sería en dicha ciudad el número esperado de enfermos que morirán diariamente de mantenerse la medida de aislamiento?

esperanza.poisson(size=120, lambda = 61)

## [1] 61

61

## [1] 61

En otas palabras, aproximadamente 61 por cada mil doscientos siete habitantes

1. Mediante la función de R

personas.enfermas <- 32610
muertes.esperadas <- 1820
valores.posibles <- 0:personas.enfermas
fx <- function(x) x*dpois(x = x, lambda = muertes.esperadas)
Ex <- sum(fx(valores.posibles))
Ex

## [1] 1820

2. Mediante simulación

muertos.simulados <- as.data.frame(rpois(n = 3261, lambda = 182))
mean(muertos.simulados$`rpois(n = 3261, lambda = 182)`)

## [1] 181.7764

3. Estudie la probabilidad de que dado un número determinado de enfermos, el actual, muera un número significativo de estos el siguiente mes si las condiciones se mantienen como hasta ahora.

personas.enfermas <- 32610
muertes.esperadas <- 1820
valores.posibles <- 0:personas.enfermas
fx.poisson <- dpois(x = valores.posibles, lambda = muertes.esperadas)
fbp = ppois(q = valores.posibles, lambda = muertes.esperadas, lower.tail = TRUE)
fx.poisson = data.frame(x = valores.posibles, "f(X)" = fx.poisson, "F(X)" = fbp)
#head(fx.poisson)

library(ggplot2)
p7 <- ggplot(fx.poisson, aes(x = x, y = f.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(ggplot2)
p8 <- ggplot(fx.poisson, aes(x = x, y = F.X.)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Distribución F(x)")

library(gridExtra)
grid.arrange(p7, p8, nrow = 2, top = "Distribución poisson")

Simulación de la distribución poisson

Una variable aleatoria discreta con distribución poisson con parámetros \(\lambda\) puede ser generada a partir de la simulación de \(size\rightarrow\infty\) variables aleatorias \(U(0,1)\) mediante el método de la transformación inversa; a través de la función de distribución de probabilidad acumulada \(F(\cdot)\).

\[ {X}{\sim}Pois(\lambda) \]

Ejemplos

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=1.5\)

library(car)
U <- runif(n = 100000, min = 0, max = 1)
CL <- "0:0.2231302=0; 
       0.2231303:0.5578254=1;
       0.5578255:0.8088468=2;
       0.8088469:0.9343575=3; 
       0.9343576:0.9814241=4;
       0.9814242:0.9955440=5;
       0.9955441:0.9990740=6;
       0.9990741:0.9998304=7;
       0.9998305:0.9999723=8;
       0.9999724:0.9999959=9;
       0.9999960:1=10"
X <- car::recode(var = U, recodes = CL)
tabla.frecuencias <- table(x = X)
simula.poisson <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.poisson, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.poisson <- as.data.frame(cbind(x = 0:max(X), Freq = dpois(x = 0:max(X), lambda = 1.5)))

library(ggplot2)
p2 <- ggplot(data = variable.poisson, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución poisson")

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=3.0\)

lambda = 3.0
U <- runif(100000)

X <- c(0:20)
P <- c()
P[1] = exp(-lambda)
F <- c()
F[1] <- P[1]
for(i in 2:length(X)){
  P[i] = (lambda*P[i-1])/X[i]
  F[i] <- F[i-1]+P[i]
}

X <- c()
for(i in 1:100000){
  X[i] <- min(which(U[i] < F))-1
}
tabla.frecuencias <- table(x = X)
simula.poisson <- as.data.frame(x = tabla.frecuencias / 100000)

library(ggplot2)
p1 <- ggplot(data = simula.poisson, aes(x = x, y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

variable.poisson <- as.data.frame(cbind(x = 0:max(X), Freq = dpois(x = 0:max(X), lambda = lambda)))

library(ggplot2)
p2 <- ggplot(data = variable.poisson, aes(x = as.factor(x), y = Freq)) +
  geom_bar(aes(color = x), stat = "identity", fill = "white") + theme(legend.position = "none") + xlab("") + ylab("Densidad f(x)")

library(gridExtra)
grid.arrange(p1, p2, ncol = 2, top = "Distribución poisson")

Ejercicios

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=50\)

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=62\)

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=42\)

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=62\)

• Simular una variable aleatoria discreta que tiene distribución poisson con parámetro \(lambda=70\)