Análisis Estadístico MEDIAN

CARGA DE DATOS Y LIBRERÍAS

library(moments)
library(MASS)
library(gt)
library(dplyr)

## 
## Adjuntando el paquete: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

datos <- read.csv("~/UNI/estadistica/2026/dataset_geologico_limpio_80.csv")
cat("Dimensiones:", nrow(datos), "observaciones,", ncol(datos), "variables\n\n")

## Dimensiones: 27784 observaciones, 58 variables

#Limpiar la variable
med_raw <- datos$MEDIAN
med <- as.numeric(as.character(med_raw))
cat("Primeros valores de MEDIANA:\n")

## Primeros valores de MEDIANA:

print(head(med))

## [1] -0.29  0.43  1.75  1.70  0.97  3.27

# Limpiar datos
med <- na.omit(med)
med <- med[med > 0]
n <- length(med)

cat("Número de datos válidos:", length(med), "\n\n")

## Número de datos válidos: 25955

GRAFICA Y TABLA DE DISTRIBUCION DE FRECUENCIA

# Parámetros
minimo <- min(med)
maximo <- max(med)
R <- maximo - minimo
k <- floor(1 + 3.3 * log10(n))
A <- R / k

cat("Número de intervalos:", k, "\n")

## Número de intervalos: 15

cat("Amplitud:", A, "\n\n")

## Amplitud: 0.8816506

# Intervalos
li <- seq(minimo, maximo - A, by = A)
ls <- seq(minimo + A, maximo, by = A)
MC <- (li + ls) / 2

# Frecuencias absolutas
ni <- numeric(length(li))

for (i in 1:length(li)) {
  ni[i] <- sum(med >= li[i] & med < ls[i])
}

# Último intervalo incluye extremo superior
ni[length(li)] <- sum(med >= li[length(li)] & med <= maximo)

# Frecuencia relativa
hi <- (ni / sum(ni)) * 100

# Tabla final
tabla_MEDIAN <- data.frame(
  Li = round(li, 2),
  Ls= round(ls, 2),
  MC = round(MC, 2),
  Ni = ni,
  Hi = round(hi, 2)
)

# Agregar TOTAL
tabla_MEDIAN[nrow(tabla_MEDIAN) + 1, ] <- 
  c(NA, NA, "TOTAL", sum(ni), sum(hi))

tabla_MEDIAN

##       Li    Ls    MC    Ni    Hi
## 1      0  0.88  0.44  2500  9.63
## 2   0.88  1.76  1.32  3938 15.17
## 3   1.76  2.65  2.21  4498 17.33
## 4   2.65  3.53  3.09  4040 15.57
## 5   3.53  4.41  3.97  2119  8.16
## 6   4.41  5.29  4.85  1526  5.88
## 7   5.29  6.17  5.73  1509  5.81
## 8   6.17  7.05  6.61  2042  7.87
## 9   7.05  7.94   7.5  2188  8.43
## 10  7.94  8.82  8.38  1099  4.23
## 11  8.82   9.7  9.26   443  1.71
## 12   9.7 10.58 10.14    45  0.17
## 13 10.58 11.46 11.02     4  0.02
## 14 11.46 12.34  11.9     3  0.01
## 15 12.34 13.23 12.79     1     0
## 16  <NA>  <NA> TOTAL 25955   100

tabla_MEDIAN %>%
  gt() %>%
  tab_header(
    title = md("**Tabla Nº: Distribución de frecuencias de la variable MEDIAN**"),
    subtitle = md("**Tamaño de grano en depósitos marinos**")
  ) %>%
  tab_source_note(
    source_note = md("**Autor: Grupo 3**")
  ) %>%
  cols_label(
    Li = "Límite inferior",
    Ls = "Límite superior",
    MC = "Marca de clase",
    Ni = "Frecuencia (ni)",
    Hi = "Frecuencia relativa (%)"
  ) %>%
  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"
  ) %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(
      rows = MC == "TOTAL"
    )
  )

