Monte Carlo

Modelo Beta-Binomial

1. Planteamiento del Modelo

Sea \(Y\) el número de éxitos en \(n\) ensayos Bernoulli independientes con probabilidad de éxito \(\theta\).

La función de verosimilitud es:

\[ Y \mid \theta \sim \text{Binomial}(n, \theta) \]

\[ p(y \mid \theta) = \binom{n}{y} \theta^y (1-\theta)^{n-y} \]

2. Distribución a priori

En el enfoque bayesiano, se asigna una distribución previa a \(\theta\).
El prior conjugado para la Binomial es la distribución Beta:

\[ \theta \sim \text{Beta}(\alpha, \beta) \]

con densidad:

\[ p(\theta) = \frac{\Gamma(\alpha+\beta)} {\Gamma(\alpha)\Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1} \]

donde:

  • \(\alpha > 0\)
  • \(\beta > 0\)

3. Distribución Posterior

Por el Teorema de Bayes:

\[ p(\theta \mid y) \propto p(y \mid \theta)p(\theta) \]

Multiplicando verosimilitud y prior:

\[ p(\theta \mid y) \propto \theta^{y+\alpha-1}(1-\theta)^{n-y+\beta-1} \]

Por lo tanto, la distribución posterior es:

\[ \theta \mid y \sim \text{Beta}(\alpha + y,\; \beta + n - y) \]

Esto muestra la propiedad de conjugación.

4. Momentos de la Distribución Beta

Para una variable aleatoria \(\theta \sim \text{Beta}(\alpha, \beta)\):

Esperanza

\[ E(\theta) = \frac{\alpha}{\alpha+\beta} \]

Varianza

\[ \text{Var}(\theta) = \frac{\alpha\beta} {(\alpha+\beta)^2(\alpha+\beta+1)} \]

Desviación estándar

\[ \text{DE}(\theta) = \sqrt{ \frac{\alpha\beta} {(\alpha+\beta)^2(\alpha+\beta+1)} } \]

5. Intervalo Creíble

Un intervalo creíble del 95% para \(\theta\) se obtiene a partir de los cuantiles de la distribución posterior:

\[ IC_{95\%} = \left[ q_{0.025},\; q_{0.975} \right] \]

donde:

\[ q_p = F^{-1}(p) \]

y \(F^{-1}\) es la función cuantil de la distribución Beta posterior.

# ==========================
# MODELO BETA-BINOMIAL
# ==========================

# Librerías
library(ggplot2)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
# --------------------------
# Datos
# --------------------------
n <- 20
y <- 12

# Prior
alpha_prior <- 2
beta_prior  <- 2

# Posterior
alpha_post <- alpha_prior + y
beta_post  <- beta_prior + n - y

# --------------------------
# Funciones estadísticas
# --------------------------

# Esperanza
mean_prior <- alpha_prior / (alpha_prior + beta_prior)
mean_post  <- alpha_post  / (alpha_post  + beta_post)

# Varianza
var_prior <- (alpha_prior * beta_prior) /
  ((alpha_prior + beta_prior)^2 * (alpha_prior + beta_prior + 1))

var_post <- (alpha_post * beta_post) /
  ((alpha_post + beta_post)^2 * (alpha_post + beta_post + 1))

# Desviación estándar
sd_prior <- sqrt(var_prior)
sd_post  <- sqrt(var_post)

# Intervalos creíbles 95%
ic_prior <- qbeta(c(0.025, 0.975), alpha_prior, beta_prior)
ic_post  <- qbeta(c(0.025, 0.975), alpha_post, beta_post)

# --------------------------
# Mostrar resultados
# --------------------------

cat("===== PRIOR =====\n")
## ===== PRIOR =====
cat("Esperanza:", mean_prior, "\n")
## Esperanza: 0.5
cat("Varianza:", var_prior, "\n")
## Varianza: 0.05
cat("Desv.Est:", sd_prior, "\n")
## Desv.Est: 0.2236068
cat("IC 95%:", ic_prior, "\n\n")
## IC 95%: 0.09429932 0.9057007
cat("===== POSTERIOR =====\n")
## ===== POSTERIOR =====
cat("Esperanza:", mean_post, "\n")
## Esperanza: 0.5833333
cat("Varianza:", var_post, "\n")
## Varianza: 0.009722222
cat("Desv.Est:", sd_post, "\n")
## Desv.Est: 0.09860133
cat("IC 95%:", ic_post, "\n")
## IC 95%: 0.385419 0.7680858
# --------------------------
# Gráficas con ggplot2
# --------------------------

theta <- seq(0, 1, length.out = 1000)

df <- data.frame(
  theta = theta,
  Prior = dbeta(theta, alpha_prior, beta_prior),
  Posterior = dbeta(theta, alpha_post, beta_post)
)

df_long <- df %>%
  tidyr::pivot_longer(-theta,
                      names_to = "Distribucion",
                      values_to = "Densidad")

ggplot(df_long, aes(x = theta, y = Densidad, color = Distribucion)) +
  geom_line(size = 1.2) +
  geom_vline(xintercept = mean_post, linetype = "dashed") +
  labs(title = "Modelo Beta-Binomial",
       subtitle = "Distribución Prior vs Posterior",
       x = expression(theta),
       y = "Densidad") +
  theme_minimal() +
  scale_color_manual(values = c("blue", "red"))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Modelo Gamma-Poisson

1. Planteamiento del Modelo

Sea \(Y\) el número de eventos ocurridos en un intervalo fijo de tiempo o espacio.

Se asume que:

\[ Y \mid \lambda \sim \text{Poisson}(\lambda) \]

con función de probabilidad:

\[ p(y \mid \lambda) = \frac{\lambda^y e^{-\lambda}}{y!}, \quad y = 0,1,2,\dots \]

donde:

  • \(\lambda > 0\) es la tasa promedio de ocurrencia.

2. Distribución a priori

En el enfoque bayesiano, se asigna una distribución previa a \(\lambda\).

El prior conjugado para la distribución Poisson es la distribución Gamma:

\[ \lambda \sim \text{Gamma}(\alpha, \beta) \]

(donde usamos la parametrización por forma–tasa)

La densidad es:

\[ p(\lambda) = \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta \lambda} \]

donde:

  • \(\alpha > 0\) es el parámetro de forma
  • \(\beta > 0\) es el parámetro de tasa

3. Distribución Posterior

Por el Teorema de Bayes:

\[ p(\lambda \mid y) \propto p(y \mid \lambda) p(\lambda) \]

