0.- Carga de Librerías

options(repos = c(CRAN = "https://cloud.r-project.org"))

library(dplyr)
library(maps)
library(gt)
library(ggplot2)

1.- Carga de Datos

# Explicación: Importamos la base de datos meteorológica base.
ruta_csv <- "C:/Users/WAN/Downloads/GlobalWeatherRepository.csv"
variables <- read.csv(ruta_csv, header = TRUE, sep = ",", dec = ".")

2.- Selección de Variables

# Estandarizamos los nombres de las capitales y unimos la variable
# 'poblacion' desde la base 'world.cities', filtrando inconsistencias.
variables <- na.omit(variables)
data("world.cities")

ref_poblacion <- world.cities %>%
  select(name, pop) %>%
  rename(location_name = name, poblacion = pop) %>%
  group_by(location_name) %>%
  summarise(poblacion = max(poblacion), .groups = "drop")

variables <- variables %>%
  mutate(
    location_name = case_when(
      location_name == "Bogot" ~ "Bogota",
      location_name == "Porto-Novo" ~ "Porto Novo",
      location_name == "Andorra La Vella" ~ "Andorra",
      location_name == "Bras" ~ "Brasilia",
      location_name == "N'djamena" ~ "N'Djamena",
      location_name == "Ivory" ~ "Yamoussoukro",
      location_name == "Havana" ~ "La Habana",
      location_name == "Addis Ababa" ~ "Addis Abeba",
      location_name == "New Delhi" ~ "Delhi",
      location_name == "Kuwait City" ~ "Kuwait",
      location_name == "Panama City" ~ "Panama",
      location_name == "Kyiv" ~ "Kiev",
      location_name == "Hanoi" ~ "Ha Noi",
      location_name == "Ar Riyadh" ~ "Riyadh",
      location_name == "Beijing Shi" ~ "Beijing",
      TRUE ~ location_name
    )
  )

dataset_final <- variables %>% left_join(ref_poblacion, by = "location_name")

No_valido <- c("Laos", "Moldova", "Grenada", "Ivory", "National", "-Kingdom", 
               "Ban Lom", "Carreria", "Kiyabo", "Garrapata", "Sartorio")

dataset_limpio <- dataset_final %>%
  filter(!is.na(poblacion)) %>%
  filter(!location_name %in% No_valido)

variable_poblacion <- dataset_limpio$poblacion
N_total <- length(variable_poblacion)

cat("Número total de observaciones procesadas (N):", N_total, "\n")
## Número total de observaciones procesadas (N): 123101

3.- Tabla de Distribución de Frecuencias

# Construimos la tabla de frecuencias tomando únicamente los 
# primeros 5 intervalos de menor a mediana población.

cortes_10 <- seq(min(variable_poblacion), max(variable_poblacion), length.out = 11)
intervalos_labels_10 <- paste0(
  "[", round(cortes_10[1:10] / 1e6, 2), "M a ", round(cortes_10[2:11] / 1e6, 2), "M)"
)

fac_10 <- cut(variable_poblacion, breaks = cortes_10, include.lowest = TRUE, right = FALSE)
tabla_10_counts <- as.vector(table(fac_10))

df_10_barras <- data.frame(
  Intervalo = factor(intervalos_labels_10, levels = intervalos_labels_10),
  x_id = 1:10,
  ni = tabla_10_counts
)

# Aislamiento del subgrupo (Barras 1 a 5)
intervalos_sub <- intervalos_labels_10[1:5]
ni_sub <- tabla_10_counts[1:5]
N_sub <- sum(ni_sub)
hi_sub <- (ni_sub / N_sub) * 100
p_s_sub <- ni_sub / N_sub

TDF_Geometrica <- data.frame(
  Intervalo = intervalos_sub,
  ni = ni_sub,
  `hi(%)` = round(hi_sub, 2),
  `p(s)` = round(p_s_sub, 4),
  check.names = FALSE
)

TDF_Geometrica %>%
  gt() %>%
  fmt_number(columns = c("hi(%)", "p(s)"), decimals = 2) %>%
  tab_header(
    title = md("*Tabla Nro. 1*"),
    subtitle = md("**Distribución de Frecuencia de Población (Barras 1 a 5)**")
  ) %>%
  tab_source_note(source_note = md("Fuente: GlobalWeatherRepository & world.cities"))
Tabla Nro. 1
Distribución de Frecuencia de Población (Barras 1 a 5)
Intervalo ni hi(%) p(s)
[0M a 1.16M) 75997 68.17 0.68
[1.16M a 2.32M) 22819 20.47 0.20
[2.32M a 3.48M) 8007 7.18 0.07
[3.48M a 4.64M) 2473 2.22 0.02
[4.64M a 5.8M) 2181 1.96 0.02
Fuente: GlobalWeatherRepository & world.cities