Tabla Nº: Distribución de frecuencias de la variable MEDIAN
Tamaño de grano en depósitos marinos
Límite inferior	Límite superior	Marca de clase	Frecuencia (ni)	Frecuencia relativa (%)
0	0.88	0.44	2500	9.63
0.88	1.76	1.32	3938	15.17
1.76	2.65	2.21	4498	17.33
2.65	3.53	3.09	4040	15.57
3.53	4.41	3.97	2119	8.16
4.41	5.29	4.85	1526	5.88
5.29	6.17	5.73	1509	5.81
6.17	7.05	6.61	2042	7.87
7.05	7.94	7.5	2188	8.43
7.94	8.82	8.38	1099	4.23
8.82	9.7	9.26	443	1.71
9.7	10.58	10.14	45	0.17
10.58	11.46	11.02	4	0.02
11.46	12.34	11.9	3	0.01
12.34	13.23	12.79	1	0
NA	NA	TOTAL	25955	100
Autor: Grupo 3

# Cortes exactos de Sturges
cortes <- c(li, maximo)
par(mfrow = c(1,1))

hist(med,
     probability = TRUE,
     breaks = cortes,
     col = "lightblue",
     border = "white",
     main = "Distribución original del tamaño de grano",
     xlab = "MEDIAN",
     ylab = "Densidad",
     xlim = range(cortes),
     ylim = c(0, 0.30))

CONJETURA DEL MODELO 1

Debido a la similitud de las barras asociamos con el modelo de probabilidad normal al tramo 1

par(mar = c(5,5,5,2)) 
# Punto de corte en la barra 6
corte_6 <- ls[6]

# Datos de las primeras 6 clases
med_normal <- med[med <= corte_6]

# Cortes de las primeras 6 barras
cortes_6 <- c(li[1:6], ls[6])

# Parámetros de la normal
media_n <- mean(med_normal)
sd_n <- sd(med_normal)

hist(med_normal,
     probability = TRUE,
     breaks = cortes_6,
     col = "lightgreen",
     border = "white",
     main = "Primeras 6 barras con ajuste Normal",
     xlab = "MEDIAN",
     ylab = "Densidad",
     xlim = range(cortes_6),
     ylim = c(0, 0.35))

axis(1,
     at = seq(floor(min(med_normal)),
              ceiling(max(med_normal)),
              by = 1))

curve(dnorm(x,
            mean = media_n,
            sd = sd_n),
      add = TRUE,
      col = "blue",
      lwd = 2)

legend("topright",
       legend = c("Normal"),
       col = c("blue"),
       lwd = 2,
       bty = "n")

TEST DE APROBACIÓN

# Histograma sin graficar
Histograma_1 <- hist(med_normal,
                     breaks = cortes_6,
                     plot = FALSE)
# Frecuencias observadas
Fo_1 <- Histograma_1$counts
n1 <- sum(Fo_1)

# Frecuencias esperadas (modelo normal)
Fe_1 <- numeric(length(Fo_1))

for (i in 1:length(Fo_1)) {
  Fe_1[i] <- n1 * (pnorm(ls[i], mean = media_n, sd = sd_n) -
                     pnorm(li[i], mean = media_n, sd = sd_n))
}

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

# Convertir a porcentaje
Fo_1 <- (Fo_1 / n1) * 100
Fe_1 <- (Fe_1 / n1) * 100

# Gráfico correlación
plot(Fo_1, Fe_1,
     main = "Gráfica Nº: Correlación de frecuencias - Modelo Normal",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     col = "blue3")

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

# Correlación
Correlacion_1 <- cor(Fo_1, Fe_1) * 100
Correlacion_1

## [1] 94.87911

cat("APRUEBA TEST DE PEARSON\n\n")

## APRUEBA TEST DE PEARSON

# ================================
# TEST CHI-CUADRADO
# ================================

grados_libertad_1 <- length(Fo_1) - 1
nivel_significancia <- 0.95

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

umbral_aceptacion_1 <- qchisq(nivel_significancia, grados_libertad_1)

cat("Chi-cuadrado:", x2_1, "\n")

## Chi-cuadrado: 6.500696

cat("Umbral:", umbral_aceptacion_1, "\n")

## Umbral: 11.0705

if(x2_1 < umbral_aceptacion_1){
  cat("APRUEBA TEST DE CHI-CUADRADO\n")
} else {
  cat("NO APRUEBA TEST DE CHI-CUADRADO\n")
}