Multiplicando verosimilitud y prior:

\[ p(\lambda \mid y) \propto \lambda^{y+\alpha-1} e^{-(\beta+1)\lambda} \]

Por lo tanto:

\[ \lambda \mid y \sim \text{Gamma}(\alpha + y,\; \beta + 1) \]

Caso General (muestra de tamaño n)

Si:

\[ Y_1,\dots,Y_n \mid \lambda \sim \text{Poisson}(\lambda) \]

y definimos:

\[ S = \sum_{i=1}^n Y_i \]

entonces la posterior es:

\[ \lambda \mid \mathbf{y} \sim \text{Gamma}(\alpha + S,\; \beta + n) \]

4. Momentos de la Distribución Gamma

Para:

\[ \lambda \sim \text{Gamma}(\alpha,\beta) \]

Esperanza

\[ E(\lambda) = \frac{\alpha}{\beta} \]

Varianza

\[ \text{Var}(\lambda) = \frac{\alpha}{\beta^2} \]

Desviación estándar

\[ \text{DE}(\lambda) = \sqrt{\frac{\alpha}{\beta^2}} \]

5. Intervalo Creíble

Un intervalo creíble del 95% se obtiene mediante los cuantiles:

\[ IC_{95\%} = \left[ q_{0.025},\; q_{0.975} \right] \]

# =====================================
# MODELO GAMMA - POISSON
# =====================================

library(ggplot2)
library(dplyr)
library(tidyr)

# -----------------------------
# Datos
# -----------------------------
y <- c(4,5,3,6,4,5,7,2,6,4)

n <- length(y)
S <- sum(y)

# -----------------------------
# Prior
# -----------------------------
alpha_prior <- 2
beta_prior  <- 1

# -----------------------------
# Posterior
# -----------------------------
alpha_post <- alpha_prior + S
beta_post  <- beta_prior + n

# -----------------------------
# Momentos
# -----------------------------
mean_prior <- alpha_prior / beta_prior
var_prior  <- alpha_prior / beta_prior^2
sd_prior   <- sqrt(var_prior)

mean_post <- alpha_post / beta_post
var_post  <- alpha_post / beta_post^2
sd_post   <- sqrt(var_post)

# Intervalo creíble 95%
ic_post <- qgamma(c(0.025, 0.975),
                  shape = alpha_post,
                  rate  = beta_post)

cat("===== POSTERIOR =====\n")
## ===== POSTERIOR =====
cat("Esperanza:", mean_post, "\n")
## Esperanza: 4.363636
cat("Varianza:", var_post, "\n")
## Varianza: 0.3966942
cat("Desv.Est:", sd_post, "\n")
## Desv.Est: 0.6298367
cat("IC 95%:", ic_post, "\n")
## IC 95%: 3.217401 5.681821
# =====================================
# GRÁFICA 1: PRIOR vs POSTERIOR
# =====================================

lambda <- seq(0, 10, length.out = 1000)

df <- data.frame(
  lambda = lambda,
  Prior = dgamma(lambda, alpha_prior, beta_prior),
  Posterior = dgamma(lambda, alpha_post, beta_post)
)

df_long <- df %>%
  pivot_longer(-lambda,
               names_to = "Distribucion",
               values_to = "Densidad")

ggplot(df_long,
       aes(x = lambda,
           y = Densidad,
           color = Distribucion,
           fill = Distribucion)) +
  geom_line(size = 1.2) +
  geom_area(alpha = 0.2, position = "identity") +
  geom_vline(xintercept = mean_post,
             linetype = "dashed",
             size = 1) +
  labs(title = "Modelo Gamma–Poisson",
       subtitle = "Distribución Prior vs Posterior",
       x = expression(lambda),
       y = "Densidad") +
  theme_minimal(base_size = 15) +
  scale_color_manual(values = c("#2C7BB6", "#D7191C")) +
  scale_fill_manual(values = c("#2C7BB6", "#D7191C")) +
  theme(legend.position = "top")

# =====================================
# GRÁFICA 2: DISTRIBUCIÓN PREDICTIVA
# =====================================

# Predictiva es Binomial Negativa
x_vals <- 0:15

pred <- data.frame(
  x = x_vals,
  Prob = dnbinom(x_vals,
                 size = alpha_post,
                 prob = beta_post / (beta_post + 1))
)

ggplot(pred,
       aes(x = x,
           y = Prob)) +
  geom_col(fill = "#FDAE61",
           color = "black",
           width = 0.7) +
  labs(title = "Distribución Predictiva",
       subtitle = "Gamma–Poisson → Binomial Negativa",
       x = "Número futuro de llamadas",
       y = "Probabilidad") +
  theme_minimal(base_size = 15)

Modelo Beta–Geométrico

1. Planteamiento del Modelo

Sea \(Y\) el número de ensayos hasta el primer éxito.

Se asume que:

\[ Y \mid p \sim \text{Geométrica}(p) \]

donde \(p\) es la probabilidad de éxito en cada ensayo independiente.

Usando la parametrización:

\[ P(Y = y \mid p) = (1-p)^{y-1} p, \quad y = 1,2,3,\dots \]

2. Distribución a priori

En el enfoque bayesiano se asigna un prior para \(p\).

El prior conjugado de la Geométrica es la distribución Beta:

\[ p \sim \text{Beta}(\alpha,\beta) \]

con densidad:

\[ f(p) = \frac{\Gamma(\alpha+\beta)} {\Gamma(\alpha)\Gamma(\beta)} p^{\alpha-1}(1-p)^{\beta-1} \]

donde:

  • \(\alpha > 0\)
  • \(\beta > 0\)

3. Distribución Posterior

Aplicando el Teorema de Bayes:

\[ f(p \mid y) \propto f(y \mid p) f(p) \]

Sustituyendo:

\[ f(p \mid y) \propto (1-p)^{y-1} p \cdot p^{\alpha-1}(1-p)^{\beta-1} \]

Agrupando términos:

\[ f(p \mid y) \propto p^{\alpha}(1-p)^{\beta+y-2} \]

Por lo tanto:

\[ p \mid y \sim \text{Beta}(\alpha+1,\; \beta+y-1) \]

Caso general (muestra de tamaño n)

Si:

\[ Y_1,\dots,Y_n \mid p \sim \text{Geométrica}(p) \]

y definimos:

\[ S = \sum_{i=1}^n Y_i \]

La función de verosimilitud es:

\[ L(p) \propto p^n (1-p)^{S-n} \]

