CARGA DE DATOS Y LIBRERÍAS

#Carga de datos
datos <- read.csv("C:/Users/Grace/OneDrive - Universidad Central del Ecuador/Documentos/dataset_geologico_limpio_80.csv", 
                  header = TRUE, 
                  sep = ",", 
                  dec = ".")

#Librerías
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
library(knitr)
library(gt)
library(moments)

TABLA DE DISTRIBUCIÓN DE CANTIDAD

# LIMPIEZA DE LA VARIABLE ARCILLA

arcilla <- as.numeric(datos$CLAY_PCT)
arcilla <- na.omit(arcilla)
arcilla <- subset(arcilla, arcilla > 0)

# SEPARAR OUTLIERS
caja <- boxplot(arcilla, plot = FALSE)

limite_sup <- caja$stats[5]
limite_inf <- caja$stats[1]

arcilla_outliers <- arcilla[arcilla < limite_inf | arcilla > limite_sup]
arcilla_sin_outliers <- arcilla[arcilla >= limite_inf & arcilla <= limite_sup]

# RESUMEN
cat("Cantidad con outliers:", length(arcilla), "\n")
## Cantidad con outliers: 25966
cat("Cantidad de outliers:", length(arcilla_outliers), "\n")
## Cantidad de outliers: 1084
cat("Cantidad sin outliers:", length(arcilla_sin_outliers), "\n")
## Cantidad sin outliers: 24882
# HISTOGRAMA (BASE DE LA TABLA)

histograma <- hist(arcilla_sin_outliers, breaks = 6, plot = FALSE)

# FRECUENCIA ABSOLUTA
ni <- histograma$counts

# FRECUENCIA RELATIVA
hi <- ni / sum(ni) * 100

# INTERVALOS
intervalos <- paste0(
  "[", round(histograma$breaks[-length(histograma$breaks)], 2),
  ", ",
  round(histograma$breaks[-1], 2),
  ")"
)

# TABLA
tabla_frecuencias <- data.frame(
  Intervalo = intervalos,
  ni = ni,
  hi = round(hi, 2)
)

tabla_frecuencias
##   Intervalo    ni    hi
## 1   [0, 10) 14602 58.68
## 2  [10, 20)  3849 15.47
## 3  [20, 30)  2475  9.95
## 4  [30, 40)  2072  8.33
## 5  [40, 50)  1285  5.16
## 6  [50, 60)   599  2.41
tabla_arcilla_gt <- tabla_frecuencias %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 1**"),
    subtitle = md("**Distribución de frecuencias del contenido de arcilla (%)**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 3")
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black",
    row.striping.include_table_body = TRUE
  )

tabla_arcilla_gt
Tabla N° 1
Distribución de frecuencias del contenido de arcilla (%)
Intervalo ni hi
[0, 10) 14602 58.68
[10, 20) 3849 15.47
[20, 30) 2475 9.95
[30, 40) 2072 8.33
[40, 50) 1285 5.16
[50, 60) 599 2.41
Autor: Grupo 3

GRÁFICA DE DISTRIBUCIÓN DE CANTIDAD

## GRÁFICA DE DISTRIBUCIÓN
histograma <- hist(arcilla_sin_outliers,
                   freq = TRUE,
                   main="Gráfica 1. Distribución de cantidad de arcilla (%) en sedimentos marinos",
                   xlab="Arcilla (%)",
                   ylab="Cantidad",
                   col="blue")

CONJETURA DEL MODELO

## CONJETURA DEL MODELO

# HISTOGRAMA CONTROLADO (CLAVE)

histograma <- hist(arcilla_sin_outliers,
                   breaks = 6,
                   freq = FALSE,
                   main = "Gráfica 2. Comparación con modelo exponencial",
                   xlab = "Arcilla (%)",
                   ylab = "Densidad de probabilidad",
                   col = "lightblue",
                   border = "black")

# PARÁMETROS
media <- mean(arcilla_sin_outliers)
lambda <- 1 / media

# CURVA
x <- seq(min(arcilla_sin_outliers), max(arcilla_sin_outliers), 0.01)
curve(dexp(x, rate = lambda), add = TRUE, col = "black", lwd = 3)

# FRECUENCIAS OBSERVADAS
Fo <- histograma$counts

# FRECUENCIAS ESPERADAS
h <- length(histograma$counts)

P <- c()
for (i in 1:h) {
  P[i] <- pexp(histograma$breaks[i+1], rate = lambda) -
          pexp(histograma$breaks[i], rate = lambda)
}

Fe <- P * length(arcilla_sin_outliers)

