0.- Carga de librerías

library(gt)
## Warning: package 'gt' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

1.- Cargar datos

datos <- read.csv("C:/Users/JOSELYN/Documents/UNIVERSIDAD 3-4 SEMESTRE/ESTADISTICA/MAPA QGIS/Datos Cambiados.csv", 
                  header = TRUE, dec = ".", sep = ",")

2.- Variable aleatoria

# Extracción de la variable Dióxido de azufre (SO2)
# Se eliminan los guiones "-" que representan valores inexistentes
dioxidoazufre <- datos$SO2[datos$SO2 != "-"]
# Conversión a numérico y limpieza de posibles NA generados
dioxidoazufre <- as.numeric(dioxidoazufre)
dioxidoazufre <- dioxidoazufre[!is.na(dioxidoazufre)]
# Justificación porque la var<iable es continua: porque puede tomar cualquier valor dentro de los números reales positivos, incluido el cero, sin presentar saltos en su dominio o la variable cambia de manera gradual y no pasa bruscamente de un valor a otro.

3.- Extraer TDF y GDF

# Forzamos exactamente 10 intervalos sin decimales (de 0 a 200 en saltos de 20)
cortes_exactos <- seq(0, 200, by = 20)

Histograma_dioxidoazufre <- hist(dioxidoazufre, breaks = cortes_exactos, plot = FALSE) 
breaks <- Histograma_dioxidoazufre$breaks
Li <- breaks[1:(length(breaks)-1)]
Ls <- breaks[2:length(breaks)]
ni <- Histograma_dioxidoazufre$counts
N <- length(dioxidoazufre)

# Cálculo de frecuencias y probabilidad
hi <- (ni / N) * 100
p_s <- ni / N  

TDF_dioxidoazufresimplificado <- data.frame(
  # Al ser cortes enteros, ya no necesitamos la función round() aquí
  Intervalo = paste0("[", Li, " - ", Ls, ")"),
  ni = ni,
  hi = round(hi, 2),
  p_s = round(p_s, 2) # Reducimos la probabilidad a solo 2 decimales
)

colnames(TDF_dioxidoazufresimplificado) <- c(
  "Intervalo",
  "ni",
  "hi(%)",
  "p(s)"
)

totales <- data.frame(
  Intervalo = "Totales",
  ni = sum(ni),           # suma total de ni
  hi = sum(hi),           # debe sumar 100
  p_s = sum(round(p_s, 2)) # Ajustado para que cuadre visualmente con la columna
)

colnames(totales) <- c(
  "Intervalo",
  "ni",
  "hi(%)",
  "p(s)"
)

# Agregar al final de la tabla
TDF_dioxidoazufresimplificado <- rbind(TDF_dioxidoazufresimplificado, totales)

TDF_dioxidoazufresimplificado %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 1*"),
    subtitle = md("**Distribución de frecuencia simplificada de concentración de dióxido de azufre, estudio calidad del aire en India entre 2015-2020**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Ariel Chiluisa\nFuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india")
  ) %>%
  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",
    
    # Líneas verticales para los márgenes exteriores
    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),
    
    # Líneas verticales para separar las columnas internamente
    column_labels.vlines.color = "black",
    column_labels.vlines.style = "solid",
    column_labels.vlines.width = px(1),
    table_body.vlines.color = "black",
    table_body.vlines.style = "solid",
    table_body.vlines.width = px(1)
  )
Tabla Nro. 1
Distribución de frecuencia simplificada de concentración de dióxido de azufre, estudio calidad del aire en India entre 2015-2020
Intervalo ni hi(%) p(s)
[0 - 20) 21244 82.74 0.83
[20 - 40) 2735 10.65 0.11
[40 - 60) 953 3.71 0.04
[60 - 80) 287 1.12 0.01
[80 - 100) 174 0.68 0.01
[100 - 120) 128 0.50 0.00
[120 - 140) 83 0.32 0.00
[140 - 160) 40 0.16 0.00
[160 - 180) 28 0.11 0.00
[180 - 200) 5 0.02 0.00
Totales 25677 100.00 1.00
Autor: Ariel Chiluisa Fuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india 
# Histograma porcentual de la variable
par(mar = c(5.1, 4.1, 4.1, 2.1)) # Ajuste de márgenes para que no se corten los textos

