Laboratorio 14: Modelos Probabilísticos Discretos en R

Objetivo General: Aprender a utilizar R y sus funciones estadísticas base (junto con ggplot2 para gráficas) para modelar, calcular y graficar las distribuciones de probabilidad discretas más importantes: Bernoulli, Binomial, Binomial Negativa, Poisson e Hipergeométrica.

Herramientas a Utilizar: - RStudio o Google Colab (kernel de R) - R 4.x - Paquetes: ggplot2, dplyr


Parte 1: Preparación del Entorno

Antes de comenzar, debemos cargar las librerías que nos permitirán realizar los cálculos y las gráficas.

Código para ejecutar:

# install.packages(c("ggplot2", "dplyr"))  # Ejecutar solo la primera vez

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
cat("¡Entorno listo para el laboratorio!\n")
## ¡Entorno listo para el laboratorio!

Parte 2: Función de Ayuda para Graficar

Para no repetir código, crearemos una función que nos permita graficar cualquier distribución discreta de manera sencilla. Esta función recibirá los valores posibles (eje x), sus probabilidades (eje y) y, opcionalmente, qué barras resaltar.

Código para ejecutar:

# @title **A. Función para Graficar Distribuciones Discretas**
plot_distribution <- function(values, pmf, title, highlight = NULL,
                              highlight_color = "red", label = NULL) {
  df <- data.frame(k = values, prob = pmf)

  p <- ggplot(df, aes(x = factor(k), y = prob)) +
    geom_bar(stat = "identity", fill = "skyblue", color = "black") +
    labs(title = title, x = "Valores (k)", y = "Probabilidad P(X=k)") +
    theme_minimal(base_size = 14) +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"),
          axis.text.x = element_text(size = 9))

  if (!is.null(highlight)) {
    df_high <- df %>% filter(k %in% highlight)
    p <- p + geom_bar(data = df_high, aes(x = factor(k), y = prob),
                      stat = "identity", fill = highlight_color, color = "black")
  }

  if (!is.null(label)) {
    p <- p + annotate("text", x = length(values)/2, y = max(pmf) * 0.9,
                      label = label, size = 5, color = "darkgreen")
  }

  print(p)
}

Parte 3: Exploración de Modelos Probabilísticos Discretos

A continuación, exploraremos cada distribución con un ejemplo práctico.

3.1. Distribución de Bernoulli

Definición: Modela un experimento con dos únicos resultados posibles: “éxito” (1) o “fracaso” (0). La probabilidad de éxito se denota como p.

  • Función de Probabilidad: \(P(X=k) = p^k (1-p)^{1-k}\), para \(k = 0, 1\).
  • Media (Esperanza): \(E(X) = p\)
  • Varianza: \(Var(X) = p(1-p)\)

Ejemplo: Lanzamiento de una moneda una vez. Sea X = 1 si sale cara, 0 si sale sello. Sea \(p=0.5\) la probabilidad de obtener cara.

Código para ejecutar:

# @title **B1. Ejemplo de Distribución de Bernoulli**
n_bernoulli <- 1   # Un solo ensayo
p_bernoulli <- 0.5 # Probabilidad de éxito (cara)

x_bernoulli <- 0:n_bernoulli
pmf_bernoulli <- dbinom(x_bernoulli, n_bernoulli, p_bernoulli)

plot_distribution(x_bernoulli, pmf_bernoulli,
                  sprintf("Distribución de Bernoulli_roiver (p=%.1f)", p_bernoulli))

cat(sprintf("Media (Esperanza): %.2f\n", n_bernoulli * p_bernoulli))
## Media (Esperanza): 0.50
cat(sprintf("Varianza: %.2f\n", n_bernoulli * p_bernoulli * (1 - p_bernoulli)))
## Varianza: 0.25

3.2. Distribución Binomial

Definición: Modela el número de éxitos en n ensayos independientes de Bernoulli. Cada ensayo tiene una probabilidad de éxito p.

  • Función de Probabilidad: \(P(X=k) = \binom{n}{k} p^k (1-p)^{n-k}\)
  • Media (Esperanza): \(E(X) = n \cdot p\)
  • Varianza: \(Var(X) = n \cdot p \cdot (1-p)\)

Ejemplo: Lanzar una moneda 10 veces. Sea X = número de caras. \(n=10\), \(p=0.5\).

Código para ejecutar:

# @title **B2. Ejemplo de Distribución Binomial**
n_binom <- 10
p_binom <- 0.5

x_binom <- 0:n_binom
pmf_binom <- dbinom(x_binom, n_binom, p_binom)

plot_distribution(x_binom, pmf_binom,
                  sprintf("Distribución Binomial_roiver (n=%d, p=%.1f)", n_binom, p_binom))