Entonces la posterior es:

\[ p \mid \mathbf{y} \sim \text{Beta}(\alpha + n,\; \beta + S - n) \]

4. Momentos de la Distribución Beta

Si:

\[ p \sim \text{Beta}(\alpha,\beta) \]

Esperanza

\[ E(p) = \frac{\alpha}{\alpha+\beta} \]

Varianza

\[ \text{Var}(p) = \frac{\alpha\beta} {(\alpha+\beta)^2(\alpha+\beta+1)} \]

Desviación estándar

\[ \text{DE}(p) = \sqrt{ \frac{\alpha\beta} {(\alpha+\beta)^2(\alpha+\beta+1)} } \]

5. Distribución Predictiva

La distribución marginal de \(Y\) (integrando \(p\)) se obtiene como:

\[ P(Y=y) = \int_0^1 (1-p)^{y-1} p \; f(p) \, dp \]

El resultado es:

\[ P(Y=y) = \frac{B(\alpha+1,\beta+y-1)} {B(\alpha,\beta)} \]

donde \(B(\cdot,\cdot)\) es la función Beta.

Esta distribución es conocida como la distribución Beta–Geométrica.

6. Interpretación Bayesiana

  • La distribución Beta modela la incertidumbre sobre la probabilidad de éxito \(p\).
  • La distribución Geométrica modela el número de intentos hasta el primer éxito.
  • La posterior actualiza la creencia sobre \(p\) después de observar datos.
  • La distribución predictiva incorpora incertidumbre paramétrica.

7. Resumen del Modelo

Si:

\[ Y_i \sim \text{Geométrica}(p) \]

y

\[ p \sim \text{Beta}(\alpha,\beta) \]

entonces:

\[ p \mid \mathbf{y} \sim \text{Beta}(\alpha+n,\; \beta+S-n) \]

Este modelo se conoce como el Modelo Beta–Geométrico, un caso clásico de conjugación bayesiana para datos de tiempo hasta el primer éxito.

# ==========================================
# MODELO BETA - GEOMÉTRICO
# Aplicación: intentos hasta primera compra
# ==========================================

library(ggplot2)
library(dplyr)
library(tidyr)

set.seed(123)

# -----------------------------
# 1. Datos simulados
# -----------------------------
p_real <- 0.30

y <- rgeom(50, prob = p_real) + 1  # Geométrica desde 1
n <- length(y)
S <- sum(y)

# -----------------------------
# 2. Prior
# -----------------------------
alpha_prior <- 2
beta_prior  <- 2

# -----------------------------
# 3. Posterior
# -----------------------------
alpha_post <- alpha_prior + n
beta_post  <- beta_prior + S - n

# -----------------------------
# 4. Momentos posteriores
# -----------------------------
mean_post <- alpha_post / (alpha_post + beta_post)

var_post <- (alpha_post * beta_post) /
  ((alpha_post + beta_post)^2 *
   (alpha_post + beta_post + 1))

sd_post <- sqrt(var_post)

# Intervalo creíble 95%
ic_post <- qbeta(c(0.025, 0.975),
                 alpha_post,
                 beta_post)

cat("===== RESULTADOS POSTERIORES =====\n")
## ===== RESULTADOS POSTERIORES =====
cat("Media:", mean_post, "\n")
## Media: 0.3095238
cat("Varianza:", var_post, "\n")
## Varianza: 0.001264608
cat("Desv.Est:", sd_post, "\n")
## Desv.Est: 0.03556133
cat("IC 95%:", ic_post, "\n")
## IC 95%: 0.2420855 0.3812577
# ==========================================
# 5. GRÁFICA PRIOR vs POSTERIOR
# ==========================================

p_seq <- seq(0, 1, length.out = 1000)

df <- data.frame(
  p = p_seq,
  Prior = dbeta(p_seq, alpha_prior, beta_prior),
  Posterior = dbeta(p_seq, alpha_post, beta_post)
)

df_long <- df %>%
  pivot_longer(-p,
               names_to = "Distribucion",
               values_to = "Densidad")

ggplot(df_long,
       aes(x = p,
           y = Densidad,
           color = Distribucion,
           fill = Distribucion)) +
  geom_line(size = 1.2) +
  geom_area(alpha = 0.2, position = "identity") +
  geom_vline(xintercept = mean_post,
             linetype = "dashed",
             size = 1) +
  labs(title = "Modelo Beta–Geométrico",
       subtitle = "Distribución Prior vs Posterior de p",
       x = "Probabilidad de conversión (p)",
       y = "Densidad") +
  theme_minimal(base_size = 15) +
  scale_color_manual(values = c("#2C7BB6", "#D7191C")) +
  scale_fill_manual(values = c("#2C7BB6", "#D7191C")) +
  theme(legend.position = "top")

# ==========================================
# 6. DISTRIBUCIÓN PREDICTIVA
# ==========================================

y_vals <- 1:15

pred_prob <- sapply(y_vals, function(k){
  beta(alpha_post + 1,
       beta_post + k - 1) /
    beta(alpha_post,
         beta_post)
})

pred_df <- data.frame(
  Intentos = y_vals,
  Probabilidad = pred_prob
)

ggplot(pred_df,
       aes(x = Intentos,
           y = Probabilidad)) +
  geom_col(fill = "#FDAE61",
           color = "black",
           width = 0.7) +
  labs(title = "Distribución Predictiva Beta–Geométrica",
       subtitle = "Probabilidad de número de intentos futuros",
       x = "Intentos hasta primera compra",
       y = "Probabilidad") +
  theme_minimal(base_size = 15)

Modelo Binomial Negativo

1. Definición

La distribución Binomial Negativa modela el número de fracasos antes de alcanzar un número fijo de éxitos \(r\) en una secuencia de ensayos Bernoulli independientes con probabilidad de éxito \(p\).

Se denota como:

\[ Y \sim \text{Binomial Negativa}(r, p) \]

2. Función de Probabilidad

Bajo la parametrización “número de fracasos antes de r éxitos”:

\[ P(Y = y) = \binom{y + r - 1}{r - 1} (1-p)^y p^r, \quad y = 0,1,2,\dots \]

donde:

  • \(r > 0\) es el número fijo de éxitos
  • \(0 < p < 1\) es la probabilidad de éxito

3. Interpretación

La variable aleatoria \(Y\) representa:

El número de fracasos que ocurren antes de observar \(r\) éxitos.

Caso especial: Si \(r = 1\), la distribución se reduce a la Geométrica.

