library(dplyr)
library(gt)

1 Carga de Datos y Librerías

datos <- read.csv("~/Estudio/TERCER SEMESTRE/Estadistica/Dataset.csv",
                   sep = ";", stringsAsFactors = FALSE)

EMIS_PM25 <- as.numeric(datos$EMIS_PM25)
EMIS_PM25 <- na.omit(EMIS_PM25)
n <- length(EMIS_PM25)
## [1] 2890

Nota: a diferencia de EMIS_CO2/EMIS_NOX/EMIS_CH4 (registradas a nivel estatal y repetidas en cada mina del mismo estado), EMIS_PM25 varía de forma única entre minas, por lo que se trabaja directamente con las n = 2890 observaciones a nivel-mina, sin reducir a nivel-estado.

2 Tabla de Distribución de Frecuencias

k <- ceiling(1 + 3.322 * log10(n))
k
## [1] 13
breaks <- seq(min(EMIS_PM25), max(EMIS_PM25) + 1e-6, length.out = k + 1)

histograma_pm25 <- hist(EMIS_PM25, breaks = breaks, plot = FALSE)

lis <- round(histograma_pm25$breaks[1:k], 2)
lss <- round(histograma_pm25$breaks[2:(k+1)], 2)
MC  <- round(histograma_pm25$mids, 2)
ni  <- histograma_pm25$counts
hi  <- round((ni/sum(ni))*100, 2)

TDFpm25 <- data.frame(lis, lss, MC, ni, hi)

fila_total <- data.frame(
  lis = "TOTAL", lss = "", MC = "",
  ni  = sum(TDFpm25$ni),
  hi  = round((sum(TDFpm25$ni) / n) * 100, 2)
)

TDFpm25_total <- rbind(TDFpm25, fila_total)

tabla_pm25 <- TDFpm25_total %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N°1**"),
    subtitle = md("Tabla de distribución de frecuencias de EMIS_PM25")
  ) %>%
  cols_label(
    lis = "Límite inferior",
    lss = "Límite superior",
    MC  = "Marca de clase",
    ni  = "Frecuencia absoluta",
    hi  = "Frecuencia relativa (%)"
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 4 - Minas")
  ) %>%
  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"
  )

tabla_pm25
Tabla N°1
Tabla de distribución de frecuencias de EMIS_PM25
Límite inferior Límite superior Marca de clase Frecuencia absoluta Frecuencia relativa (%)
16.4 17.82 17.11 22 0.76
17.82 19.24 18.53 59 2.04
19.24 20.66 19.95 144 4.98
20.66 22.09 21.38 207 7.16
22.09 23.51 22.8 353 12.21
23.51 24.93 24.22 416 14.39
24.93 26.35 25.64 499 17.27
26.35 27.77 27.06 440 15.22
27.77 29.19 28.48 346 11.97
29.19 30.62 29.9 210 7.27
30.62 32.04 31.33 116 4.01
32.04 33.46 32.75 53 1.83
33.46 34.88 34.17 25 0.87
TOTAL 2890 100.00
Autor: Grupo 4 - Minas

3 Gráfica de Distribución de Frecuencias

hist(EMIS_PM25,
     breaks = breaks,
     main = "Gráfica N°1: Distribución de EMIS_PM25",
     xlab = "Emisiones de PM2.5 (toneladas)",
     ylab = "Frecuencia",
     col  = "gray")

4 Conjetura del Modelo

media_pm25 <- mean(EMIS_PM25)
sd_pm25    <- sd(EMIS_PM25)

media_pm25
## [1] 25.60317
sd_pm25
## [1] 3.326145
densidad_hist <- hist(EMIS_PM25, plot = FALSE)$density
densidad_curva <- dnorm(seq(min(EMIS_PM25), max(EMIS_PM25), length.out = 500),
                         mean = media_pm25, sd = sd_pm25)

ylim_max <- max(densidad_hist, densidad_curva) * 1.05

