##Modelo Log-normal COD

# ============================================
# 3.X MODELO LOG-NORMAL (BASE CHINA - COD)
# Variable: COD_mg_L
# ============================================
library(gt)
library(dplyr)
## 
## 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. Carga de datos (CSV)
datos <- read.csv("china_water_pollution_data.csv", header = TRUE, sep = ",")

# 1. Extracción de la variable y justificación
# Justificación: La DQO (COD, mg/L) es una variable cuantitativa continua porque representa una concentración que puede tomar valores reales positivos dentro de un rango y no está limitada a valores enteros específicos.Es decir r, Su dominio se asocia al conjunto de los números reales.

COD <- as.numeric(datos$COD_mg_L)

# Limpieza básica (lognormal requiere valores positivos y finitos)
COD <- COD[is.finite(COD)]
COD <- COD[COD > 0]

# 2. TDF simplificada (como tu ejemplo)
Histograma_COD <- hist(COD, plot = FALSE)

breaks <- Histograma_COD$breaks
Li <- breaks[1:(length(breaks)-1)]
Ls <- breaks[2:length(breaks)]
ni <- Histograma_COD$counts

n <- length(COD)
hi <- (ni / n) * 100

TDF_COD <- data.frame(
  Intervalo = paste0("[", round(Li, 2), " - ", round(Ls, 2), ")"),
  ni = ni,
  hi = round(hi, 2)
)

colnames(TDF_COD) <- c("Intervalo", "ni", "hi(%)")

totaless <- data.frame(
  Intervalo = "Totales",
  ni = sum(ni),
  hi = round(sum(hi), 2)
)

colnames(totaless) <- c("Intervalo", "ni", "hi(%)")

# Agregar al final de la tabla
TDF_COD <- rbind(TDF_COD, totaless)

# TABLA CREACIÓN (igual estilo)
TDF_COD %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 1*"),
    subtitle = md("**Distribución de frecuencia simplificada de la DQO (COD) en la base de contaminación de agua**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 1")
  ) %>%
  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 Nro. 1
Distribución de frecuencia simplificada de la DQO (COD) en la base de contaminación de agua
Intervalo ni hi(%)
[2 - 4) 1 0.03
[4 - 6) 8 0.27
[6 - 8) 17 0.57
[8 - 10) 41 1.37
[10 - 12) 106 3.53
[12 - 14) 169 5.63
[14 - 16) 296 9.87
[16 - 18) 407 13.57
[18 - 20) 474 15.80
[20 - 22) 462 15.40
[22 - 24) 362 12.07
[24 - 26) 310 10.33
[26 - 28) 183 6.10
[28 - 30) 103 3.43
[30 - 32) 35 1.17
[32 - 34) 17 0.57
[34 - 36) 7 0.23
[36 - 38) 2 0.07
Totales 3000 100.00
Autor: Grupo 1 
# 3. GDF - Histograma porcentual (barplot hi%)
TDF_COD_graf <- TDF_COD[TDF_COD$Intervalo != "Totales", ]

par(mar = c(9, 5, 4, 2))

post <- barplot(
  TDF_COD_graf$`hi(%)`,
  space = 0,
  col = "salmon",
  main = "Gráfica N°1: Distribución porcentual de DQO (COD),\nEstudio de contaminación de agua en China en 2023",
  xlab = "COD (mg/L)",
  ylab = "Porcentaje (%)",
  ylim = c(0, 100),
  xaxt = "n"
)

# Eje x con cortes (breaks) 
limites <- c(post[1] - diff(post)[1]/2, post + diff(post)[1]/2)

axis(
  side = 1,
  at = limites,
  labels = round(breaks, 2),
  las = 2,
  cex.axis = 0.8,
  tck = -0.02
)

# 4. CONJETURA
# Al observar el histograma porcentual, se aprecia que los valores de COD se concentran
# en rangos bajos y se extienden hacia valores altos con una cola a la derecha.
# Esta asimetría positiva y el hecho de que COD sea estrictamente positivo sugieren
# que la variable puede aproximarse a un modelo Log-normal.

# 5. CÁLCULO DE PARÁMETROS DISTRIBUCIÓN LOG-NORMAL
min(COD)
## [1] 2.09
log_COD <- log(COD)
mulog <- mean(log_COD)
sigmalog <- sd(log_COD)

mulog
## [1] 2.960155
sigmalog
## [1] 0.2762394
# 6. GDF - HISTOGRAMA DENSIDAD DE PROBABILIDAD + CURVA LOGNORMAL
Histograma_COD2 <- hist(
  COD,
  breaks = breaks,
  col = "salmon",
  freq = FALSE,
  main = "Gráfica N°2: Comparación de la Realidad y el Modelo Log-normal\nde la DQO(COD) en el estudio de la contaminacion del agua en China 2023",
  xlab = "COD (mg/L)",
  ylab = "Densidad de probabilidad",
  cex.main = 0.9,
  xaxt = "n"
)