4. Momentos

Para:

\[ Y \sim \text{Binomial Negativa}(r,p) \]

Esperanza

\[ E(Y) = \frac{r(1-p)}{p} \]

Varianza

\[ \text{Var}(Y) = \frac{r(1-p)}{p^2} \]

Desviación estándar

\[ \text{DE}(Y) = \sqrt{\frac{r(1-p)}{p^2}} \]

5. Propiedad de Sobredispersión

La distribución Binomial Negativa presenta:

\[ \text{Var}(Y) > E(Y) \]

lo que la hace útil para modelar datos de conteo con sobredispersión, donde la varianza excede la media.

6. Parametrización alternativa (media-dispersión)

Frecuentemente se parametriza en términos de:

\[ \mu = E(Y) \]

y un parámetro de dispersión \(k\), donde:

\[ \text{Var}(Y) = \mu + \frac{\mu^2}{k} \]

Esta forma es común en modelos de regresión para datos de conteo.

7. Relación con el Modelo Gamma–Poisson

Si:

\[ Y \mid \lambda \sim \text{Poisson}(\lambda) \]

y

\[ \lambda \sim \text{Gamma}(r, \beta) \]

entonces la distribución marginal de \(Y\) es:

\[ Y \sim \text{Binomial Negativa} \]

Por lo tanto, la Binomial Negativa surge como una mezcla Gamma–Poisson y modela heterogeneidad no observada en la tasa \(\lambda\).

8. Aplicaciones

La distribución Binomial Negativa se utiliza en:

  • Modelado de accidentes
  • Epidemiología
  • Conteo de reclamos en seguros
  • Análisis de demanda
  • Modelos de regresión para datos de conteo con sobredispersión

9. Resumen

Si:

\[ Y \sim \text{Binomial Negativa}(r,p) \]

entonces:

\[ E(Y) = \frac{r(1-p)}{p} \]

\[ \text{Var}(Y) = \frac{r(1-p)}{p^2} \]

La distribución es apropiada cuando los datos de conteo presentan mayor variabilidad que la distribución Poisson.

# ==========================================
# MODELO BINOMIAL NEGATIVO - APLICACIÓN
# Ejemplo: Número de reclamos diarios
# ==========================================

library(ggplot2)
library(MASS)
## 
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
set.seed(123)

# -----------------------------
# 1. Simulación de datos
# -----------------------------
r <- 5        # número de éxitos (parámetro tamaño)
p <- 0.4      # probabilidad de éxito

datos <- rnbinom(200, size = r, prob = p)

# -----------------------------
# 2. Estadísticos descriptivos
# -----------------------------
media <- mean(datos)
varianza <- var(datos)

cat("Media:", media, "\n")
## Media: 7.295
cat("Varianza:", varianza, "\n")
## Varianza: 17.61606
# Verificamos sobredispersión
cat("¿Varianza > Media?:", varianza > media, "\n")
## ¿Varianza > Media?: TRUE
# -----------------------------
# 3. Gráfico de frecuencias
# -----------------------------
ggplot(data.frame(x = datos), aes(x = x)) +
  geom_histogram(aes(y = ..density..),
                 bins = 20,
                 fill = "#FDAE61",
                 color = "black",
                 alpha = 0.7) +
  stat_function(fun = dnbinom,
                args = list(size = r, prob = p),
                color = "#D7191C",
                size = 1.2) +
  labs(title = "Distribución Binomial Negativa",
       subtitle = "Histograma + Densidad teórica",
       x = "Número de reclamos",
       y = "Densidad") +
  theme_minimal(base_size = 15)
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.720000

# -----------------------------
# 4. Comparación con Poisson
# -----------------------------
lambda_poisson <- media

ggplot(data.frame(x = datos), aes(x = x)) +
  geom_histogram(aes(y = ..density..),
                 bins = 20,
                 fill = "#2C7BB6",
                 color = "black",
                 alpha = 0.6) +
  stat_function(fun = dpois,
                args = list(lambda = lambda_poisson),
                color = "darkgreen",
                size = 1.2) +
  stat_function(fun = dnbinom,
                args = list(size = r, prob = p),
                color = "red",
                size = 1.2) +
  labs(title = "Comparación: Poisson vs Binomial Negativa",
       subtitle = "Rojo = Binomial Negativa | Verde = Poisson",
       x = "Número de reclamos",
       y = "Densidad") +
  theme_minimal(base_size = 15)
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 0.280000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 0.560000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 0.840000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 1.120000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 1.400000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 1.680000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 1.960000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 2.240000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 2.520000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 2.800000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 3.080000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 3.360000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 3.640000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 3.920000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 4.200000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 4.480000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 4.760000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 5.040000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 5.320000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 5.600000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 5.880000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 6.160000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 6.440000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 6.720000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 7.280000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 7.560000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 7.840000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 8.120000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 8.400000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 8.680000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 8.960000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 9.240000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 9.520000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 9.800000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 10.080000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 10.360000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 10.640000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 10.920000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 11.200000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 11.480000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 11.760000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 12.040000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 12.320000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 12.600000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 12.880000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 13.160000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 13.440000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 13.720000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 14.280000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 14.560000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 14.840000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 15.120000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 15.400000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 15.680000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 15.960000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 16.240000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 16.520000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 16.800000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 17.080000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 17.360000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 17.640000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 17.920000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 18.200000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 18.480000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 18.760000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 19.040000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 19.320000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 19.600000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 19.880000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 20.160000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 20.440000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 20.720000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 21.280000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 21.560000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 21.840000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 22.120000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 22.400000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 22.680000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 22.960000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 23.240000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 23.520000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 23.800000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 24.080000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 24.360000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 24.640000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 24.920000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 25.200000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 25.480000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 25.760000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 26.040000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 26.320000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 26.600000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 26.880000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 27.160000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 27.440000
## Warning in fun(x_trans, lambda = 7.295): non-integer x = 27.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 0.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 1.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 2.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 3.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 4.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 5.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 6.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 7.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 8.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 9.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 10.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 11.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 12.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 13.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 14.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 15.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 16.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 17.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 18.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 19.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 20.720000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.280000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.560000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 21.840000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.120000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.400000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.680000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 22.960000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.240000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.520000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 23.800000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.080000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.360000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.640000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 24.920000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.200000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.480000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 25.760000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.040000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.320000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.600000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 26.880000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.160000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.440000
## Warning in fun(x_trans, size = 5, prob = 0.4): non-integer x = 27.720000