# Se crea el lienzo base (ocultando las barras por defecto con col=NA y border=NA)
plot(Histograma_dioxidoazufre,
     freq = FALSE,
     col = NA, 
     border = NA,
     main = "Gráfica N°1: Distribución de la Concentración de Dióxido de Azufre\nen estudio calidad del aire en India, 2015-2020",
     xlab = "Dióxido de azufre (µg/m³)",
     ylab = "Porcentaje",
     xaxt = "n",
     yaxt = "n",
     cex.main = 1,
     ylim = c(0, max(hi) * 1.1))

# Cuadrícula para facilitar la lectura
abline(v = breaks, col = "gray70", lty = 2, lwd = 0.8)
abline(h = pretty(c(0, max(hi))), col = "gray70", lty = 2, lwd = 0.8)

# Barras porcentuales superpuestas (usando hi)
rect(
  breaks[-length(breaks)],
  0,
  breaks[-1],
  hi,
  col = "orange", # Color adaptado a tu variable SO2
  border = "black"
)

# Eje X con los límites de los intervalos enteros
axis(1,
     at = breaks,
     labels = breaks, # Se eliminó el round() para mostrar los enteros limpios
     las = 1,
     cex.axis = 0.9)

# Eje Y en porcentaje
axis(2,
     at = pretty(c(0, max(hi))),
     labels = pretty(c(0, max(hi))),
     las = 1)

# Borde exterior de la gráfica
box()

TDF Y GDF AJUSTADO

# TDF
# Filtramos los datos para eliminar los valores atípicos (mayores a 140)
# Creamos una nueva variable para no dañar tus datos originales
dioxidoazufre_ajustado <- dioxidoazufre[dioxidoazufre <= 140]

# Definimos los nuevos cortes exactos de 0 a 140 en saltos de 20
cortes_ajustados <- seq(0, 140, by = 20)

# Creamos el histograma interno con los datos filtrados
Histograma_ajustado <- hist(dioxidoazufre_ajustado, breaks = cortes_ajustados, plot = FALSE) 

breaks_ajust <- Histograma_ajustado$breaks
Li_ajust <- breaks_ajust[1:(length(breaks_ajust)-1)]
Ls_ajust <- breaks_ajust[2:length(breaks_ajust)]
ni_ajust <- Histograma_ajustado$counts
N_ajustado <- length(dioxidoazufre_ajustado)

# Cálculo de frecuencias y probabilidad con la nueva cantidad de datos
hi_ajust <- (ni_ajust / N_ajustado) * 100
p_s_ajust <- ni_ajust / N_ajustado  

TDF_ajustada <- data.frame(
  Intervalo = paste0("[", Li_ajust, " - ", Ls_ajust, ")"),
  ni = ni_ajust,
  hi = round(hi_ajust, 2),
  p_s = round(p_s_ajust, 2)
)

colnames(TDF_ajustada) <- c(
  "Intervalo",
  "ni",
  "hi(%)",
  "p(s)"
)

totales_ajust <- data.frame(
  Intervalo = "Totales",
  ni = sum(ni_ajust),
  hi = sum(hi_ajust),
  p_s = sum(round(p_s_ajust, 2))
)

colnames(totales_ajust) <- c(
  "Intervalo",
  "ni",
  "hi(%)",
  "p(s)"
)

# Agregar la fila de totales al final
TDF_ajustada <- rbind(TDF_ajustada, totales_ajust)

# Imprimir la tabla con el mismo formato profesional
TDF_ajustada %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 2*"),
    subtitle = md("**Distribución de frecuencia especifica de concentración de dióxido de azufre**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Ariel Chiluisa\nFuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india")
  ) %>%
  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),
    
    column_labels.vlines.color = "black",
    column_labels.vlines.style = "solid",
    column_labels.vlines.width = px(1),
    table_body.vlines.color = "black",
    table_body.vlines.style = "solid",
    table_body.vlines.width = px(1)
  )
