Emisiones de NOx (EMIS_NOX) 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

# La columna de NOx viene con formato inconsistente: algunos valores usan el
# punto como separador de miles (ej. "40.103.865") y otros como separador
# decimal (ej. "14859.23"). Se limpia contando los puntos: si hay mas de uno,
# son separadores de miles y se eliminan; si hay uno solo, es decimal y se
# conserva.
limpiar_numero <- function(x) {
  x <- trimws(as.character(x))
  x[x %in% c("N/A", "N//A", "NA", "")] <- NA
  n_puntos <- lengths(regmatches(x, gregexpr("\\.", x)))
  x <- ifelse(n_puntos > 1, gsub("\\.", "", x), x)
  as.numeric(x)
}

EMIS_NOX <- limpiar_numero(datos$NOx.total.emissions..tons..by.state)
EMIS_NOX <- na.omit(EMIS_NOX)
n <- length(EMIS_NOX)
n

## [1] 2993

4 Tabla de Distribución de Frecuencias

k <- ceiling(1 + 3.322 * log10(n))
k

## [1] 13

breaks <- seq(min(EMIS_NOX), max(EMIS_NOX) + 1e-6, length.out = k + 1)

histograma_nox <- hist(EMIS_NOX, breaks = breaks, plot = FALSE)

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

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

fila_total <- data.frame(
  lis = "TOTAL", lss = "", MC = "",
  ni  = sum(TDFnox$ni),
  hi  = (sum(TDFnox$hi))
)

TDFnox_total <- rbind(TDFnox, fila_total)

tabla_nox <- TDFnox_total %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N°1**"),
    subtitle = md("Tabla de distribución de frecuencias de EMIS_NOX (nivel-mina)")
  ) %>%
  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_nox

Tabla N°1
Tabla de distribución de frecuencias de EMIS_NOX (nivel-mina)
Límite inferior	Límite superior	Marca de clase	Frecuencia absoluta	Frecuencia relativa (%)
252.01	6857936.08	3429094.05	432	14.43
6857936.08	13715620.16	10286778.12	126	4.21
13715620.16	20573304.24	17144462.2	271	9.05
20573304.24	27430988.31	24002146.27	259	8.65
27430988.31	34288672.39	30859830.35	155	5.18
34288672.39	41146356.47	37717514.43	618	20.65
41146356.47	48004040.54	44575198.5	911	30.44
48004040.54	54861724.62	51432882.58	159	5.31
54861724.62	61719408.69	58290566.66	0	0.00
61719408.69	68577092.77	65148250.73	0	0.00
68577092.77	75434776.85	72005934.81	0	0.00
75434776.85	82292460.92	78863618.89	0	0.00
82292460.92	89150145	85721302.96	62	2.07
TOTAL			2993	99.99
Autor: Luis Cruz

5 Gráfica de Distribución de Frecuencias

hist(EMIS_NOX,
     breaks = breaks,
     main = "Gráfica N°1: Distribución de EMIS_NOX (nivel-mina)",
     xlab = "Emisiones de NOx (toneladas)",
     ylab = "Frecuencia",
     col  = "gray")

6 Tratamiento de Datos

Las emisiones de NOx se registran a nivel estatal (no por mina), y cada una de las 2993 minas del dataset hereda el valor de su estado. Esto significa que, aunque la tabla y el histograma anteriores usan todas las filas (igual que en la fase Descriptiva), la muestra real e independiente para conjeturar y probar un modelo de probabilidad son los valores únicos por estado. Usar n = 2993 en el test de aprobación infla artificialmente el tamaño muestral y sobre-potencia el Chi-cuadrado (mismo criterio ya documentado para EMIS_CO2).

Por esta razón, desde aquí en adelante se trabaja con EMIS_NOX_edo.

EMIS_NOX_edo <- unique(EMIS_NOX)
n_edo <- length(EMIS_NOX_edo)
n_edo

## [1] 47

7 Conjetura del Modelo

Se conjeturó primero un modelo Normal, siguiendo el mismo tratamiento usado en el resto de variables continuas del proyecto. Sin embargo, ese modelo fue rechazado por el test de Chi-cuadrado, debido a la asimetría positiva de los datos —la mayoría de los estados presenta emisiones bajas de NOx y unos pocos estados concentran emisiones muy altas. Por eso se reemplazó por un modelo Exponencial, adecuado para variables continuas positivas con una cola larga a la derecha como ésta.

8 Parámetros

# Calculo del parametro Exponencial (metodo de momentos)
media_edo <- mean(EMIS_NOX_edo)
lambda_1  <- 1 / media_edo

media_edo

## [1] 21944217

lambda_1

## [1] 0.00000004557009

9 Sobreposición de la Realidad con el Modelo

hist(EMIS_NOX_edo,
     freq = FALSE,
     ylim = c(0, max(hist(EMIS_NOX_edo, plot = FALSE)$density,
                      dexp(seq(min(EMIS_NOX_edo), max(EMIS_NOX_edo), length.out = 200),
                           rate = lambda_1)) * 1.05),
     main = "Gráfica N°2: Comparación de la realidad con el modelo Exponencial\nde EMIS_NOX (nivel-estado)",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de NOx (toneladas)",
     col = "lightgray",
     border = "black")