# -----------------------------
# 5. Ajuste del modelo
# -----------------------------
modelo_nb <- fitdistr(datos, "Negative Binomial")

cat("\n===== Parámetros estimados =====\n")
## 
## ===== Parámetros estimados =====
print(modelo_nb)
##      size         mu    
##   5.2371043   7.2950699 
##  (0.9259280) (0.2954409)

Modelo Uniforme–Pareto

1. Planteamiento del Modelo

Sea \(X_1, \dots, X_n\) una muestra aleatoria tal que:

\[ X_i \mid \theta \sim \text{Uniforme}(0,\theta) \]

donde \(\theta > 0\) es un parámetro desconocido que representa el límite superior.

La función de densidad es:

\[ f(x \mid \theta) = \begin{cases} \frac{1}{\theta}, & 0 < x < \theta \\ 0, & \text{en otro caso} \end{cases} \]

2. Función de Verosimilitud

La función de verosimilitud para una muestra es:

\[ L(\theta \mid \mathbf{x}) = \prod_{i=1}^n \frac{1}{\theta} \mathbb{I}(0 < x_i < \theta) \]

Esto implica que:

\[ L(\theta \mid \mathbf{x}) = \frac{1}{\theta^n} \mathbb{I}(\theta > x_{(n)}) \]

donde:

\[ x_{(n)} = \max\{x_1,\dots,x_n\} \]

3. Distribución a priori

Se asigna como prior una distribución Pareto:

\[ \theta \sim \text{Pareto}(\alpha, \beta) \]

con densidad:

\[ \pi(\theta) = \alpha \beta^\alpha \theta^{-(\alpha+1)} \quad \text{para } \theta > \beta \]

donde:

  • \(\alpha > 0\) es el parámetro de forma
  • \(\beta > 0\) es el parámetro mínimo

4. Distribución Posterior

Aplicando el Teorema de Bayes:

\[ \pi(\theta \mid \mathbf{x}) \propto L(\theta \mid \mathbf{x}) \pi(\theta) \]

Sustituyendo:

\[ \pi(\theta \mid \mathbf{x}) \propto \theta^{-n} \theta^{-(\alpha+1)} \]

\[ \pi(\theta \mid \mathbf{x}) \propto \theta^{-(\alpha+n+1)} \]

con soporte:

\[ \theta > \max(\beta, x_{(n)}) \]

Por lo tanto:

\[ \theta \mid \mathbf{x} \sim \text{Pareto}(\alpha + n,\; \max(\beta, x_{(n)})) \]

Esto muestra que la distribución Pareto es conjugada para el modelo Uniforme(0,θ).

5. Momentos de la Distribución Pareto

Si:

\[ \theta \sim \text{Pareto}(\alpha,\beta) \]

Esperanza (si \(\alpha > 1\))

\[ E(\theta) = \frac{\alpha \beta}{\alpha - 1} \]

Varianza (si \(\alpha > 2\))

\[ \text{Var}(\theta) = \frac{\alpha \beta^2} {(\alpha - 1)^2 (\alpha - 2)} \]

6. Interpretación Bayesiana

  • La distribución Uniforme modela observaciones acotadas superiormente.
  • La Pareto representa incertidumbre sobre el límite máximo.
  • La posterior sigue siendo Pareto (propiedad de conjugación).
  • El estimador bayesiano bajo pérdida cuadrática es la media posterior.

7. Resumen del Modelo

Si:

\[ X_i \sim \text{Uniforme}(0,\theta) \]

y

\[ \theta \sim \text{Pareto}(\alpha,\beta) \]

entonces:

\[ \theta \mid \mathbf{x} \sim \text{Pareto}(\alpha+n,\; \max(\beta, x_{(n)})) \]

El modelo Uniforme–Pareto es un ejemplo clásico de conjugación bayesiana para parámetros de escala con soporte restringido.

# ==========================================
# MODELO UNIFORME - PARETO BAYESIANO
# ==========================================

library(ggplot2)

set.seed(123)

# -----------------------------
# 1. Datos simulados
# -----------------------------
theta_real <- 10
n <- 20

datos <- runif(n, 0, theta_real)

x_max <- max(datos)

# -----------------------------
# 2. Prior Pareto
# -----------------------------
alpha_prior <- 3
beta_prior  <- 5   # mínimo posible para theta

# -----------------------------
# 3. Posterior
# -----------------------------
alpha_post <- alpha_prior + n
beta_post  <- max(beta_prior, x_max)

cat("Máximo observado:", x_max, "\n")
## Máximo observado: 9.568333
cat("Alpha posterior:", alpha_post, "\n")
## Alpha posterior: 23
cat("Beta posterior:", beta_post, "\n")
## Beta posterior: 9.568333
# -----------------------------
# 4. Momentos posteriores
# -----------------------------
# Esperanza (si alpha_post > 1)
mean_post <- (alpha_post * beta_post) / (alpha_post - 1)

# Intervalo creíble 95%
q_pareto <- function(p, alpha, beta){
  beta / (1 - p)^(1/alpha)
}

ic_post <- c(
  q_pareto(0.025, alpha_post, beta_post),
  q_pareto(0.975, alpha_post, beta_post)
)

cat("Media posterior:", mean_post, "\n")
## Media posterior: 10.00326
cat("IC 95%:", ic_post, "\n")
## IC 95%: 9.578872 11.23288
# -----------------------------
# 5. Gráfica Prior vs Posterior
# -----------------------------
dpareto <- function(x, alpha, beta){
  ifelse(x >= beta,
         alpha * beta^alpha / x^(alpha+1),
         0)
}

theta_seq <- seq(beta_prior, 20, length.out = 1000)

df <- data.frame(
  theta = theta_seq,
  Prior = dpareto(theta_seq, alpha_prior, beta_prior),
  Posterior = dpareto(theta_seq, alpha_post, beta_post)
)

library(tidyr)

df_long <- pivot_longer(df,
                        -theta,
                        names_to = "Distribucion",
                        values_to = "Densidad")

ggplot(df_long,
       aes(x = theta,
           y = Densidad,
           color = Distribucion,
           fill = Distribucion)) +
  geom_line(size = 1.2) +
  geom_area(alpha = 0.2, position = "identity") +
  geom_vline(xintercept = x_max,
             linetype = "dashed",
             size = 1) +
  labs(title = "Modelo Uniforme–Pareto Bayesiano",
       subtitle = "Prior vs Posterior de θ",
       x = expression(theta),
       y = "Densidad") +
  theme_minimal(base_size = 15) +
  scale_color_manual(values = c("#2C7BB6", "#D7191C")) +
  scale_fill_manual(values = c("#2C7BB6", "#D7191C")) +
  theme(legend.position = "top")