Tabla Nro. 2
Distribución de frecuencia especifica de concentración de dióxido de azufre
Intervalo ni hi(%) p(s)
[0 - 20) 21244 82.97 0.83
[20 - 40) 2735 10.68 0.11
[40 - 60) 953 3.72 0.04
[60 - 80) 287 1.12 0.01
[80 - 100) 174 0.68 0.01
[100 - 120) 128 0.50 0.00
[120 - 140) 83 0.32 0.00
Totales 25604 100.00 1.00
Autor: Ariel Chiluisa Fuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india 
# GDF Especifica

par(mar = c(5.1, 4.1, 4.1, 2.1)) # Ajuste de márgenes para que no se corten los textos

# Se crea el lienzo base ocultando las barras por defecto con col=NA y border=NA
plot(Histograma_ajustado,
     freq = FALSE,
     col = NA, 
     border = NA,
     main = "Gráfica N°2: Distribución porcentual específica de la Concentración\nde Dióxido de Azufre, India (2015-2020)",
     xlab = "Dióxido de azufre (µg/m³)",
     ylab = "Porcentaje",
     xaxt = "n",
     yaxt = "n",
     cex.main = 1,
     ylim = c(0, max(hi_ajust) * 1.1))

# Cuadrícula para facilitar la lectura, ajustada a los nuevos límites
abline(v = cortes_ajustados, col = "gray70", lty = 2, lwd = 0.8)
abline(h = pretty(c(0, max(hi_ajust))), col = "gray70", lty = 2, lwd = 0.8)

# Barras porcentuales superpuestas (usando la nueva hi_ajust)
rect(
  cortes_ajustados[-length(cortes_ajustados)],
  0,
  cortes_ajustados[-1],
  hi_ajust,
  col = "orange", # Mantenemos el color adaptado a la variable SO2
  border = "black"
)

# Eje X con los límites de los intervalos exactos (0 a 140)
axis(1,
     at = cortes_ajustados,
     labels = cortes_ajustados, 
     las = 1,
     cex.axis = 0.9)

# Eje Y en porcentaje, autoajustado a las nuevas proporciones
axis(2,
     at = pretty(c(0, max(hi_ajust))),
     labels = pretty(c(0, max(hi_ajust))),
     las = 1)

# Borde exterior de la gráfica
box()

4.- Conjetura

#Se conjetura que la concentración de dióxido de azufre (SO2) sigue un modelo  de probabilidad exponencial, ya que la variable en su histograma muestra que la mayor parte de los datos está concentrada en el extremo izquierdo (bajos valores) y la frecuencia disminuye rápidamente conforme aumentan los valores.

5.- Gráfica de sobreposición

# Parámetros exponenciales
media <- mean(dioxidoazufre_ajustado)
lambda <- 1 / media

# Histograma porcentual (lienzo vacío)
hist(dioxidoazufre_ajustado,
     breaks = cortes_ajustados,
     freq = FALSE,
     col = NA,
     border = NA,
     main = "Gráfica Nro 3\nComparación de la Realidad y el Modelo Exponencial\nde la Concentración de Dióxido de Azufre en India 2015–2020",
     xlab = "Dióxido de azufre (µg/m³)",
     ylab = "Porcentaje",
     xaxt = "n",
     yaxt = "n",
     ylim = c(0, max(hi_ajust) * 1.25))

# Cuadrícula
abline(v = cortes_ajustados, col = "gray80", lty = 2)
abline(h = pretty(c(0, max(hi_ajust))), col = "gray80", lty = 2)

# Barras 
rect(
  cortes_ajustados[-length(cortes_ajustados)],
  0,
  cortes_ajustados[-1],
  hi_ajust,
  col = "orange", # Manteniendo el color para SO2
  border = "black"
)