## APRUEBA TEST DE CHI-CUADRADO

# ================================
# TABLA RESUMEN
# ================================

Variable <- c("MEDIAN (Normal)")
tabla_resumen <- data.frame(
  Variable,
  round(Correlacion_1, 2),
  round(x2_1, 2),
  round(umbral_aceptacion_1, 2)
)

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

tabla_resumen

##          Variable Test Pearson (%) Chi Cuadrado Umbral
## 1 MEDIAN (Normal)            94.88          6.5  11.07

tabla_resumen %>%
  gt() %>%
  tab_header(
    title = md("**Tabla Nº: Test de bondad de ajuste (Modelo Normal)**"),
    subtitle = md("**Variable MEDIAN**")
  ) %>%
  tab_source_note(
    source_note = md("**Autor: Grupo 3**")
  ) %>%
  cols_label(
    Variable = "Variable",
    `Test Pearson (%)` = "Test Pearson (%)",
    `Chi Cuadrado` = "Chi Cuadrado",
    Umbral = "Umbral de aceptación"
  ) %>%
  fmt_number(
    columns = c(`Test Pearson (%)`, `Chi Cuadrado`, Umbral),
    decimals = 2
  ) %>%
  cols_align(
    align = "center",
    everything()
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray"
  )

Tabla Nº: Test de bondad de ajuste (Modelo Normal)
Variable MEDIAN
Variable	Test Pearson (%)	Chi Cuadrado	Umbral de aceptación
MEDIAN (Normal)	94.88	6.50	11.07
Autor: Grupo 3

CÁLCULO DE PROBABILIDADES

¿Cuál es la probabilidad de que MEDIAN esté entre 2 y 4?

# Parámetros del modelo normal ya calculados antes
u_1 <- media_n
sigma_1 <- sd_n
med_1 <- med_normal

# Cálculo 
Probabilidad_1 <- (pnorm(4, u_1, sigma_1) - pnorm(2, u_1, sigma_1)) * 100
Probabilidad_1

## [1] 52.22483

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

plot(x, dnorm(x, u_1, sigma_1), 
     col = "skyblue3",
     lwd = 2,
     main = "Gráfica N°: Cálculo de probabilidades del modelo normal",
     ylab = "Densidad de probabilidad",
     xlab = "MEDIAN")

# Definir el rango de la sección que quieres pintar
x_section <- seq(2, 4, 0.001)
y_section <- dnorm(x_section, u_1, sigma_1)

# Pintar la sección de la curva
lines(x_section, y_section, col = "red", lwd = 2)

# Pintar el área debajo de la línea roja
polygon(c(x_section, rev(x_section)), 
        c(y_section, rep(0, length(y_section))),
        col = rgb(1, 0, 0, 0.6),
        border = NA)

¿De 300 observaciones, cuántas tendrían MEDIAN entre 2 y 4?

cantidad_1 <- (pnorm(4, u_1, sigma_1) - pnorm(2, u_1, sigma_1)) * 300
cantidad_1

## [1] 156.6745

CONJETURA DEL MODELO 2

Debido a la similitud de las barras asociamos con el modelo de probabilidad lognormal al tramo 2

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

# Datos de las barras restantes
med_lognorm <- med[med > corte_6]

# Cortes de las barras restantes
li_ln <- li[7:length(li)]
ls_ln <- ls[7:length(ls)]
cortes_ln <- c(li_ln, maximo)

# Ajuste log-normal
fit_ln <- fitdistr(med_lognorm, "lognormal")
meanlog_ln <- fit_ln$estimate["meanlog"]
sdlog_ln <- fit_ln$estimate["sdlog"]

# Histograma
hist(med_lognorm,
     probability = TRUE,
     breaks = cortes_ln,
     col = "plum",
     border = "white",
     main = "Barras restantes con ajuste Log-normal",
     xlab = "MEDIAN",
     ylab = "Densidad",
     xlim = range(cortes_ln),
     ylim = c(0, 0.4),
     xaxt = "n")

axis(1,
     at = seq(floor(min(med_lognorm)),
              ceiling(max(med_lognorm)),
              by = 1))