Prior de Jeffreys (No Informativa)

1. Motivación

En estadística bayesiana, cuando no se tiene información previa sobre un parámetro \(\theta\), se busca una prior no informativa que represente neutralidad.
Jeffreys propuso una prior basada en la información de Fisher, que es invariante frente a reparametrizaciones.

2. Información de Fisher

La información de Fisher para un parámetro \(\theta\) es:

\[ I(\theta) = \mathbb{E} \left[ \left( \frac{\partial}{\partial \theta} \log f(Y \mid \theta) \right)^2 \Big| \theta \right] \]

donde \(f(Y \mid \theta)\) es la función de verosimilitud.

  • \(I(\theta)\) mide cuánta información sobre \(\theta\) contienen los datos.

3. Definición de la prior de Jeffreys

Jeffreys propone la prior no informativa:

\[ \pi_J(\theta) \propto \sqrt{I(\theta)} \]

Propiedades clave

  1. Invarianza: Si \(\phi = g(\theta)\), entonces

\[ \pi_J(\phi) \propto \sqrt{I(\phi)} \]

  1. No depende de la parametrización específica.
  2. Representa información mínima sobre \(\theta\).

4. Ejemplo 1: Modelo Bernoulli

Sea \(Y \sim \text{Bernoulli}(p)\), \(0 < p < 1\).
La verosimilitud es:

\[ f(y \mid p) = p^y (1-p)^{1-y}, \quad y=0,1 \]

Paso 1: Derivada del log-likelihood

\[ \log f(y \mid p) = y \log p + (1-y) \log (1-p) \]

\[ \frac{\partial}{\partial p} \log f(y \mid p) = \frac{y}{p} - \frac{1-y}{1-p} \]

Paso 2: Información de Fisher

\[ I(p) = \mathbb{E}\left[ \left( \frac{\partial}{\partial p} \log f(Y \mid p) \right)^2 \right] = \frac{1}{p(1-p)} \]

Paso 3: Prior de Jeffreys

\[ \pi_J(p) \propto \sqrt{I(p)} = \frac{1}{\sqrt{p(1-p)}} \]

Esto es equivalente a una Beta(1/2,1/2):

\[ p \sim \text{Beta}\left(\frac{1}{2},\frac{1}{2}\right) \]

5. Ejemplo 2: Modelo Normal con media desconocida y varianza conocida

Sea \(Y \sim N(\mu, \sigma^2)\), con \(\sigma^2\) conocida.

\[ \log f(y \mid \mu) = -\frac{1}{2\sigma^2} (y-\mu)^2 + \text{constante} \]

\[ \frac{\partial}{\partial \mu} \log f(y \mid \mu) = \frac{y-\mu}{\sigma^2} \]

\[ I(\mu) = \mathbb{E}\left[ \left( \frac{y-\mu}{\sigma^2} \right)^2 \right] = \frac{1}{\sigma^2} \]

\[ \pi_J(\mu) \propto \sqrt{I(\mu)} = \text{constante} \]

La prior de Jeffreys para \(\mu\) es uniforme: \(\pi_J(\mu) \propto 1\), reflejando completa no información.

6. Interpretación

  • La prior de Jeffreys no agrega información artificial.
  • Es especialmente útil en problemas de estimación de parámetros escalares o invariantes.
  • Garantiza que la posterior sea consistente bajo transformaciones de parámetros.

7. Resumen

  1. Información de Fisher:
    \[ I(\theta) = \mathbb{E}\Big[ (\partial_\theta \log f(Y \mid \theta))^2 \Big] \]

  2. Prior de Jeffreys:
    \[ \pi_J(\theta) \propto \sqrt{I(\theta)} \]

  3. Ejemplos:

    • Bernoulli: \(\pi_J(p) = \text{Beta}(1/2,1/2)\)
    • Normal (\(\sigma^2\) conocida): \(\pi_J(\mu) \propto 1\)
# ==========================================
# COMPARACIÓN BINOMIAL vs POISSON
# ==========================================

library(ggplot2)
library(dplyr)

set.seed(123)

# -----------------------------
# 1. DATOS SIMULADOS
# -----------------------------

# Binomial: número de éxitos en n ensayos
n_ensayos <- 10
p_real <- 0.4
datos_binom <- rbinom(50, size = n_ensayos, prob = p_real)

# Poisson: número de eventos por unidad
lambda_real <- 4
datos_pois <- rpois(50, lambda = lambda_real)

# -----------------------------
# 2. ESTADÍSTICOS DESCRIPTIVOS
# -----------------------------

# Binomial
media_binom <- mean(datos_binom)
var_binom <- var(datos_binom)

# Poisson
media_pois <- mean(datos_pois)
var_pois <- var(datos_pois)

cat("===== ESTADÍSTICOS =====\n")
## ===== ESTADÍSTICOS =====
cat("Binomial: media =", media_binom, ", var =", var_binom, "\n")
## Binomial: media = 4.14 , var = 2.775918
cat("Poisson: media =", media_pois, ", var =", var_pois, "\n")
## Poisson: media = 3.92 , var = 3.217959
# -----------------------------
# 3. INTERVALOS FREQUENTISTAS
# -----------------------------

# Binomial exacto (Clopper-Pearson)
binom_test <- binom.test(sum(datos_binom), n_ensayos*length(datos_binom))
cat("\n===== INTERVALO BINOMIAL =====\n")
## 
## ===== INTERVALO BINOMIAL =====
print(binom_test$conf.int)
## [1] 0.3704469 0.4585760
## attr(,"conf.level")
## [1] 0.95
# Poisson exacto
poisson_test <- poisson.test(sum(datos_pois))
cat("\n===== INTERVALO POISSON =====\n")
## 
## ===== INTERVALO POISSON =====
print(poisson_test$conf.int)
## [1] 169.5196 225.4440
## attr(,"conf.level")
## [1] 0.95
# -----------------------------
# 4. GRÁFICAS
# -----------------------------