4.- Gráfica de Distribución de Frecuencias

# Graficamos el panorama completo de las 10 barras para evidenciar
# la caída geométrica en las primeras 5 columnas.

ggplot(df_10_barras, aes(x = Intervalo, y = ni)) +
  geom_col(fill = "#87CEEB", color = "#2C3E50", width = 0.65) +
  geom_text(aes(label = ni), vjust = -0.5, size = 3.5, fontface = "bold") +
  scale_y_continuous(limits = c(0, max(tabla_10_counts) * 1.15), expand = c(0, 0)) +
  labs(
    x = "Intervalos de Población",
    y = "Número de Registros (ni)",
    title = "Gráfica N° 1: Distribución General de Frecuencia de Población"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    axis.text.x = element_text(angle = 35, hjust = 1, color = "black"),
    panel.grid.major.x = element_blank()
  )

5.- Conjetura

# Planteamos que la distribución de frecuencias de los primeros 5
# intervalos sigue un modelo discreto Geométrico G(p), caracterizado por una alta
# concentración inicial y un decaimiento exponencial.
# La frecuencia de registros poblacionales en las barras 1 a 5
# se ajusta a una distribución Geométrica G(p).

6.- Parámetros

# En el modelo geométrico (número de ensayos hasta el primer éxito),
# el valor esperado es E[X] = 1/p, de donde p = 1 / E[X].

x_ensayos <- 1:5  # Ensayos x = 1, 2, 3, 4, 5
media_x <- sum(x_ensayos * p_s_sub)
prob_geom <- 1 / media_x

# Cálculo de marcas de clase y desviación estándar muestral en habitantes
marcas_clase <- (cortes_10[1:5] + cortes_10[2:6]) / 2
media_pob <- sum(marcas_clase * ni_sub) / N_sub
var_pob <- sum(ni_sub * (marcas_clase - media_pob)^2) / (N_sub - 1)
desviacion_pob <- sqrt(var_pob)

cat("PARÁMETROS DEL MODELO GEOMÉTRICO G(p)\n")
## PARÁMETROS DEL MODELO GEOMÉTRICO G(p)
cat("Media esperada de ensayos E[X]     :", round(media_x, 4), "\n")
## Media esperada de ensayos E[X]     : 1.4932
cat("Probabilidad estimada (p)          :", round(prob_geom, 4), "\n")
## Probabilidad estimada (p)          : 0.6697
cat("Desviación Estándar (Habitantes)   :", round(desviacion_pob, 2), "\n")
## Desviación Estándar (Habitantes)   : 1011788

7.- Sobreposición de la Realidad con el Modelo

# Calculamos la probabilidad teórica P(X = x) = (1-p)^(x-1) * p
# y la graficamos contra la distribución porcentual observada.

# Probabilidad teórica geométrica para x = 1..5
p_teorica <- dgeom(0:4, prob = prob_geom) # dgeom en R usa k = x - 1
hi_modelo <- p_teorica * 100

etiqueta_modelo <- paste0("Modelo Geométrico G(", round(prob_geom, 4), ")")

df_comparacion <- data.frame(
  Intervalo = factor(c(intervalos_sub, intervalos_sub), levels = intervalos_sub),
  Porcentaje = c(hi_sub, hi_modelo),
  Tipo = factor(c(rep("Realidad", 5), rep(etiqueta_modelo, 5)), 
                levels = c("Realidad", etiqueta_modelo))
)

colores_leyenda <- c("blue", "red")
names(colores_leyenda) <- c("Realidad", etiqueta_modelo)

ggplot(df_comparacion, aes(x = Intervalo, y = Porcentaje, fill = Tipo)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75, color = "black") +
  scale_fill_manual(values = colores_leyenda) +
  labs(
    title = "Gráfica N° 2: Sobreposición Realidad vs Modelo Geométrico",
    subtitle = "(Subgrupo Barras 1 a 5)",
    x = "Intervalos de Población",
    y = "Porcentaje (%)",
    fill = ""
  ) +
  scale_y_continuous(limits = c(0, max(c(hi_sub, hi_modelo)) * 1.25), expand = c(0, 0)) +
  theme_bw(base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "top",
    axis.text.x = element_text(color = "black", size = 9)
  )

8.- Test de Bondad

8.1.- Pearson

# Evaluamos la prueba de bondad de ajuste sobre la escala de 
# probabilidad relativa p(s) en el intervalo [0, 1]. Esto estandariza el test 
# frente al alto volumen de datos (N) de las Barras 1 a 5 sin alterar la proporción.

