Modelo Geométrico y Poisson

0. Carga de librerías

library(gt)
library(e1071)
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

library(ggplot2)

## Warning: Can't find generic `as.gtable` in package gtable to register S3 method.
## ℹ This message is only shown to developers using devtools.
## ℹ Do you need to update gtable to the latest version?

1.Carga de datos

setwd("/cloud/project/")
datos<-read.csv("DerramesEEUU.csv", header = TRUE, sep=";" , dec=",",na.strings ="-")

2. Selección de la variable aleatoria

# Filtramos directamente la data antes de extraer la columna
datos_filtrados <- datos %>% filter(AnioAccidente <= 2016)
zona_anio <- datos_filtrados$AnioAccidente

3. Tabla de Distribución de Frecuencia

library(dplyr)
library(knitr)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

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

TDF_anio <- datos_filtrados %>%
  count(AnioAccidente, name = "ni") %>%
  arrange(AnioAccidente) %>%  
  mutate(hi_exacto = ni / sum(ni))

# Convertimos el año a carácter para evitar conflictos en la tabla
TDF_anio$AnioAccidente <- as.character(TDF_anio$AnioAccidente)

# Crear la fila de Sumatoria (TOTAL)
Sumatoria_anio <- data.frame(
  AnioAccidente = "TOTAL",
  ni = sum(TDF_anio$ni),
  hi_exacto = sum(TDF_anio$hi_exacto) 
) %>%
  mutate(
    N = "", 
    hi_porc = sprintf("%.2f", round(hi_exacto * 100, 2)), 
    hi = sprintf("%.4f", round(hi_exacto, 4))             
  ) %>%
  select(N, AnioAccidente, ni, hi_porc, hi)

TDF_anio <- TDF_anio %>%
  mutate(
    N = as.character(row_number()), 
    hi_porc = sprintf("%.2f", round(hi_exacto * 100, 2)),
    hi = sprintf("%.4f", round(hi_exacto, 4))
  ) %>%
  select(N, AnioAccidente, ni, hi_porc, hi)

TDF_final_anio <- rbind(TDF_anio, Sumatoria_anio)
colnames(TDF_final_anio) <- c("N", "Año", "ni", "hi (%)", "hi")

# --- MOSTRAR RESULTADO CON ESTRUCTURA KABLEEXTRA ---
titulo_formal_anio <- "Tabla N°1 <br/> Distribución de frecuencias de accidentes según el Año [2010 - 2016]"

kable(TDF_final_anio, 
      align = 'c',
      row.names = FALSE, 
      col.names = c("N°", "Año del Accidente", "ni", "hi (%)", "hi")) %>% 
  kable_styling(full_width = FALSE, position = "center", 
                bootstrap_options = c("striped", "hover", "condensed", "bordered")) %>%
  add_header_above(c(" " = 3, "Frecuencia relativa" = 2), bold = TRUE, background = "#D5D8DC") %>%
  add_header_above(setNames(5, titulo_formal_anio), align = "center", escape = FALSE, bold = FALSE, background = "white") %>%
  row_spec(0, bold = TRUE) %>%
  row_spec(nrow(TDF_final_anio), bold = TRUE, background = "#f2f2f2")

## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

Tabla N°1 Distribución de frecuencias de accidentes según el Año [2010 - 2016]
			Frecuencia relativa
N°	Año del Accidente	ni	hi (%)	hi
1	2010	346	12.55	0.1255
2	2011	336	12.18	0.1218
3	2012	362	13.13	0.1313
4	2013	400	14.50	0.1450
5	2014	447	16.21	0.1621
6	2015	453	16.42	0.1642
7	2016	414	15.01	0.1501
	TOTAL	2758	100.00	1.0000

4. Gráfica de Distribución de Frecuencia

# 1. PREPARACIÓN DE DATOS PARA EL GRÁFICO
datos_grafico_anio <- TDF_final_anio %>%
  filter(Año != "TOTAL") %>%              # Excluir la fila del total
  mutate(ni = as.numeric(ni),             # Asegurar que 'ni' sea número
         Año = factor(Año))               # Convertir a factor para eje X discreto

