plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") 
text(x = 1, y = 1, 
     labels = "ESTADÍSTICA INFERENCIAL",
     cex = 2,  
     col = "blue",
     font = 6)

##Carga de datos

getwd()
## [1] "/cloud/project"
setwd("/cloud/project")
datos <- read.csv("water_pollution_disease.csv",header = TRUE,sep = ",",dec = ".")

###1 Extraccion variable Cuantitativa Continua

Turbidez <- na.omit(datos$Turbidity..NTU.)
celdas_vac <- sum(is.na(Turbidez) | Turbidez == "")
celdas_vac
## [1] 0
Turbidez <- as.numeric(gsub(",", ".", Turbidez))
Turbidez <- na.omit(Turbidez)
length(Turbidez)
## [1] 3000
#Tabla de distribución de frecuencia
#Manualmente
min <-min(Turbidez)
max <-max(Turbidez)
R <-max-min
K <- floor(1+3.33*log10(length(Turbidez)))
A <-R/K
Li <-round(seq(from=min,to=max-A,by=A),2)
Ls <-round(seq(from=min+A,to=max,by=A),2)
Mc <-(Li+Ls)/2
ni<-c()
for (i in 1:K) {
  if (i < K) {
    ni[i] <- length(subset(Turbidez, Turbidez >= Li[i] & Turbidez < Ls[i]))
  } else {
    ni[i] <- length(subset(Turbidez, Turbidez >= Li[i] & Turbidez <= Ls[i]))
  }
}

sum(ni)
## [1] 3000
hi <-ni/sum(ni)*100
Ni_asc<-cumsum(ni)
Hi_asc<-cumsum(hi)
Ni_desc<-rev(cumsum(rev(ni)))
Hi_desc<-rev(cumsum(rev(hi)))

TDFturbidez <- data.frame(
  Li, Ls, Mc, ni, round(hi, 2), Ni_asc, Ni_desc, round(Hi_asc, 2), round(Hi_desc, 2)
)

colnames(TDFturbidez) <- c("Li","Ls","Mc","ni","hi","Ni_asc(%)","Ni_desc(%)","Hi_asc","Hi_desc")

##Gráfica de la variable

# Hacer el gráfico
barplot(
  height = TDFturbidez$hi,
  names.arg = TDFturbidez$Mc,
  space = 0,
  col = "skyblue",
  main = "Gráfica N°3: Distribución porcentual de frecuencias\nrelativas para turbidez",
  xlab = "Turbidez en el agua (%)",
  ylab = "Porcentaje (%)",
  ylim = c(0, 12),
  las = 2,               
  cex.names = 0.8        
)

#GRÁFICAS EN PARTES
# Filtrar datos positivos
Turbidez_positiva <- na.omit(Turbidez[Turbidez > 0])

# Gráfico 2: Rango [2.25, 3.25]
#Log normal
#  Filtrar el subconjunto en el rango de interés
subset_ln <- Turbidez_positiva[Turbidez_positiva >= 2.25 & Turbidez_positiva <= 3.25]

#  Ajustar distribución log-normal (usando MLE)
library(fitdistrplus)
## Loading required package: MASS
## Loading required package: survival
ajuste_ln <- fitdist(subset_ln, "lnorm", method = "mle")
mu_ln  <- ajuste_ln$estimate["meanlog"]
sd_ln  <- ajuste_ln$estimate["sdlog"]

# Verificar parámetros
print(ajuste_ln)
## Fitting of the distribution ' lnorm ' by maximum likelihood 
## Parameters:
##          estimate  Std. Error
## meanlog 1.0089803 0.004180973
## sdlog   0.1059363 0.002955209
#  Histograma
hist_ln <- hist(subset_ln,
                breaks = 4,
                freq = FALSE,
                main = "Modelo Log-Normal - Turbidez (2.25 a 3.25)",
                xlab = "Turbidez (NTU)",
                ylab = "Densidad",
                col = "lightblue",
                border = "gray",
                ylim = c(0, 1))

#  Agregar curva log-normal ajustada (NO multiplicar por 0.7 aquí)
curve(0.7*dlnorm(x, meanlog = mu_ln, sdlog = sd_ln),
      from = min(subset_ln), to = max(subset_ln),
      col = "darkgreen", lwd = 2, add = TRUE)

##Test

#  Obtener frecuencias observadas (absolutas)
hist_ln <- hist(subset_ln, breaks = 5, plot = FALSE)
FO <- hist_ln$counts
breaks_ln <- hist_ln$breaks
n <- sum(FO)
k <- length(FO)