curve(dlnorm(x,
             meanlog = meanlog_ln,
             sdlog = sdlog_ln),
      add = TRUE,
      col = "red",
      lwd = 2)

legend("topright",
       legend = c("Log-normal"),
       col = c("red"),
       lwd = 2,
       bty = "n")

TEST DE APROBACIÓN

# Histograma sin graficar
Histograma_2 <- hist(med_lognorm,
                     breaks = cortes_ln,
                     plot = FALSE)

# Frecuencias observadas
Fo_2 <- Histograma_2$counts
n2 <- sum(Fo_2)

# Frecuencias esperadas (modelo log-normal)
Fe_2 <- numeric(length(Fo_2))

for (i in 1:length(Fo_2)) {
  Fe_2[i] <- n2 * (plnorm(ls_ln[i],
                          meanlog = meanlog_ln,
                          sdlog = sdlog_ln) -
                     plnorm(li_ln[i],
                            meanlog = meanlog_ln,
                            sdlog = sdlog_ln))
}

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

# Convertir a porcentaje
Fo_2 <- (Fo_2 / n2) * 100
Fe_2 <- (Fe_2 / n2) * 100

plot(Fo_2, Fe_2,
     main = "Gráfica Nº: Correlación de frecuencias - Modelo Log-normal",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     col = "purple3")

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

Correlacion_2 <- cor(Fo_2, Fe_2) * 100
Correlacion_2

## [1] 98.25249

cat("APRUEBA TEST DE PEARSON\n\n")

## APRUEBA TEST DE PEARSON

# ================================
# TEST CHI-CUADRADO
# ================================

grados_libertad_2 <- length(Fo_2) - 1
nivel_significancia <- 0.95

x2_2 <- sum((Fe_2 - Fo_2)^2 / Fe_2)

umbral_aceptacion_2 <- qchisq(nivel_significancia, grados_libertad_2)

cat("Chi-cuadrado:", x2_2, "\n")

## Chi-cuadrado: 2.910801

cat("Umbral:", umbral_aceptacion_2, "\n")

## Umbral: 15.50731

if(x2_2 < umbral_aceptacion_2){
  cat("APRUEBA TEST DE CHI-CUADRADO\n")
} else {
  cat("NO APRUEBA TEST DE CHI-CUADRADO\n")
}

## APRUEBA TEST DE CHI-CUADRADO

# ================================
# TABLA RESUMEN
# ================================

Variable <- c("MEDIAN (Log-normal)")
tabla_resumen_2 <- data.frame(
  Variable,
  round(Correlacion_2, 2),
  round(x2_2, 2),
  round(umbral_aceptacion_2, 2)
)

colnames(tabla_resumen_2) <- c(
  "Variable",
  "Test Pearson (%)",
  "Chi Cuadrado",
  "Umbral"
)

tabla_resumen_2

##              Variable Test Pearson (%) Chi Cuadrado Umbral
## 1 MEDIAN (Log-normal)            98.25         2.91  15.51

tabla_resumen_2 %>%
  gt() %>%
  tab_header(
    title = md("**Tabla Nº: Test de bondad de ajuste (Modelo Log-normal)**"),
    subtitle = md("**Variable MEDIAN**")
  ) %>%
  tab_source_note(
    source_note = md("**Autor: Grupo 3**")
  ) %>%
  cols_label(
    Variable = "Variable",
    `Test Pearson (%)` = "Test Pearson (%)",
    `Chi Cuadrado` = "Chi Cuadrado",
    Umbral = "Umbral de aceptación"
  ) %>%
  fmt_number(
    columns = c(`Test Pearson (%)`, `Chi Cuadrado`, Umbral),
    decimals = 2
  ) %>%
  cols_align(
    align = "center",
    everything()
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray"
  )

Tabla Nº: Test de bondad de ajuste (Modelo Log-normal)
Variable MEDIAN
Variable	Test Pearson (%)	Chi Cuadrado	Umbral de aceptación
MEDIAN (Log-normal)	98.25	2.91	15.51
Autor: Grupo 3

CÁLCULO DE PROBABILIDADES

¿Cuál es la probabilidad de que MEDIAN esté entre 6 y 9?