cat(sprintf("Media (Esperanza): %.2f\n", n_binom * p_binom))
## Media (Esperanza): 5.00
cat(sprintf("Varianza: %.2f\n", n_binom * p_binom * (1 - p_binom)))
## Varianza: 2.50

Ejemplo 1 (del PDF): Supongamos que lanzamos una moneda justa (probabilidad de éxito p=0.5), n=10 veces. ¿Cuál es la probabilidad de obtener exactamente 6 caras?

# @title **Ejemplo 1. Probabilidad de obtener exactamente 6 caras**
n <- 10; p <- 0.5; k <- 6

prob_binom_1 <- dbinom(k, n, p)
cat(sprintf("Probabilidad de obtener exactamente 6 caras: P(X = 6) = %.5f\n", prob_binom_1))
## Probabilidad de obtener exactamente 6 caras: P(X = 6) = 0.20508
plot_distribution(0:n, dbinom(0:n, n, p),
                  "Distribución Binomial -roiver P(X=6)",
                  highlight = k,
                  label = sprintf("P(X = 6) = %.5f", prob_binom_1))

Ejemplo 2 (del PDF): La probabilidad de que una pieza fabricada por una empresa sea defectuosa es 0.1. Halla la probabilidad de que en una muestra de 25 piezas se encuentren a lo sumo cuatro piezas defectuosas. Haz un gráfico y sombrea P(X <= 4) cada rectángulo con rojo y dentro del gráfico coloca P(X <= 4).

# @title **Ejemplo 2. Probabilidad de a lo sumo 4 piezas defectuosas**
n <- 25; p <- 0.1

probabilidades <- dbinom(0:n, n, p)
probabilidad_acumulada <- pbinom(4, n, p)

plot_distribution(0:n, probabilidades,
                  "Distribución Binomial - P(X <= 4)",
                  highlight = 0:4,
                  label = sprintf("P(X <= 4) = %.5f", probabilidad_acumulada))

cat(sprintf("La probabilidad de que en una muestra de 25 piezas se encuentren a lo sumo 4 defectuosas es: %.5f\n",
            probabilidad_acumulada))
## La probabilidad de que en una muestra de 25 piezas se encuentren a lo sumo 4 defectuosas es: 0.90201

3.3. Distribución Binomial Negativa

Definición: Modela el número de ensayos de Bernoulli necesarios para obtener un número fijo de éxitos (r).

  • Función de Probabilidad: \(P(X=k) = \binom{k-1}{r-1} p^r (1-p)^{k-r}\), para \(k = r, r+1, \dots\)
  • Media (Esperanza): \(E(X) = r / p\)
  • Varianza: \(Var(X) = r(1-p) / p^2\)

Nota de R: dnbinom(x, size=r, prob=p) devuelve la probabilidad de obtener x fracasos antes del r-ésimo éxito (misma convención que scipy.stats.nbinom).

Ejemplo 1 (del PDF): ¿Cuál es la probabilidad de obtener 5 éxitos antes de 3 fracasos si la probabilidad de éxito es p = 0.7?

# @title **Ejemplo 1. Probabilidad de 5 éxitos antes de 3 fracasos**
r <- 5     # Número de éxitos
p <- 0.7   # Probabilidad de éxito
k <- 3     # Número de fracasos antes de obtener 5 éxitos

prob_nbinom_1 <- dnbinom(k, size = r, prob = p)
cat(sprintf("Probabilidad de obtener 5 éxitos antes de 3 fracasos: P(X=5) = %.5f\n",
            prob_nbinom_1))
## Probabilidad de obtener 5 éxitos antes de 3 fracasos: P(X=5) = 0.15883
# Momentos
esperanza_nbinom  <- r * (1 - p) / p
varianza_nbinom   <- r * (1 - p) / (p^2)
desviacion_nbinom <- sqrt(varianza_nbinom)

cat("\nMomentos en el Modelo Binomial Negativo\n")
## 
## Momentos en el Modelo Binomial Negativo
cat(sprintf("Esperanza: E(X) = %.2f\n", esperanza_nbinom))
## Esperanza: E(X) = 2.14
cat(sprintf("Varianza: V(X) = %.2f\n", varianza_nbinom))
## Varianza: V(X) = 3.06
cat(sprintf("Desviación estándar: %.2f\n", desviacion_nbinom))
## Desviación estándar: 1.75

Ejemplo 2 (del PDF): La probabilidad de un niño expuesto al COVID-19 se contagie es de 0.6. ¿Cuál es la probabilidad de que el 10° niño expuesto sea el tercero en contraerla?

# @title **Ejemplo 2. Probabilidad de que el 10° niño sea el tercero en contagiarse**
n <- 10; r <- 3; p <- 0.6
k <- n - r  # Número de fracasos antes de obtener 3 éxitos