TEST DE APROBACIÓN

# =========================
# TEST DE PEARSON
# =========================

n <- length(arcilla_sin_outliers)

Fo <- (Fo/n)*100
Fe <- (Fe/n)*100

plot(Fo, Fe,
     main = "Gráfica 3: Correlación (modelo exponencial - Arcilla)",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     pch = 19,
     col = "blue3")

abline(a = 0, b = 1, col = "red", lwd = 2)

Correlacion <- cor(Fo, Fe) * 100
Correlacion
## [1] 97.213
# =========================
# TEST DE CHI-CUADRADO
# =========================
# GRADOS DE LIBERTAD
gl <- length(histograma$counts)-1

# ESTADÍSTICO
x2 <- sum((Fe - Fo)^2 / Fe)

# UMBRAL
umbral <- qchisq(0.97, gl)

x2
## [1] 11.73207
umbral
## [1] 12.37462
x2 < umbral
## [1] TRUE
# =========================
# TABLA RESUMEN
# =========================
Variable <- c("Arcilla (%)")

tabla_resumen <- data.frame(
  Variable,
  round(Correlacion, 2),
  round(x2, 2),
  round(umbral, 2)
)

colnames(tabla_resumen) <- c(
  "Variable",
  "Test Pearson (%)",
  "Chi Cuadrado",
  "Umbral de aceptación"
)

kable(tabla_resumen)
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Arcilla (%) 97.21 11.73 12.37

CÁLCULO DE PROBABILIDADES

# =========================
# PROBABILIDAD ENTRE 10% y 30%
# =========================

probabilidad_arcilla <- pexp(30, lambda) -
  pexp(10, lambda)

# En porcentaje
probabilidad_arcilla * 100
## [1] 35.88084
# =========================
# GRÁFICA DE PROBABILIDAD
# =========================

# Rango para la curva
x <- seq(min(arcilla_sin_outliers), max(arcilla_sin_outliers), 0.01)

# Curva exponencial
plot(x, dexp(x, lambda),
     col = "skyblue3",
     lwd = 2,
     type = "l",
     main = "Gráfica 4. Cálculo de probabilidades del contenido de arcilla (%)",
     ylab = "Densidad de probabilidad",
     xlab = "Arcilla (%)")

# Área de probabilidad (10–30%)
x_area <- seq(10, 30, 0.01)
y_area <- dexp(x_area, lambda)

# Línea
lines(x_area, y_area, col = "red", lwd = 2)

# Área sombreada
polygon(c(x_area, rev(x_area)),
        c(y_area, rep(0, length(y_area))),
        col = rgb(1, 0, 0, 0.5),
        border = NA)

# Leyenda
legend("topright",
       legend = c("Modelo exponencial", "Área de Probabilidad"),
       col = c("skyblue3", "red"),
       lwd = 2,
       cex = 0.7)

# Texto
texto_prob <- paste0("Probabilidad = ",
                     round(probabilidad_arcilla*100, 2), " %")

text(x = max(arcilla_sin_outliers)*0.6,
     y = max(dexp(x, lambda)) * 0.7,
     labels = texto_prob,
     col = "black",
     cex = 0.9,
     font = 2)

# =========================
# CANTIDAD EN 300 MUESTRAS
# =========================

cantidad_muestras <- probabilidad_arcilla * 300

cantidad_muestras
## [1] 107.6425

INTERVALOS DE CONFIANZA

# MEDIA
x <- mean(arcilla_sin_outliers)
x
## [1] 12.51745
# DESVIACIÓN ESTÁNDAR
sigma <- sd(arcilla_sin_outliers)
sigma
## [1] 14.83708
# TAMAÑO MUESTRAL
n <- length(arcilla_sin_outliers)
n
## [1] 24882
# ERROR ESTÁNDAR
e <- sigma / sqrt(n)
e
## [1] 0.09406016
# INTERVALO 95%
li <- x - 2*e
ls <- x + 2*e

li
## [1] 12.32933
ls
## [1] 12.70557
# Tabla
tabla_media <- data.frame(
  "Límite inferior" = round(li,2),
  "Media poblacional" = round(x,2),
  "Límite superior" = round(ls,2),
  "Error estándar" = round(e,2)
)

tabla_media %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 3**"),
    subtitle = md("**Intervalo de confianza de la media (95%)**")
  )
Tabla N° 3
Intervalo de confianza de la media (95%)
Límite.inferior Media.poblacional Límite.superior Error.estándar
12.33 12.52 12.71 0.09

CONCLUSIÓN