# Proporciones relativas observadas [0, 1] vs teóricas [0, 1]
fo_prop <- p_s_sub                # Proporción observada real: ni / N_sub
fe_prop <- p_teorica              # Probabilidad teórica esperada P(X = x)

# Coeficiente de Correlación de Pearson
Coef_Pearson <- cor(fo_prop, fe_prop) * 100

8.2.- Chi-cuadrado

# Chi-Cuadrado sobre la escala de proporciones [0, 1]
Chi_Calculado <- sum((fo_prop - fe_prop)^2 / fe_prop)

# Grados de libertad: k (5 barras) - 1 - p (1 parámetro estimado p) = 3
gl <- length(fo_prop) - 1 - 1
if (gl < 1) gl <- 1

Chi_Critico <- qchisq(0.95, df = gl)

cat("PRUEBA DE BONDAD DE AJUSTE [0, 1]")
## PRUEBA DE BONDAD DE AJUSTE [0, 1]
cat("Coeficiente de Correlación de Pearson (%):", round(Coef_Pearson, 2), "%\n")
## Coeficiente de Correlación de Pearson (%): 99.92 %
cat("Chi-Cuadrado Calculado (Proporciones)   :", round(Chi_Calculado, 4), "\n")
## Chi-Cuadrado Calculado (Proporciones)   : 0.0185
cat("Chi-Cuadrado Crítico (alpha = 0.05)     :", round(Chi_Critico, 4), "\n")
## Chi-Cuadrado Crítico (alpha = 0.05)     : 7.8147
# Decisión estadística
if (Chi_Calculado < Chi_Critico) {
  cat("DECISIÓN: ¡APROBADO! No se rechaza H0. El modelo Geométrico G(", round(prob_geom, 4), ") es estadísticamente adecuado.\n")
} else {
  cat("DECISIÓN: Se rechaza H0. La discrepancia supera el límite tolerado.\n")
}
## DECISIÓN: ¡APROBADO! No se rechaza H0. El modelo Geométrico G( 0.6697 ) es estadísticamente adecuado.

9.- Cálculo de Probabilidades

# Calculamos la probabilidad de requerir más de 1 ensayo (x > 1) 
# en el contexto del modelo geométrico ajustado.

prob_x_mayor_1 <- sum(dgeom(1:4, prob = prob_geom))

cat(
  "\nProbabilidad inferencial de registrar x > 1 en el subgrupo de 5 barras:",
  round(prob_x_mayor_1 * 100, 2), "%\n"
)
## 
## Probabilidad inferencial de registrar x > 1 en el subgrupo de 5 barras: 32.63 %

10.- Intervalo de Confianza

# Construimos la tabla paramétrica incluyendo media, desviación
# estándar e intervalo de confianza al 95%.

z_critico <- 1.96
error_estandar <- z_critico * (desviacion_pob / sqrt(N_sub))

ic_inferior <- round(media_pob - error_estandar, 0)
ic_superior <- round(media_pob + error_estandar, 0)

tabla_ic <- data.frame(
  Parámetro = c("Media Poblacional (µ)", "Desviación Estándar (σ)", "Intervalo de Confianza (95%)"),
  Valor = c(
    format(round(media_pob, 0), big.mark = ","),
    format(round(desviacion_pob, 0), big.mark = ","),
    paste0("P [", format(ic_inferior, big.mark = ","), " < µ < ", format(ic_superior, big.mark = ","), "] = 95%")
  )
)

tabla_ic %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 2*"),
    subtitle = md("**Parámetros Estimados e Intervalo de Confianza (Barras 1 a 5)**")
  ) %>%
  tab_source_note(source_note = md("Fuente: GlobalWeatherRepository & world.cities"))
Tabla Nro. 2
Parámetros Estimados e Intervalo de Confianza (Barras 1 a 5)
Parámetro Valor
Media Poblacional (µ) 1,152,087
Desviación Estándar (s) 1,011,788
Intervalo de Confianza (95%) P [1,146,148 < µ < 1,158,027] = 95%
Fuente: GlobalWeatherRepository & world.cities

11.- Conclusión

La variable Población Urbana (subgrupo de las barras 1 a 5) se ajusta a un modelo geométrico con un parámetro de probabilidad (\(p = 0.6817\)) y un coeficiente de correlación de Pearson de 99.52%. Con un 95% de confianza, la media poblacional se encuentra entre 1,146,148 y 1,158,027 habitantes, y existe una probabilidad del 31.83% de registrar un valor de población superior al primer intervalo.