prob_nbinom_2 <- dnbinom(k, size = r, prob = p)
cat(sprintf("Probabilidad de que el 10° niño sea el tercero en contagiarse: %.5f\n",
            prob_nbinom_2))
## Probabilidad de que el 10° niño sea el tercero en contagiarse: 0.01274

3.4. Distribución de Poisson

Definición: Modela el número de eventos que ocurren en un intervalo de tiempo o espacio, dado un promedio de ocurrencia (\(\lambda\)).

  • Función de Probabilidad: \(P(X=k) = \frac{\lambda^k e^{-\lambda}}{k!}\), para \(k = 0, 1, 2, \dots\)
  • Media (Esperanza): \(E(X) = \lambda\)
  • Varianza: \(Var(X) = \lambda\)

Ejemplo 1 (del PDF): Una central telefónica recibe en promedio 2 llamadas por minuto. ¿Cuál es la probabilidad de que en un minuto se reciban exactamente 3 llamadas?

# @title **Ejemplo 1. Probabilidad de recibir exactamente 3 llamadas**
lambda_poisson <- 2
k <- 3

prob_poisson_1 <- dpois(k, lambda_poisson)
cat(sprintf("Probabilidad de recibir exactamente 3 llamadas: P(X=3) = %.5f\n", prob_poisson_1))
## Probabilidad de recibir exactamente 3 llamadas: P(X=3) = 0.18045
plot_distribution(0:14, dpois(0:14, lambda_poisson),
                  "Distribución de Poisson -roiver  P(X=3)",
                  highlight = k,
                  label = sprintf("P(X=3) = %.5f", prob_poisson_1))

Ejemplo 2 (del PDF): En la zona que rodea un punto kilométrico de una autovía hay una población de erizos. Se ha realizado un estudio por el que se ha deducido que los erizos cruzan la autovía a razón de 4 cada 30 minutos. La probabilidad de que 6 o más erizos crucen la autovía en una hora es: Dibuja la distribución de Poisson y sombrea el rectángulo con P(X>=6).

# @title **Ejemplo 2. Probabilidad de 6 o más erizos en una hora**
lambda_poisson <- 8  # 4 erizos/30 min * 2 = 8 por hora
prob_6_o_mas <- 1 - ppois(5, lambda_poisson)  # P(X >= 6)

cat(sprintf("La probabilidad de que 6 o más erizos crucen la autovía en una hora es: %.5f\n",
            prob_6_o_mas))
## La probabilidad de que 6 o más erizos crucen la autovía en una hora es: 0.80876
plot_distribution(0:19, dpois(0:19, lambda_poisson),
                  "Distribución de Poisson - P(X >= 6)",
                  highlight = 6:19,
                  label = sprintf("P(X >= 6) = %.5f", prob_6_o_mas))


3.5. Distribución Hipergeométrica

Definición: Modela el número de éxitos en una muestra de tamaño n, extraída sin reemplazo de una población finita de tamaño N, que contiene K éxitos.

  • Función de Probabilidad: \(P(X=k) = \frac{\binom{K}{k} \binom{N-K}{n-k}}{\binom{N}{n}}\)
  • Media (Esperanza): \(E(X) = n \cdot \frac{K}{N}\)
  • Varianza: \(Var(X) = n \cdot \frac{K}{N} \cdot \left(1 - \frac{K}{N}\right) \cdot \frac{N-n}{N-1}\)

Nota de R: dhyper(x, m=K, n=N-K, k=n) — cuidado con los nombres de los argumentos (m = éxitos en la población, n = fracasos en la población, k = tamaño de la muestra).

Ejemplo 1 (del PDF): Una urna contiene 10 bolas rojas y 20 bolas azules. Se seleccionan 5 bolas sin reemplazo. ¿Cuál es la probabilidad de seleccionar exactamente 3 bolas rojas?

# @title **Ejemplo 1. Probabilidad de seleccionar exactamente 3 bolas rojas**
N <- 30   # Población total
K <- 10   # Éxitos en la población (bolas rojas)
n <- 5    # Tamaño de la muestra
k <- 3    # Éxitos en la muestra

prob_hyper_1 <- dhyper(k, m = K, n = N - K, k = n)
cat(sprintf("Probabilidad de seleccionar exactamente 3 bolas rojas: P(X=3) = %.5f\n",
            prob_hyper_1))
## Probabilidad de seleccionar exactamente 3 bolas rojas: P(X=3) = 0.15999
# Momentos
esperanza_hyper  <- n * (K / N)
varianza_hyper   <- n * (K/N) * (1 - K/N) * ((N - n) / (N - 1))
desviacion_hyper <- sqrt(varianza_hyper)