x <- seq(min(EMIS_NOX_edo), max(EMIS_NOX_edo), length.out = 200)
curve(dexp(x, rate = lambda_1),
      col = "blue",
      lwd = 2,
      add = TRUE)

10 Test de Bondad

k_edo <- ceiling(1 + 3.322 * log10(n_edo))
k_edo

## [1] 7

breaks_edo <- seq(min(EMIS_NOX_edo), max(EMIS_NOX_edo) + 1, length.out = k_edo + 1)

cortes_edo <- cut(EMIS_NOX_edo, breaks = breaks_edo, include.lowest = TRUE, right = FALSE)

Fo_1 <- as.vector(table(cortes_edo))
Fo_1

## [1] 17 13  9  6  1  0  1

Fe_1 <- n_edo * diff(pexp(breaks_edo, rate = lambda_1))
Fe_1

## [1] 20.6941945 11.5823944  6.4825843  3.6282566  2.0307096  1.1365738  0.6361323

# Expresar Fo y Fe en porcentaje para comparar
Fo_pct <- (Fo_1/n_edo)*100
Fe_pct <- (Fe_1/n_edo)*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 Exponencial de EMIS_NOX",
     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] 95.29

# Chi-cuadrado
grados_libertad_1 <- length(Fo_1) - 1 - 1   # 1 parametro estimado (rate)
grados_libertad_1

## [1] 5

nivel_significancia <- 0.95

x2_1 <- sum((Fo_1 - Fe_1)^2 / Fe_1)
x2_1

## [1] 5.228802

umbral_aceptacion_1 <- qchisq(nivel_significancia, grados_libertad_1)
umbral_aceptacion_1

## [1] 11.0705

# Tabla resumen
Variable <- c("Emisiones de NOx (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 Exponencial)")) %>%
  tab_source_note(source_note = md("Autor: Luis Cruz"))

Tabla N°2
Resumen del test de bondad de ajuste (modelo Exponencial)
Variable	Test Pearson (%)	Chi Cuadrado	Umbral de aceptación
Emisiones de NOx (ton)	95.29	5.23	11.07
Autor: Luis Cruz

Con Chi² = 5.23 por debajo del umbral de aceptación 11.07 (95% de confianza, gl = 5), el modelo Exponencial se acepta.

11 Cálculo de Probabilidades

¿Cuál es la probabilidad de que un estado productor tenga emisiones de NOx entre el primer y el tercer intervalo de la distribución observada?

q1 <- round(quantile(EMIS_NOX_edo, 0.25), 0)
q3 <- round(quantile(EMIS_NOX_edo, 0.75), 0)

Probabilidad_1 <- (pexp(q3, rate = lambda_1) - pexp(q1, rate = lambda_1)) * 100
Probabilidad_1

##      75% 
## 52.46669

x <- seq(min(EMIS_NOX_edo), max(EMIS_NOX_edo), length.out = 300)
plot(x,
     dexp(x, rate = lambda_1),
     col = "skyblue3",
     lwd = 1,
     type = "l",
     main = "Gráfica N°4: Cálculo de probabilidades",
     ylab = "Densidad de probabilidad",
     xlab = "Emisiones de NOx (toneladas)")

x_section <- seq(q1, q3, length.out = 200)
y_section <- dexp(x_section, rate = lambda_1)

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 Exponencial", "Área de Probabilidad"),
       col = c("skyblue3", "red"),
       lwd = 2,
       cex = 0.7)

texto_prob <- paste0("Probabilidad = ", round(Probabilidad_1, 2), " %")
text(x = median(EMIS_NOX_edo), y = max(dexp(x, rate = lambda_1))*0.9,
     labels = texto_prob,
     col = "black",
     cex = 0.8,
     font = 2)

# De los estados productores con dato de NOx, cuantos tendrian emisiones en ese rango
cantidad_1 <- (pexp(q3, rate = lambda_1) - pexp(q1, rate = lambda_1)) * n_edo
cantidad_1

##      75% 
## 24.65934

12 Intervalo de Confianza

x_media <- mean(EMIS_NOX_edo)
x_media

## [1] 21944217

sigma <- sd(EMIS_NOX_edo)
sigma

## [1] 18995469

n_ic <- length(EMIS_NOX_edo)
n_ic

## [1] 47

e <- (sigma/sqrt(n_ic)) * qt(0.975, n_ic - 1)
e

## [1] 5577279

li <- x_media - e
li

## [1] 16366939

ls <- x_media + e
ls

## [1] 27521496

tabla_media <- data.frame(Intervalo = sprintf("P [%s< μ <%s] = 95%%", format(round(li,0), big.mark=","), format(round(ls,0), big.mark=",")))

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 [16,366,939< μ <27,521,496] = 95%
Autor: Luis Cruz

13 Conclusiones

El comportamiento de EMIS_NOX se explica con un modelo Exponencial de parámetro rate = 4.56e-08. Podemos afirmar con un 95% de confianza que la media aritmética real de EMIS_NOX se encuentra entre 16,366,939 y 27,521,496 toneladas, y una desviación estándar de 18,995,469 toneladas.