# 2. GENERAR EL GRÁFICO
ggplot(datos_grafico_anio, aes(x = Año, y = ni)) +
  geom_bar(stat = "identity", fill = "skyblue", color = "black", width = 0.7) +
  scale_y_continuous(limits = c(0, max(datos_grafico_anio$ni) * 1.15)) +
  
  # Títulos y Etiquetas
  labs(
    title = "Gráfica Temporal: Distribución de Incidentes [2010 - 2016]",
    x = "Año del Accidente",
    y = "Cantidad de Incidentes"
  ) +
  
  # Estilo limpio
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text.x = element_text(angle = 0, hjust = 0.5, color = "black"), 
    axis.text.y = element_text(color = "black")
  )

5. Conjetura

A partir de la Gráfica Temporal: Distribución de Incidentes [2010 - 2016], se propone la siguiente conjetura metodológica para el análisis de los datos:

Primer período (2010 - 2011): Las dos primeras barras correspondientes a los años 2010 y 2011 se abordarán mediante un modelo geométrico, adecuado para capturar la dinámica inicial de ocurrencia de los incidentes en este intervalo específico.

Segundo período (2012 - 2016): Las barras restantes correspondientes al período de 2012 a 2016 se trabajarán utilizando un modelo de Poisson, permitiendo modelar de forma óptima la tasa de eventos discretos en este tramo donde se observa una tendencia de crecimiento y posterior estabilización.

6. Parámetros

6.1 Modelo Geométrico

# 1. Selección y preparación de datos (Agrupación 1 Temporal)
library(dplyr)

# Asumiendo que TDF_final_anio ya está creado (y filtramos la fila TOTAL si existe)
grupo_1_anios <- c("2010", "2011")
tdf_geom1 <- TDF_final_anio %>%
  filter(Año %in% grupo_1_anios) %>%
  mutate(ni = as.numeric(ni)) %>%
  arrange(Año)

# Definimos los niveles del experimento (0 para 2010, 1 para 2011)
tdf_geom1$x <- 0:1  

# 2. Calcular parámetros Geométricos (Media y p)
# Calculamos la media ponderada de la agrupación
media_geom <- sum(tdf_geom1$x * tdf_geom1$ni) / sum(tdf_geom1$ni)

# Para dgeom, la media es (1-p)/p. Despejando p obtenemos:
p_geom <- 1 / (media_geom + 1)

# 3. Calcular Distribución Geométrica
# dgeom calcula la probabilidad geométrica para cada x dado p
P_Geometrica_Raw <- dgeom(tdf_geom1$x, prob = p_geom)

# NORMALIZACIÓN: Dado que evaluamos un subgrupo truncado (x=0 y x=1),
# normalizamos para que el modelo sume 100% en esta agrupación local.
P_Geometrica_Norm <- P_Geometrica_Raw / sum(P_Geometrica_Raw)

# Probabilidad Observada (Realidad)
prob_observada <- tdf_geom1$ni / sum(tdf_geom1$ni)

# 4. Crear Data Frame de resultados
Resultados_Geometrica <- data.frame(
  Anio = tdf_geom1$Año,
  Media_Obs = round(media_geom, 4),
  P_Geom = round(p_geom, 4),
  Prob_Realidad = round(prob_observada, 4),
  Prob_Modelo = round(P_Geometrica_Norm, 4)
)

# Mostrar resultados tabulares
print(Resultados_Geometrica)

##   Anio Media_Obs P_Geom Prob_Realidad Prob_Modelo
## 1 2010    0.4927 0.6699        0.5073      0.7518
## 2 2011    0.4927 0.6699        0.4927      0.2482

6.2 Modelo Poisson

# 1. Selección y preparación de datos (Agrupación 2 Temporal)
library(dplyr)
library(tidyr)
library(ggplot2)

# Asumiendo que TDF_final_anio ya está creado
grupo_2_anios <- c("2012", "2013", "2014", "2015", "2016")
tdf_pois2 <- TDF_final_anio %>%
  filter(Año %in% grupo_2_anios) %>%
  mutate(ni = as.numeric(ni)) %>%
  arrange(Año)

# Definimos los niveles cronológicos (0 para 2012, hasta 4 para 2016)
tdf_pois2$x <- 0:4  

# 2. Calcular parámetros del Modelo de Poisson (Lambda)
# Lambda (λ) en Poisson representa el número medio de eventos en el intervalo
lambda_pois <- sum(tdf_pois2$x * tdf_pois2$ni) / sum(tdf_pois2$ni)

# 3. Calcular Distribución de Poisson
# dpois calcula la probabilidad de ocurrencia para cada x dado λ
P_Poisson_Raw <- dpois(tdf_pois2$x, lambda = lambda_pois)