# Curva exponencial convertida a porcentaje
# x va desde 0 hasta el máximo de tus cortes (140)
x <- seq(0, 140, length = 1000)

# Multiplicamos por 20 (amplitud del intervalo) y por 100 (porcentaje)
y <- dexp(x, rate = lambda) * 20 * 100

lines(x, y,
      col = "darkred", # Color oscuro para que contraste bien con la barra
      lwd = 3)

# Ejes
axis(1,
     at = cortes_ajustados,
     labels = cortes_ajustados)

axis(2,
     at = pretty(c(0, max(hi_ajust))),
     las = 1)

# Añadir leyenda al gráfico
legend("topright", 
       legend = c("Datos reales (FO)", "Modelo Exponencial (FE)"), 
       col = c("orange", "darkred"), 
       lty = c(NA, 1),           # NA para que no dibuje línea en el cuadro del histograma
       pch = c(22, NA),          # 22 es un cuadrado (representa la barra)
       pt.bg = c("orange", NA), 
       pt.cex = 2,               # Tamaño del símbolo
       lwd = c(1, 3),            # Grosor de la línea
       bty = "n")                # "n" para quitar el borde del recuadro de la leyenda
box()

6.- Test de Bondad

## TEST DE PEARSON

fo_pearson <- ni_ajust
# N_ajustado ya está calculado en el paso de la TDF ajustada

# Frecuencias esperadas usando 'lambda' (el nombre que creaste arriba)
fe_pearson <- N_ajustado * (pexp(Ls_ajust, rate = lambda) - pexp(Li_ajust, rate = lambda))

Coef_Pearson <- cor(fo_pearson, fe_pearson) * 100

# Resultado de Pearson
cat("Coeficiente de Pearson (%):", Coef_Pearson, "\n")
## Coeficiente de Pearson (%): 99.30358
##TEST CHI-CUADRADO

# Total de datos ajustados
N_chi <- sum(ni_ajust)

# Frecuencia relativa observada
fo <- ni_ajust / N_chi

# Número de intervalos
k <- length(fo)

# Probabilidades teóricas por intervalo usando 'lambda'
P <- c(0)

for (i in 1:k) {
  P[i] <- pexp(Ls_ajust[i], rate = lambda) - pexp(Li_ajust[i], rate = lambda)
}

fe <- P  # Probabilidades esperadas

# Chi-cuadrado
Chi_Calculado <- sum((fo - fe)^2 / fe)

# Grados de libertad
gl <- k - 1

# Valor crítico (α = 0.05)
Chi_Critico <- qchisq(0.95, df = gl)

# Resultados
cat("Chi Calculado:", Chi_Calculado, "\n")
## Chi Calculado: 0.1387734
cat("Chi Crítico:", Chi_Critico, "\n")
## Chi Crítico: 12.59159
# Decisión
if (Chi_Calculado < Chi_Critico) {
  print("Evalua H0: El modelo exponencial es adecuado.")
} else {
  print("Se rechaza H0: El modelo exponencial no es adecuado.")
}
## [1] "Evalua H0: El modelo exponencial es adecuado."

7.- Calculo de probabilidad

#--------------------------------------
# PREGUNTA DE PORCENTAJE
# ¿Cuál es la probabilidad de que la concentración de SO2 no supere los 40 µg/m³?

# PROBABILIDAD: P(SO2 ≤ 40)

prob_40 <- pexp(40, rate = lambda)
prob_40_porcentaje <- prob_40 * 100

cat("Probabilidad de que el SO2 no supere los 40 µg/m³:",
    round(prob_40_porcentaje, 2), "%\n")
## Probabilidad de que el SO2 no supere los 40 µg/m³: 94.12 %
# PREGUNTA DE CANTIDAD
# ¿Cuántos días en 10 años se espera que el SO2 no supere los 40 µg/m³?

dias_10_anios <- prob_40 * (10 * 365)

cat("Cantidad esperada de días:", round(dias_10_anios), "días\n")
## Cantidad esperada de días: 3435 días
# PREGUNTA DE INTERVALO
# ¿Cuál es la probabilidad de que la concentración de SO2 se encuentre entre 20 y 40 µg/m³?