# Parámetros del modelo log-normal ya calculados antes
med_2 <- med_lognorm
meanlog_2 <- meanlog_ln
sdlog_2 <- sdlog_ln

# Cálculos
Probabilidad_2 <- (plnorm(9, meanlog_2, sdlog_2) - plnorm(6, meanlog_2, sdlog_2)) * 100
Probabilidad_2

## [1] 81.99887

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

plot(x, dlnorm(x, meanlog_2, sdlog_2), 
     col = "darkorchid3",
     lwd = 2,
     main = "Gráfica N°: Cálculo de probabilidades del modelo log-normal",
     ylab = "Densidad de probabilidad",
     xlab = "MEDIAN")

# Definir el rango de la sección que quieres pintar
x_section <- seq(6, 9, 0.001)
y_section <- dlnorm(x_section, meanlog_2, sdlog_2)

# Pintar la sección de la curva
lines(x_section, y_section, col = "red", lwd = 2)

# Pintar el área debajo de la línea roja
polygon(c(x_section, rev(x_section)), 
        c(y_section, rep(0, length(y_section))),
        col = rgb(1, 0, 0, 0.6),
        border = NA)

¿De 300 observaciones, cuántas tendrían MEDIAN entre 6 y 9?

# Cálculos 
cantidad_2 <- (plnorm(9, meanlog_2, sdlog_2) - plnorm(6, meanlog_2, sdlog_2)) * 300
cantidad_2

## [1] 245.9966

INTERVALOS DE CONFIANZA

# Media aritmética
x_barra <- mean(med)
x_barra

## [1] 3.746772

# Desviación estándar
sigma <- sd(med)
sigma

## [1] 2.442155

# Tamaño muestral
n <- length(med)
n

## [1] 25955

# Error estándar
e <- sigma / sqrt(n)
e

## [1] 0.01515873

# Intervalo de confianza aproximado al 95%
li_ic <- x_barra - 2 * e
li_ic

## [1] 3.716454

ls_ic <- x_barra + 2 * e
ls_ic

## [1] 3.777089

# Tabla resumen
tabla_media <- data.frame(
  Variable = "MEDIAN",
  Limite_inferior = round(li_ic, 2),
  Limite_superior = round(ls_ic, 2),
  Desviación_estandar_poblacional = round(e, 6)
)

tabla_media

##   Variable Limite_inferior Limite_superior Desviación_estandar_poblacional
## 1   MEDIAN            3.72            3.78                        0.015159

tabla_media %>%
  gt() %>%
  tab_header(
    title = md("**Tabla Nº: Intervalo de confianza de la media**"),
    subtitle = md("**Variable MEDIAN (95% de confianza)**")
  ) %>%
  tab_source_note(
    source_note = md("**Autor: Grupo 3**")
  ) %>%
  cols_label(
    Variable = "Variable",
    Limite_inferior = "Límite inferior",
    Limite_superior = "Límite superior",
    Desviación_estandar_poblacional = "Desviación estándar poblacional"
  ) %>%
  fmt_number(
    columns = c(Limite_inferior, Limite_superior, Desviación_estandar_poblacional),
    decimals = 2
  ) %>%
  cols_align(
    align = "center",
    everything()
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray"
  )

Tabla Nº: Intervalo de confianza de la media
Variable MEDIAN (95% de confianza)
Variable	Límite inferior	Límite superior	Desviación estándar poblacional
MEDIAN	3.72	3.78	0.02
Autor: Grupo 3

CONCLUSIONES

La variable MEDIAN se explica mediante un modelo mixto en dos grupos: un modelo normal en las primeras 6 barras y un modelo log-normal en las barras restantes. La media aritmética de la variable es 3.75 y su desviación estándar es 2.44

De esta manera, se lograron calcular probabilidades. Por ejemplo, bajo el modelo normal, la probabilidad de que la variable MEDIAN se ubique entre 2 y 4 es de “52.22%”; mientras que, bajo el modelo log-normal, la probabilidad de que MEDIAN se encuentre entre 6 y 9 es de 82%.

Mediante el teorema del límite central, se establece que la media poblacional de la variable MEDIAN se encuentra entre 3.72 y 3.78 con un 95% de confianza.**