hist(EMIS_PM25,
     freq = FALSE,
     ylim = c(0, ylim_max),
     main = "Gráfica N°2: Comparación de la realidad con el modelo Normal\nde EMIS_PM25",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de PM2.5 (toneladas)",
     col = "lightgray",
     border = "black")

x <- seq(min(EMIS_PM25), max(EMIS_PM25), length.out = 500)
curve(dnorm(x, mean = media_pm25, sd = sd_pm25),
      from = min(EMIS_PM25),
      to = max(EMIS_PM25),
      n = 500,
      col = "blue",
      lwd = 2,
      add = TRUE)

EMIS_PM25 presenta una distribución prácticamente simétrica en torno a su media , por lo que un modelo Normal es adecuado desde el planteamiento inicial, sin necesidad de recurrir a un modelo alternativo.

Las emisiones de PM2.5 son el resultado de múltiples procesos mineros y factores ambientales. Esta combinación genera una variabilidad natural que puede aproximarse mediante una distribución Normal

5 Test de Aprobación

cortes <- cut(EMIS_PM25, breaks = breaks, include.lowest = TRUE, right = FALSE)

Fo_1 <- as.vector(table(cortes))
Fo_1
##  [1]  22  59 144 207 353 416 499 440 346 210 116  53  25
Fe_1 <- n * diff(pnorm(breaks, mean = media_pm25, sd = sd_pm25))
Fe_1
##  [1]  19.72269  52.81343 118.12676 220.69616 344.42692 449.01916 488.99442
##  [8] 444.85252 338.06427 214.60880 113.80230  50.40763  18.64946
# Expresar Fo y Fe en porcentaje para comparar
Fo_pct <- (Fo_1/n)*100
Fe_pct <- (Fe_1/n)*100

plot(Fo_pct, Fe_pct,
     xlim = c(0, max(Fo_pct, Fe_pct)),
     ylim = c(0, max(Fo_pct, Fe_pct)),
     main = "Gráfica N°3: Correlación de frecuencias observadas y esperadas\ndel modelo Normal de EMIS_PM25",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     col  = "blue3",
     pch  = 19)

abline(a = 0, b = 1, col = "red", lwd = 2)

Correlacion_1 <- cor(Fo_pct, Fe_pct)*100
Correlacion_1
## [1] 99.70658
# Chi-cuadrado
grados_libertad_1 <- length(Fo_1) - 1 - 2   # 2 parametros estimados (media, sd)
grados_libertad_1
## [1] 10
nivel_significancia <- 0.95

x2_1 <- sum((Fo_1 - Fe_1)^2 / Fe_1)
x2_1
## [1] 13.02729
umbral_aceptacion_1 <- qchisq(nivel_significancia, grados_libertad_1)
umbral_aceptacion_1
## [1] 18.30704
# Tabla resumen
Variable <- c("Emisiones de PM2.5 (ton)")

tabla_resumen_1 <- data.frame(
  Variable,
  round(Correlacion_1, 2),
  round(x2_1, 2),
  round(umbral_aceptacion_1, 2)
)

colnames(tabla_resumen_1) <- c("Variable", "Test Pearson (%)", "Chi Cuadrado", "Umbral de aceptación")

tabla_resumen_1 %>%
  gt() %>%
  tab_header(title = md("**Tabla N°2**"),
             subtitle = md("Resumen del test de bondad de ajuste (modelo Normal)")) %>%
  tab_source_note(source_note = md("Autor: Grupo 4 - Minas"))
Tabla N°2
Resumen del test de bondad de ajuste (modelo Normal)
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Emisiones de PM2.5 (ton) 99.71 13.03 18.31
Autor: Grupo 4 - Minas

Con Chi² = 13.03 por debajo del umbral de aceptación 18.31 (95% de confianza, gl = 10), el modelo Normal se acepta.

6 Cálculo de Probabilidades

¿Cuál es la probabilidad de que una mina tenga emisiones de PM2.5 entre el primer y el tercer cuartil de la distribución observada?