# NORMALIZACIÓN: Ajustamos la probabilidad teórica al 100% de esta ventana de 5 años
P_Poisson_Norm <- P_Poisson_Raw / sum(P_Poisson_Raw)

# Probabilidad Observada (Realidad empírica)
prob_observada <- tdf_pois2$ni / sum(tdf_pois2$ni)

# 4. Crear Data Frame de resultados numéricos
Resultados_Poisson <- data.frame(
  Anio = tdf_pois2$Año,
  Media_Lambda = round(lambda_pois, 4),
  Prob_Realidad = round(prob_observada, 4),
  Prob_Modelo = round(P_Poisson_Norm, 4)
)

# Mostrar resultados tabulares
print(Resultados_Poisson)

##   Anio Media_Lambda Prob_Realidad Prob_Modelo
## 1 2012       2.0756        0.1744      0.1334
## 2 2013       2.0756        0.1927      0.2770
## 3 2014       2.0756        0.2153      0.2875
## 4 2015       2.0756        0.2182      0.1989
## 5 2016       2.0756        0.1994      0.1032

7. Sobreposición de la realidad con el modelo

7.1 Modelo Geométrico

# 3. Crear DataFrame asegurando exactamente 2 filas para el período 2010-2011
df_comparativo_geom <- data.frame(
  Anio = factor(grupo_1_anios, levels = grupo_1_anios),
  Probabilidad_Observada = as.numeric(head(prob_observada, 2)),
  Modelo_Geometrico = as.numeric(head(P_Geometrica_Norm, 2))
) %>%
  pivot_longer(
    cols = c(Probabilidad_Observada, Modelo_Geometrico), 
    names_to = "Tipo", 
    values_to = "Probabilidad"
  )

# 4. Generar la Gráfica
ggplot(df_comparativo_geom, aes(x = Anio, y = Probabilidad, fill = Tipo)) +
  geom_bar(stat = "identity", position = position_dodge(), color = "black", width = 0.7) +
  scale_fill_manual(
    values = c("Modelo_Geometrico" = "#1f78b4", "Probabilidad_Observada" = "#a6cee3"),
    labels = c("Modelo Geométrico", "Probabilidad Observada")
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(df_comparativo_geom$Probabilidad, na.rm = TRUE) * 1.2)) +
  labs(
    title = "Gráfica Temporal: Relación entre el modelo geométrico y la realidad",
    subtitle = paste("Agrupación 1 (2010 - 2011) | Parámetro: p =", round(p_geom, 4)),
    x = "Año del Accidente",
    y = "Probabilidad",
    fill = ""
  ) +
  theme_bw() + 
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = 0, hjust = 0.5, color = "black")
  )

7.2 Modelo Poisson

# 5. GRÁFICA COMPARATIVA: MODELO VS REALIDAD
# Preparar DataFrame para formato de barras agrupadas (long format)
df_comparativo_pois <- data.frame(
  Año = factor(tdf_pois2$Año, levels = grupo_2_anios),
  `Probabilidad Observada` = prob_observada,
  `Modelo de Poisson` = P_Poisson_Norm
) %>%
  pivot_longer(cols = c(`Probabilidad.Observada`, `Modelo.de.Poisson`), 
               names_to = "Tipo", 
               values_to = "Probabilidad")

# Generar la Gráfica
ggplot(df_comparativo_pois, aes(x = Año, y = Probabilidad, fill = Tipo)) +
  # Barras agrupadas (position_dodge)
  geom_bar(stat = "identity", position = position_dodge(), color = "black", width = 0.7) +
  
  # Paleta de colores consistente
  scale_fill_manual(values = c("Modelo.de.Poisson" = "#1f78b4", 
                               "Probabilidad.Observada" = "#a6cee3"),
                    labels = c("Modelo", "Realidad")) +
  
  # Configuración de ejes y etiquetas
  scale_y_continuous(expand = c(0, 0), limits = c(0, max(df_comparativo_pois$Probabilidad) * 1.2)) +
  labs(
    title = "Gráfica Temporal: Relación entre el Modelo de Poisson y la Realidad",
    # Mostramos el parámetro Lambda en el subtítulo
    subtitle = paste("Agrupación 2 (2012 - 2016) | Parámetro: λ =", round(lambda_pois, 4)),
    x = "Año del Accidente",
    y = "Probabilidad",
    fill = ""
  ) +
  
  # Estilo visual limpio
  theme_bw() + 
  theme(
    legend.position = "top",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    axis.text.x = element_text(angle = 0, hjust = 0.5, color="black")
  )

