Emisiones de PM2.5 (EMIS_PM25) Inferencial

1 Carga de Librerías

library(dplyr)
library(gt)

2 Carga de Datos

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

3 Selección de Variable Aleatoria

EMIS_PM25 <- as.numeric(datos$EMIS_PM25)
EMIS_PM25 <- na.omit(EMIS_PM25)
n <- length(EMIS_PM25)

4 Tratamiento de Datos

## [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.

5 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: Luis Cruz")
  ) %>%
  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: Luis Cruz

6 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")

7 Conjetura del Modelo

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.

8 Parámetros

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

media_pm25

## [1] 25.60317

sd_pm25

## [1] 3.326145

9 Sobreposición de la Realidad con el Modelo

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)

10 Test de Bondad

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: Luis Cruz"))

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: Luis Cruz

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

11 Cálculo de Probabilidades

¿Cuál es la probabilidad de que una mina tenga emisiones de PM2.5 entre el primer y el tercer intervalo 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

12 Intervalo 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(Intervalo = sprintf("P [%.2f< μ <%.2f] = 95%%", li, ls))

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: Luis Cruz")) %>%
  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 N°3
Intervalo de confianza de la media poblacional (95%)
Intervalo
P [25.48< μ <25.72] = 95%
Autor: Luis Cruz

13 Conclusiones

El comportamiento de EMIS_PM25 se explica con un modelo Normal de parámetros media = 25.60 y desviación estándar = 3.33. Podemos afirmar con un 95% de confianza que la media aritmética real de EMIS_PM25 se encuentra entre 25.48 y 25.72 toneladas, y una desviación estándar de 3.33 toneladas.