prob_20_40 <- (pexp(40, rate = lambda) -
                 pexp(20, rate = lambda)) * 100

cat("Probabilidad entre 20 y 40 µg/m³:",
    round(prob_20_40,2), "%\n")
## Probabilidad entre 20 y 40 µg/m³: 18.37 %
# DEMOSTRACIÓN GRÁFICA

x <- seq(0, 140, by = 0.01)

y <- dexp(x, rate = lambda)
y <- y * 100 * 20     # 20 = amplitud del intervalo

plot(x, y,
     col = "orange",
     lwd = 2,
     type = "l",
     xlim = c(0,140),
     ylim = c(0,max(y)),
     main = "Gráfica N°4: Cálculo de Probabilidad",
     ylab = "Densidad de probabilidad",
     xlab = "Dióxido de azufre (µg/m³)",
     xaxt = "n")

axis(1,
     at = seq(0,140,20),
     las = 1)
# Área de probabilidad P(SO2 ≤ 40)

x_section <- seq(0,40,0.01)

y_section <- dexp(x_section, rate=lambda)
y_section <- y_section*100*20

lines(x_section,
      y_section,
      col="darkgreen",
      lwd=2)

polygon(c(x_section,rev(x_section)),
        c(y_section,rep(0,length(y_section))),
        col=rgb(0,0.6,0,0.35),
        border=NA)
# Área de probabilidad P(20 ≤ SO2 ≤ 40)

x_section2 <- seq(20,40,0.01)

y_section2 <- dexp(x_section2, rate=lambda)
y_section2 <- y_section2*100*20

lines(x_section2,
      y_section2,
      col="blue",
      lwd=2)

polygon(c(x_section2,rev(x_section2)),
        c(y_section2,rep(0,length(y_section2))),
        col=rgb(0,0,1,0.35),
        border=NA)
legend("topright",
       legend=c("Modelo exponencial",
                "Área P(SO2 ≤ 40)",
                "Área P(20 ≤ SO2 ≤ 40)"),
       col=c("orange",
             "darkgreen",
             "blue"),
       lwd=2,
       bty="o",
       bg="white")

grid()

# Intervalo de Confianza (Tabla)
# Calculamos con los datos ajustados para mantener la precisión del modelo
media_ic <- mean(dioxidoazufre_ajustado)
desviacion_ic <- sd(dioxidoazufre_ajustado)
n_ic <- length(dioxidoazufre_ajustado)

error <- 1.96 * (desviacion_ic / sqrt(n_ic)) 
limite_inferior <- round(media_ic - error, 2)
limite_superior <- round(media_ic + error, 2)

tabla_intervalo <- data.frame(Intervalo = paste0("P [", limite_inferior, " < µ < ", limite_superior, "] = 95%"))

# Tabla con el mismo formato estético que las anteriores
tabla_intervalo %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 3*"),
    subtitle = md("**Intervalo de confianza de la concentración de dióxido de azufre**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Ariel Chiluisa\nFuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india")
  ) %>%
  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",
    
    # Bordes laterales para igualar tus Tablas 1 y 2
    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. 3
Intervalo de confianza de la concentración de dióxido de azufre
Intervalo
P [13.92 < µ < 14.32] = 95%
Autor: Ariel Chiluisa Fuente: https://www.kaggle.com/datasets/rohanrao/air-quality-data-in-india

8.- Conclusión

En conclusión:

El análisis de la concentración de SO2 (µg/m³) revela un modelo exponencial (λ = 0.069), validado mediante las pruebas de Pearson y Chi-cuadrado. Con un nivel de confianza del 95%, se estima que el promedio real de concentración se sitúa en el rango [13.92, 14.32] µg/m³. Si bien existe una probabilidad del 6.26% de superar el límite legal de 40 µg/m³, la distribución demuestra que las concentraciones predominantes son bajas, confirmando una calidad del aire generalmente estable en la región durante el periodo 2015-2020.