8. Test de Bondad

8.1 Test de Pearson (Modelo Geométrico)

# 1. Preparación de variables (Agrupación 1)
grupo_1_anios <- c("2010", "2011")
tdf_g1 <- TDF_final_anio[TDF_final_anio$Año %in% grupo_1_anios, ]
tdf_g1$ni <- as.numeric(tdf_g1$ni)
tdf_g1$x <- 0:1  # 2010 = 0, 2011 = 1

# 2. Frecuencia Observada (Fo1)
Total_N1 <- sum(tdf_g1$ni)
Fo1 <- tdf_g1$ni / Total_N1

# 3. Cálculo de parámetros del Modelo Geométrico
media_obs1 <- sum(tdf_g1$x * tdf_g1$ni) / Total_N1
p_geom <- 1 / (media_obs1 + 1)

# 4. Frecuencia Esperada (Fe1) - Modelo Geométrico
Fe1_raw <- dgeom(tdf_g1$x, prob = p_geom)
Fe1 <- Fe1_raw / sum(Fe1_raw) # Normalización

# 5. Cálculo de la Correlación
Correlacion1 <- cor(Fo1, Fe1) * 100

# --- RESULTADOS ---
cat("--- COMPARATIVA MODELO GEOMÉTRICO (2010-2011) ---\n")

## --- COMPARATIVA MODELO GEOMÉTRICO (2010-2011) ---

cat("Correlación de Pearson:", round(Correlacion1, 2), "%\n")

## Correlación de Pearson: 100 %

8.2 Test de Chi-cuadrado

# 6. Test de Chi-cuadrado de Pearson para Modelo Geométrico
# Calcular el estadístico Chi-cuadrado (x2)
x2_1 <- sum(((Fo1 - Fe1)^2) / Fe1)

# Calcular el Valor Crítico (vc)
# Confianza del 95% (0.95), 1 grado de libertad (2 años - 1)
vc1 <- qchisq(0.95, 1)

cat("Chi-Cuadrado Calculado:", round(x2_1, 4), "\n")

## Chi-Cuadrado Calculado: 0.3205

cat("Valor Crítico (Tabla):", round(vc1, 4), "\n")

## Valor Crítico (Tabla): 3.8415

if (x2_1 < vc1) {
  cat("ESTADO: APRUEBA (La diferencia no es significativa)\n")
} else {
  cat("ESTADO: NO APRUEBA\n")
}

## ESTADO: APRUEBA (La diferencia no es significativa)

8.3 Test de Pearson (Modelo Poisson)

# 1. Preparación de variables (Agrupación 2)
grupo_2_anios <- c("2012", "2013", "2014", "2015", "2016")
tdf_g2 <- TDF_final_anio[TDF_final_anio$Año %in% grupo_2_anios, ]
tdf_g2$ni <- as.numeric(tdf_g2$ni)
tdf_g2$x <- 0:4  

# 2. Frecuencia Observada (Fo2)
Total_N2 <- sum(tdf_g2$ni)
Fo2 <- tdf_g2$ni / Total_N2
lambda_opt <- 3.04

Fe2_raw <- dpois(tdf_g2$x, lambda = lambda_opt)
Fe2 <- Fe2_raw / sum(Fe2_raw) # Normalización al 100%

# 4. Cálculo de la Correlación (Pearson R)
Correlacion2 <- cor(Fo2, Fe2) * 100

# --- IMPRESIÓN DE RESULTADOS ---
cat("\n--- RESULTADOS PRUEBA DE BONDAD DE AJUSTE (POISSON 2012-2016) ---\n")

## 
## --- RESULTADOS PRUEBA DE BONDAD DE AJUSTE (POISSON 2012-2016) ---

cat("1. Correlación de Pearson:  ", round(Correlacion2, 2), "%\n")

## 1. Correlación de Pearson:   98.71 %

8.4 Test de Chi-cuadrado

# 5. Test de Chi-cuadrado
# Sumatoria de ((Fo - Fe)^2 / Fe)
x2_2 <- sum(((Fo2 - Fe2)^2) / Fe2)

# Valor Crítico (vc) al 95% de confianza con 4 grados de libertad (5 años - 1)
vc2 <- qchisq(0.95, df = 4) 

cat("2. Chi-Cuadrado Calculado:  ", round(x2_2, 5), "\n")