axis(1, at = breaks)

x <- seq(min(COD), max(COD), by = 0.001)

curve(
  dlnorm(x, meanlog = mulog, sdlog = sigmalog),
  col = "darkblue",
  lwd = 3,
  add = TRUE
)

# 7. TEST DE BONDAD

# 7.1 "Test de Pearson" 
fo <- hist(COD, breaks = breaks, plot = FALSE)$counts

fe <- numeric(length(fo))
for(i in 1:length(fo)){
  fe[i] <- n * (plnorm(breaks[i+1], meanlog = mulog, sdlog = sigmalog) -
                  plnorm(breaks[i], meanlog = mulog, sdlog = sigmalog))
}

Correlación <- cor(fo, fe) * 100
Correlación
## [1] 96.70084
# 7.2 Test Chi-cuadrado
fe_frac <- fe / n
fo_frac <- fo / n

x2 <- sum((fo_frac - fe_frac)^2 / fe_frac)
x2
## [1] 18.96594
k <- length(fo_frac)
gl <- k - 1 - 2  # -2 por estimar mulog y sigmalog
gl
## [1] 15
# Umbral de aceptación 
umbral_aceptacion <- qchisq(0.95, df = gl)
umbral_aceptacion
## [1] 24.99579
x2 < umbral_aceptacion  # TRUE = aprueba
## [1] TRUE
# 8. CÁLCULO DE PROBABILIDAD 
# Ejemplo 1: Probabilidad de que COD esté entre 20 y 30 mg/L
plnorm(30, mulog, sigmalog) - plnorm(20, mulog, sigmalog)
## [1] 0.3935829
# Ejemplo 2: Probabilidad de que COD supere 40 mg/L
1 - plnorm(40, mulog, sigmalog)
## [1] 0.004169585
# Demostración gráfica: cálculo de probabilidad (Lognormal)

x <- seq(min(COD), max(COD), by=0.001)
plot(x, dlnorm(x, mulog, sigmalog), col = "blue", lwd = 2, xaxt = "n",
     main="Gráfica N°3: Cálculo de probabilidad (COD - Lognormal)\nen el estudio de contaminación de agua en China, 2023",
     ylab="Densidad de probabilidad", xlab="COD (mg/L)")

axis(1, at = pretty(range(COD), 10))

# Área roja: COD entre 20 y 30
x_red <- seq(20, 30, by = 0.001)
y_red <- dlnorm(x_red, mulog, sigmalog)
lines(x_red, y_red, col = "red", lwd = 2)
polygon(c(x_red, rev(x_red)), c(y_red, rep(0, length(y_red))),
        col = rgb(1, 0, 0, 0.5), border = NA)

# Área verde: cola alta (automática, percentil 90)
umbral_verde <- as.numeric(quantile(COD, 0.90))

if(max(COD) > umbral_verde){
  x_green <- seq(umbral_verde, max(COD), by = 0.001)
  y_green <- dlnorm(x_green, mulog, sigmalog)
  lines(x_green, y_green, col = "green", lwd = 2)
  polygon(c(x_green, rev(x_green)), c(y_green, rep(0, length(y_green))),
          col = rgb(0, 1, 0, 0.4), border = NA)
}

legend("topright",
       legend = c("Modelo", "P(20 < COD < 30)", paste0("P(COD ≥ P90 = ", round(umbral_verde,2), ")")),
       col = c("blue", "red", "green"),
       lwd = 2,
       pch = c(NA, 15, 15),
       cex = 0.8,
       pt.cex = 1)

# 9. INTERVALO DE CONFIANZA 
media <- mean(COD)
sigma <- sd(COD)
n <- length(COD)

error <- 2 * (sigma / sqrt(n))  # aproximación 

limite_inferior <- round(media - error, 2)
limite_superior <- round(media + error, 2)

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

tabla_intervalo %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nro. 2*"),
    subtitle = md("**Intervalo de confianza de la media de COD (DQO)**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 1")
  ) %>%
  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 Nro. 2
Intervalo de confianza de la media de COD (DQO)
Intervalo
P [19.81 < µ < 20.17] = 95%
Autor: Grupo 1 
# 10. CONCLUSIÓN COD (LOGNORMAL)
# La variable DQO (COD, mg/L) sigue o se explica   con un modelo Log-normal estandar con parámetros µlog = 2.96 y σlog = 0.276 y podemos afirmar con un 95% de confianza que  la media aritmetica de esta variable se encuentra entre 19.81 y 20.17 (mg/L) con una   desviación estandar 5.01 (mg/L).