q1 <- round(quantile(EMIS_PM25, 0.25), 2)
q3 <- round(quantile(EMIS_PM25, 0.75), 2)

Probabilidad_1 <- (pnorm(q3, media_pm25, sd_pm25) - pnorm(q1, media_pm25, sd_pm25)) * 100
Probabilidad_1
##      75% 
## 51.16856
x <- seq(min(EMIS_PM25), max(EMIS_PM25), length.out = 300)
plot(x,
     dnorm(x, media_pm25, sd_pm25),
     col = "skyblue3",
     lwd = 1,
     type = "l",
     main = "Gráfica N°4: Cálculo de probabilidades",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de PM2.5 (toneladas)")

x_section <- seq(q1, q3, length.out = 200)
y_section <- dnorm(x_section, media_pm25, sd_pm25)

lines(x_section, y_section, col = "red", lwd = 2)
polygon(c(x_section, rev(x_section)),
        c(y_section, rep(0, length(y_section))),
        col = rgb(1, 0, 0, 0.6))

legend("topright",
       legend = c("Modelo Normal", "Área de Probabilidad"),
       col = c("skyblue3", "red"),
       lwd = 2,
       cex = 0.7)

texto_prob <- paste0("Probabilidad = ", round(Probabilidad_1, 2), " %")
text(x = median(EMIS_PM25), y = max(dnorm(x, media_pm25, sd_pm25))*0.9,
     labels = texto_prob,
     col = "black",
     cex = 0.8,
     font = 2)

# De las n minas del dataset, cuantas tendrian emisiones en ese rango
cantidad_1 <- (pnorm(q3, media_pm25, sd_pm25) - pnorm(q1, media_pm25, sd_pm25)) * n
cantidad_1
##      75% 
## 1478.772

7 Intervalos de Confianza

x_media <- mean(EMIS_PM25)
x_media
## [1] 25.60317
sigma <- sd(EMIS_PM25)
sigma
## [1] 3.326145
n_ic <- length(EMIS_PM25)
n_ic
## [1] 2890
e <- (sigma/sqrt(n_ic)) * qt(0.975, n_ic - 1)
e
## [1] 0.1213172
li <- x_media - e
li
## [1] 25.48185
ls <- x_media + e
ls
## [1] 25.72448
tabla_media <- data.frame(
  round(li, 2),
  Variable,
  round(ls, 2),
  round(e, 2)
)

colnames(tabla_media) <- c("Límite inferior", "Media poblacional", "Límite superior", "Error estándar")

tabla_media %>%
  gt() %>%
  tab_header(title = md("**Tabla N°3**"),
             subtitle = md("Intervalo de confianza de la media poblacional (95%)")) %>%
  tab_source_note(source_note = md("Autor: Grupo 4 - Minas"))
Tabla N°3
Intervalo de confianza de la media poblacional (95%)
Límite inferior Media poblacional Límite superior Error estándar
25.48 Emisiones de PM2.5 (ton) 25.72 0.12
Autor: Grupo 4 - Minas

8 Conclusiones

La variable EMIS_PM25 no esta duplicada a nivel estatal, sino que varia de forma unica entre minas, por lo que se trabaja directamente con las n = 2890 observaciones disponibles. Dado que la distribucion es practicamente simetrica (asimetria cercana a 0, media y mediana casi identicas), se conjeturo directamente un modelo Normal(media = 25.60, sd = 3.33), el cual ajusta muy bien a los datos observados: una correlacion de Pearson de 99.71% y un Chi-cuadrado de 13.03, muy por debajo del umbral de aceptacion de 18.31 (95% de confianza). La probabilidad de que una mina emita entre 23.29 y 27.90 toneladas de PM2.5 es de 51.17%, equivalente a aproximadamente 1479 de las 2890 minas del dataset. Finalmente, la media poblacional de emisiones se estima, con 95% de confianza, entre 25.48 y 25.72 toneladas.