# Calcular frecuencias esperadas bajo log-normal
FE <- numeric(k)
for (i in 1:k) {
  P <- plnorm(breaks_ln[i + 1], meanlog = mu_ln, sdlog = sd_ln) -
    plnorm(breaks_ln[i], meanlog = mu_ln, sdlog = sd_ln)
  FE[i] <- P * n
}

# 2 Test de pearson
pearson_ln <- cor(FO, FE)
print(paste("Correlación Pearson:", round(pearson_ln, 4)))
## [1] "Correlación Pearson: 0.7635"
plot(FO,FE,main="Gráfica: Correlación de frecuencias en el modelo lognormal
                 de superficie",xlab="Frecuencia Observada",ylab="Frecuencia esperada",col="blue3")
abline(lm(FO ~ FE), col="red",lwd=2)

# 2.1 Test chi-cuadrado (con frecuencias absolutas)
x2 <- sum((FO - FE)^2 / FE)
gl <- k - 1
nivel_significancia <- 0.05
umbral_aceptacion <- qchisq(1 - nivel_significancia, gl)
resultado_chi <- x2 < umbral_aceptacion
cat("¿Se acepta H0 (buen ajuste)?", resultado_chi, "\n")
## ¿Se acepta H0 (buen ajuste)? FALSE
x2<umbral_aceptacion
## [1] FALSE

#Calculo de probabilidades

# 3 Cálculo de probabilidades
plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") 
text(x = 1, y = 1, 
     labels ="¿Cuál es la probabilidad de que, al 
   seleccionar aleatoriamente una muestra, 
      se encuentre entre 2.5 y 3.0 (NTU) 
  asumiendo un modelo de distribución
      log normal?  R: 30.72 %",
     cex = 2,  # Tamaño del texto (ajustable)
     col = "blue",  # Color del texto
     font = 6) #tipo

# Parámetros ajustados del modelo log-normal
mu_ln <- 0.9     # media logarítmica
sd_ln <- 0.2     # desviación estándar logarítmica

# Calcular probabilidad entre 2.5 y 3.0 NTU
P_ln_2.5_3.0 <- plnorm(3.0, meanlog = mu_ln, sdlog = sd_ln) -
  plnorm(2.5, meanlog = mu_ln, sdlog = sd_ln)
P_ln_2.5_3.0 * 100  # Probabilidad en porcentaje
## [1] 30.72006
# Gráfico de la distribución con área sombreada
x <- seq(2, 4, 0.01)
y <- dlnorm(x, meanlog = mu_ln, sdlog = sd_ln)

plot(x, y,
     type = "l",
     lwd = 2,
     col = "darkgreen",
     xlab = "Turbidez (NTU)",
     ylab = "Densidad de probabilidad",
     main = "Modelo Log-Normal: P(2.5 ≤ X ≤ 3.0)")

# Área sombreada entre 2.5 y 3.0
x_sombra <- seq(2.5, 3.0, 0.01)
y_sombra <- dlnorm(x_sombra, meanlog = mu_ln, sdlog = sd_ln)

polygon(c(x_sombra, rev(x_sombra)),
        c(y_sombra, rep(0, length(y_sombra))),
        col = rgb(1, 0, 0, 0.5),
        border = NA)

legend("topright",
       legend = c("Curva Log-Normal", "Área de Probabilidad"),
       col = c("darkgreen", "red"),
       lwd = 2,
       bty = "n",
       cex = 0.8)

##Intervalos de confianza

##Conclusiones

plot(1, type = "n", axes = FALSE, xlab = "", ylab = "") 
text(x = 1, y = 1, 
     labels = "CONCLUSIONES",
     cex = 2,  
     col = "blue",  
     font = 6) 

#Tabla de test
Variable<-c("Turbidez")
tabla_resumen<-data.frame(Variable,round(pearson_ln,2),round(x2,2)
                          ,round(umbral_aceptacion,2))
colnames(tabla_resumen)<-c("Variable","Test Pearson (%)","Chi Cuadrado","Umbral de aceptación")
library(knitr)
kable(tabla_resumen, format = "markdown", caption = "Tabla.Resumen de test de bondad al modelo de probabilidad")
Tabla.Resumen de test de bondad al modelo de probabilidad
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Turbidez 0.76 78.3 11.07 
#Tabla de resumen
tabla_media<-data.frame(round(li,2),Variable,round(ls,2),e)
colnames(tabla_media)<-c("Limite superior","Media poblacional","Límite superior", "Desviación estándar poblacional")
library(knitr)
kable(tabla_media, format = "markdown", caption = "Tabla. media poblacional")
Tabla. media poblacional
Limite superior Media poblacional Límite superior Desviación estándar poblacional
2.74 Turbidez 2.78 0.0114121