# 4a: Histograma Binomial
df_binom <- data.frame(x = datos_binom)
ggplot(df_binom, aes(x = x)) +
  geom_histogram(aes(y = ..density..), bins = 10,
                 fill = "#2C7BB6", color = "black", alpha = 0.7) +
  stat_function(fun = function(x) dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos),
                color = "red", size = 1.2) +
  labs(title = "Distribución Binomial",
       subtitle = "Rojo = Ajuste teórico (media observada)",
       x = "Número de éxitos", y = "Densidad") +
  theme_minimal(base_size = 14)
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.070000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.140000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.210000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.280000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.350000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.420000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.490000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.560000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.630000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.700000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.770000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.840000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.910000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 1.980000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.050000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.120000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.190000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.260000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.330000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.400000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.470000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.540000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.610000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.680000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.750000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.820000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.890000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 2.960000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.030000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.100000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.170000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.240000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.310000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.380000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.450000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.520000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.590000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.660000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.730000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.800000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.870000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 3.940000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.010000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.080000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.150000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.220000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.290000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.360000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.430000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.500000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.570000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.640000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.710000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.780000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.850000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.920000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 4.990000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.060000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.130000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.200000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.270000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.340000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.410000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.480000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.550000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.620000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.690000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.760000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.830000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.900000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 5.970000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.040000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.110000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.180000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.250000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.320000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.390000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.460000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.530000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.600000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.670000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.740000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.810000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.880000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 6.950000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.020000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.090000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.160000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.230000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.300000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.370000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.440000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.510000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.580000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.650000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.720000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.790000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.860000
## Warning in dbinom(x, size = n_ensayos, prob = media_binom/n_ensayos):
## non-integer x = 7.930000

# 4b: Histograma Poisson
df_pois <- data.frame(x = datos_pois)
ggplot(df_pois, aes(x = x)) +
  geom_histogram(aes(y = ..density..), bins = 10,
                 fill = "#FDAE61", color = "black", alpha = 0.7) +
  stat_function(fun = function(x) dpois(x, lambda = media_pois),
                color = "darkgreen", size = 1.2) +
  labs(title = "Distribución Poisson",
       subtitle = "Verde = Ajuste teórico (media observada)",
       x = "Número de eventos", y = "Densidad") +
  theme_minimal(base_size = 14)
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.090000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.180000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.270000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.360000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.450000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.540000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.630000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.720000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.810000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.900000
## Warning in dpois(x, lambda = media_pois): non-integer x = 0.990000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.080000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.170000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.260000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.350000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.440000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.530000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.620000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.710000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.800000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.890000
## Warning in dpois(x, lambda = media_pois): non-integer x = 1.980000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.070000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.160000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.250000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.340000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.430000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.520000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.610000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.700000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.790000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.880000
## Warning in dpois(x, lambda = media_pois): non-integer x = 2.970000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.060000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.150000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.240000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.330000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.420000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.510000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.600000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.690000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.780000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.870000
## Warning in dpois(x, lambda = media_pois): non-integer x = 3.960000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.050000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.140000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.230000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.320000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.410000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.500000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.590000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.680000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.770000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.860000
## Warning in dpois(x, lambda = media_pois): non-integer x = 4.950000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.040000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.130000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.220000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.310000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.400000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.490000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.580000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.670000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.760000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.850000
## Warning in dpois(x, lambda = media_pois): non-integer x = 5.940000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.030000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.120000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.210000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.300000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.390000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.480000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.570000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.660000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.750000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.840000
## Warning in dpois(x, lambda = media_pois): non-integer x = 6.930000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.020000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.110000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.200000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.290000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.380000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.470000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.560000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.650000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.740000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.830000
## Warning in dpois(x, lambda = media_pois): non-integer x = 7.920000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.010000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.100000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.190000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.280000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.370000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.460000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.550000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.640000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.730000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.820000
## Warning in dpois(x, lambda = media_pois): non-integer x = 8.910000

# -----------------------------
# 5. COMPARACIÓN BINOMIAL vs POISSON
# -----------------------------

# Resumen en un solo gráfico
df_comp <- data.frame(
  Tipo = rep(c("Binomial", "Poisson"), each = 50),
  Valor = c(datos_binom, datos_pois)
)

ggplot(df_comp, aes(x = Valor, fill = Tipo)) +
  geom_histogram(aes(y = ..density..),
                 bins = 10, alpha = 0.6, position = "identity", color = "black") +
  labs(title = "Comparación Binomial vs Poisson",
       x = "Valor", y = "Densidad") +
  scale_fill_manual(values = c("#2C7BB6", "#FDAE61")) +
  theme_minimal(base_size = 14)

Teoría de la Decisión

1. Definición

La Teoría de la Decisión es un marco formal para elegir la mejor acción entre varias alternativas posibles, teniendo en cuenta:

  1. Las acciones posibles que un tomador de decisiones puede realizar.
  2. Los eventos o estados de la naturaleza que pueden ocurrir y que afectan los resultados.
  3. Los resultados o consecuencias asociados a cada combinación de acción y evento.
  4. Un criterio de preferencia o utilidad, que permite cuantificar la conveniencia de cada resultado.

Formalmente, un problema de decisión se define mediante un cuadro de decisión:

\[ \text{Decisión} = (A, S, u) \]

donde:

  • \(A = \{a_1, a_2, \dots, a_m\}\) es el conjunto de acciones posibles.
  • \(S = \{s_1, s_2, \dots, s_n\}\) es el conjunto de estados del mundo (eventos).
  • \(u(a_i, s_j)\) es la utilidad o pérdida asociada a elegir \(a_i\) cuando ocurre \(s_j\).

2. Matriz de Decisión

La información se organiza en una tabla de decisiones, donde:

Acción \(a_i\) Evento \(s_1\) Evento \(s_2\) … Evento \(s_n\)
\(a_1\) \(u(a_1,s_1)\) \(u(a_1,s_2)\) … \(u(a_1,s_n)\)
\(a_2\) \(u(a_2,s_1)\) \(u(a_2,s_2)\) … \(u(a_2,s_n)\)
… … … … …
\(a_m\) \(u(a_m,s_1)\) \(u(a_m,s_2)\) … \(u(a_m,s_n)\)
  • Cada celda representa el resultado de tomar acción \(a_i\) si ocurre el evento \(s_j\).
  • Los valores de utilidad pueden ser ganancias, pérdidas o cualquier medida cuantitativa de preferencia.

3. Principales Criterios de Decisión

3.1 Criterio de Maximax (optimista)

Selecciona la acción con máxima ganancia posible:

\[ a^* = \arg\max_{a_i} \left( \max_j u(a_i, s_j) \right) \]

3.2 Criterio Maximin (pesimista)

Selecciona la acción con mayor valor mínimo:

\[ a^* = \arg\max_{a_i} \left( \min_j u(a_i, s_j) \right) \]

3.3 Criterio de Laplace

Asume que todos los estados son igualmente probables:

\[ a^* = \arg\max_{a_i} \frac{1}{n} \sum_{j=1}^n u(a_i, s_j) \]

3.4 Criterio de Savage (minimización de arrepentimiento)

Se basa en la matriz de arrepentimiento, calculando la pérdida respecto a la mejor acción en cada estado:

\[ R(a_i, s_j) = \max_k u(a_k, s_j) - u(a_i, s_j) \]

Luego se selecciona:

\[ a^* = \arg\min_{a_i} \left( \max_j R(a_i, s_j) \right) \]

4. Ejemplo

Supongamos un agricultor que puede sembrar trigo o maíz, y el clima puede ser lluvioso o seco. La ganancia esperada (en miles de $) se muestra en la siguiente tabla:

Acción Lluvioso Seco
Trigo 30 10
Maíz 20 25

Paso 1: Maximax (optimista)

  • Trigo: max(30,10) = 30
  • Maíz: max(20,25) = 25
  • Decisión: Trigo

Paso 2: Maximin (pesimista)

  • Trigo: min(30,10) = 10
  • Maíz: min(20,25) = 20
  • Decisión: Maíz

Paso 3: Laplace

  • Trigo: (30+10)/2 = 20
  • Maíz: (20+25)/2 = 22.5
  • Decisión: Maíz

Paso 4: Savage

  1. Matriz de arrepentimiento:
Acción Lluvioso Seco
Trigo 0 15
Maíz 10 0
  1. Máximo arrepentimiento por acción:

  • Trigo: max(0,15) = 15

  • Maíz: max(10,0) = 10

  • Decisión: Maíz

  • La Teoría de la Decisión formaliza la elección racional frente a incertidumbre.
  • Los diferentes criterios reflejan distintas actitudes frente al riesgo:
    • Optimista → Maximax
    • Pesimista → Maximin
    • Neutral → Laplace
    • Minimización de arrepentimiento → Savage
  • Las tablas de decisión son esenciales para organizar la información y aplicar los criterios.
  • En un contexto probabilístico, si conocemos las probabilidades de cada estado \(P(s_j)\), el criterio de decisión esperado es:

\[ a^* = \arg\max_{a_i} \sum_{j=1}^{n} u(a_i, s_j) P(s_j) \]

  • Esto combina la Teoría de la Decisión con el Enfoque Bayesiano.
# ==========================================
# TEORÍA DE LA DECISIÓN
# Cálculo de criterios: Maximax, Maximin, Laplace, Savage
# ==========================================

library(dplyr)
library(ggplot2)
library(tidyr)

# -----------------------------
# 1. Tabla de decisiones
# -----------------------------
# Filas = Acciones, Columnas = Estados del mundo
# Ejemplo: Agricultor y clima
decisiones <- matrix(c(
  30, 10,   # Trigo: Lluvioso, Seco
  20, 25    # Maíz: Lluvioso, Seco
), nrow = 2, byrow = TRUE)

rownames(decisiones) <- c("Trigo", "Maíz")
colnames(decisiones) <- c("Lluvioso", "Seco")

decisiones_df <- as.data.frame(decisiones)
decisiones_df$Accion <- rownames(decisiones_df)
decisiones_df
##       Lluvioso Seco Accion
## Trigo       30   10  Trigo
## Maíz        20   25   Maíz
# ==========================================
# TEORÍA DE LA DECISIÓN 
# ==========================================

library(dplyr)
library(tidyr)
library(ggplot2)

# -----------------------------
# 1. Tabla de decisiones
# -----------------------------
# Filas = Acciones, Columnas = Estados
decisiones <- matrix(c(
  30, 10,  # Trigo: Lluvioso, Seco
  20, 25   # Maíz: Lluvioso, Seco
), nrow = 2, byrow = TRUE)

rownames(decisiones) <- c("Trigo", "Maíz")
colnames(decisiones) <- c("Lluvioso", "Seco")

decisiones_df <- as.data.frame(decisiones)
decisiones_df$Accion <- rownames(decisiones_df)

# -----------------------------
# 2. Cálculo de criterios
# -----------------------------
resultados <- decisiones_df %>%
  rowwise() %>%
  mutate(
    Maximax = max(c_across(Lluvioso:Seco)),
    Maximin = min(c_across(Lluvioso:Seco)),
    Laplace = mean(c_across(Lluvioso:Seco))
  ) %>%
  ungroup()

# Criterio Savage (minimización de arrepentimiento)
# 1. Mejor utilidad por columna (evento)
mejor_por_evento <- apply(decisiones[,1:2], 2, max)

# 2. Matriz de arrepentimiento
arrepentimiento <- t(t(mejor_por_evento - decisiones[,1:2]))

# 3. Máximo arrepentimiento por acción
resultados$Savage <- apply(arrepentimiento, 1, max)

# -----------------------------
# 3. Imprimir decisiones óptimas
# -----------------------------
decision_maximax <- resultados$Accion[which.max(resultados$Maximax)]
decision_maximin <- resultados$Accion[which.max(resultados$Maximin)]
decision_laplace <- resultados$Accion[which.max(resultados$Laplace)]
decision_savage <- resultados$Accion[which.min(resultados$Savage)]

cat("Criterio Maximax (optimista):", decision_maximax, "\n")
## Criterio Maximax (optimista): Trigo
cat("Criterio Maximin (pesimista):", decision_maximin, "\n")
## Criterio Maximin (pesimista): Maíz
cat("Criterio Laplace (neutral):", decision_laplace, "\n")
## Criterio Laplace (neutral): Maíz
cat("Criterio Savage (minimizar arrepentimiento):", decision_savage, "\n")
## Criterio Savage (minimizar arrepentimiento): Maíz
# -----------------------------
# 4. Visualización
# -----------------------------
df_long <- decisiones_df %>%
  pivot_longer(cols = Lluvioso:Seco,
               names_to = "Evento",
               values_to = "Utilidad")

ggplot(df_long, aes(x = Evento, y = Utilidad, fill = Accion)) +
  geom_col(position = "dodge") +
  labs(title = "Tabla de Decisiones",
       subtitle = "Utilidades de cada acción según el estado del mundo",
       x = "Evento / Estado del mundo",
       y = "Utilidad") +
  scale_fill_manual(values = c("#2C7BB6", "#FDAE61")) +
  theme_minimal(base_size = 14)