## 2. Chi-Cuadrado Calculado:   0.25067

cat("3. Valor Crítico (Tabla):   ", round(vc2, 5), "\n")

## 3. Valor Crítico (Tabla):    9.48773

decision2 <- ifelse(x2_2 < vc2, "APRUEBA (El error es menor al límite)", "NO APRUEBA")
cat("   -> TEST CHI-CUADRADO:    ESTADO", decision2, "\n")

##    -> TEST CHI-CUADRADO:    ESTADO APRUEBA (El error es menor al límite)

9. Cálculo de Probabilidades

PREGUNTA N 1: ¿Cuál es la probabilidad acumulada de que un accidente haya ocurrido durante los años más recientes de alta incidencia (2015 y 2016) dentro del segundo grupo?

# Sumamos las frecuencias de los años 2015 y 2016 y dividimos por el total del Grupo 2
frec_recientes <- sum(tdf_g2$ni[tdf_g2$Año %in% c("2015", "2016")])
prob_recientes <- frec_recientes / sum(tdf_g2$ni)

cat("Probabilidad de ocurrencia en periodo reciente (2015-2016):", round(prob_recientes, 4), "\n")

## Probabilidad de ocurrencia en periodo reciente (2015-2016): 0.4176

PREGUNTA N 2: ¿Cuál es el peso estadístico exclusivo del año con el mayor pico de derrames (2015) frente al total de incidentes ocurridos entre 2012 y 2016?

# Calculamos el porcentaje que representa el año 2015 sobre el total de la Agrupación 2
peso_pico_2015 <- (tdf_g2$ni[tdf_g2$Año == "2015"] / sum(tdf_g2$ni)) * 100

cat("Peso estadístico del pico de accidentes (2015):", round(peso_pico_2015, 2), "%\n")

## Peso estadístico del pico de accidentes (2015): 21.82 %

10. Intervalo de confianza

# ==============================================================================
# 10. INTERVALO DE CONFIANZA
# ==============================================================================

# Calculamos la media, desviación estándar y tamaño de muestra usando tu tabla de frecuencias (ni)
media_incidentes <- mean(TDF_anio$ni)
desviacion_incidentes <- sd(TDF_anio$ni)
n_incidentes <- length(TDF_anio$ni)

# Nivel de confianza del 95% (Z = 1.96)
error_incidentes <- 1.96 * (desviacion_incidentes / sqrt(n_incidentes))

# Límites del intervalo de confianza (redondeado a entero por ser conteo de incidentes)
limite_inferior_inc <- round(media_incidentes - error_incidentes, 0)
limite_superior_inc <- round(media_incidentes + error_incidentes, 0)

# Creamos la tabla con el formato de intervalo
tabla_intervalo_inc <- data.frame(
  Intervalo = paste0(
    "P [",
    format(limite_inferior_inc, big.mark = ","),
    " < \u03bc < ",
    format(limite_superior_inc, big.mark = ","),
    "] = 95%"
  )
)

# Presentación de la tabla con formato elegante (gt) alineado a la UCE
library(gt)
library(dplyr)

tabla_intervalo_inc %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 2*"),
    subtitle = md("**Intervalo de confianza de la cantidad de accidentes (2010-2016)**")
  ) %>%
  tab_source_note(
    source_note = md(
      "Autor : Grupo 1"
    )
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    table.border.top.style = "solid",
    table.border.bottom.style = "solid",
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    row.striping.include_table_body = TRUE,
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black",
    table.border.left.color = "black",
    table.border.left.style = "solid",
    table.border.left.width = px(1),
    table.border.right.color = "black",
    table.border.right.style = "solid",
    table.border.right.width = px(1)
  )

Tabla Nro. 2
Intervalo de confianza de la cantidad de accidentes (2010-2016)
Intervalo
P [359 < μ < 429] = 95%
Autor : Grupo 1

11. Conclusión

El comportamiento de los accidentes por año se explica con un modelo geométrico de parámetro \(p = 0.6699\). Podemos afirmar con un 95% de confianza que la media aritmética real de la cantidad de accidentes se encuentra entre [359 y 429] y una desviación estándar de \(\sigma = 49.37\).

El comportamiento de los accidentes por año se explica con un modelo de Poisson de parámetro \(\lambda = 2.0756\). Podemos afirmar con un 95% de confianza que la media aritmética real de la cantidad de accidentes se encuentra entre [359 y 429] y una desviación estándar de \(\sigma = 49.37\).