cat("\nMomentos en el Modelo Hipergeométrico\n")
## 
## Momentos en el Modelo Hipergeométrico
cat(sprintf("Esperanza: E(X) = %.2f\n", esperanza_hyper))
## Esperanza: E(X) = 1.67
cat(sprintf("Varianza: V(X) = %.2f\n", varianza_hyper))
## Varianza: V(X) = 0.96
cat(sprintf("Desviación estándar: %.2f\n", desviacion_hyper))
## Desviación estándar: 0.98
plot_distribution(0:n, dhyper(0:n, m = K, n = N - K, k = n),
                  "Distribución Hipergeométrica - P(X=3)",
                  highlight = k,
                  label = sprintf("P(X=3) = %.5f", prob_hyper_1))

Parte 4: Ejercicios Propuestos

Resuelve los siguientes problemas utilizando el código que has aprendido. Para cada uno, calcula la probabilidad y grafica la distribución.

  1. Binomial: Un estudiante responde un examen de 15 preguntas de opción múltiple (4 opciones cada una) al azar. ¿Cuál es la probabilidad de que acierte exactamente 6 preguntas? ¿Y de que acierte 5 o menos?
  2. Poisson: En un sitio web, se reciben en promedio 3 visitas por minuto. ¿Cuál es la probabilidad de recibir exactamente 5 visitas en un minuto? ¿Y de recibir 2 o menos?
  3. Hipergeométrica: Un lote de 50 piezas contiene 10 defectuosas. Se seleccionan 8 piezas al azar sin reemplazo. ¿Cuál es la probabilidad de encontrar exactamente 2 defectuosas en la muestra?
  4. Binomial Negativa: Una máquina produce piezas con una probabilidad de defecto del 5%. ¿Cuál es la probabilidad de que la tercera pieza defectuosa sea la décima producida?

Solución Ejercicio 1 (Binomial):

# @title **Ejercicio 1. Examen de 15 preguntas al azar**
n <- 15; p <- 1/4

prob_6 <- dbinom(6, n, p)
prob_5_o_menos <- pbinom(5, n, p)

cat(sprintf("P(X = 6) = %.5f\n", prob_6))
## P(X = 6) = 0.09175
cat(sprintf("P(X <= 5) = %.5f\n", prob_5_o_menos))
## P(X <= 5) = 0.85163
plot_distribution(0:n, dbinom(0:n, n, p),
                  "Binomial - Examen (n=15, p=0.25)",
                  highlight = 6,
                  label = sprintf("P(X=6)=%.5f", prob_6))

Solución Ejercicio 2 (Poisson):

# @title **Ejercicio 2. Visitas por minuto en un sitio web**
lambda <- 3

prob_5 <- dpois(5, lambda)
prob_2_o_menos <- ppois(2, lambda)

cat(sprintf("P(X = 5) = %.5f\n", prob_5))
## P(X = 5) = 0.10082
cat(sprintf("P(X <= 2) = %.5f\n", prob_2_o_menos))
## P(X <= 2) = 0.42319
plot_distribution(0:14, dpois(0:14, lambda),
                  "Poisson - Visitas por minuto (λ=3)",
                  highlight = 5,
                  label = sprintf("P(X=5)=%.5f", prob_5))

Solución Ejercicio 3 (Hipergeométrica):

# @title **Ejercicio 3. Piezas defectuosas en un lote**
N <- 50; K <- 10; n <- 8; k <- 2

prob_2_defectuosas <- dhyper(k, m = K, n = N - K, k = n)
cat(sprintf("P(X = 2 defectuosas) = %.5f\n", prob_2_defectuosas))
## P(X = 2 defectuosas) = 0.32172
plot_distribution(0:n, dhyper(0:n, m = K, n = N - K, k = n),
                  "Hipergeométrica - 2 defectuosas en muestra de 8",
                  highlight = k,
                  label = sprintf("P(X=2) = %.5f", prob_2_defectuosas))

Solución Ejercicio 4 (Binomial Negativa):

# @title **Ejercicio 4. Tercera pieza defectuosa en la décima producida**
r <- 3       # 3 defectuosas (éxitos)
p <- 0.05    # Probabilidad de defecto
k <- 10 - 3  # Fracasos antes del 3er éxito

prob_3ra_defectuosa <- dnbinom(k, size = r, prob = p)
cat(sprintf("P(la 3ª defectuosa sea la 10ª pieza) = %.5f\n", prob_3ra_defectuosa))
## P(la 3ª defectuosa sea la 10ª pieza) = 0.00314
plot_distribution(0:60 + r, dnbinom(0:60, size = r, prob = p),
                  "Binomial Negativa - Posición de la 3ª defectuosa (p=0.05)",
                  highlight = 10,
                  label = sprintf("P(X=10) = %.5f", prob_3ra_defectuosa))