##### UNIVERSIDAD CENTRAL DEL ECUADOR #####
#### AUTOR: Leonardo Ruiz ####
### CARRERA: INGENIERÍA EN PETROLEOS #####
#### VARIABLE RESERVORIO ####
## DATOS ###
# 1. CARGAR DATOS
library(readxl)
datos <- read_excel("C:/Users/LEO/Documents/Producción Campo Sacha.csv.xlsx")
# 2. EXTRAER Salinidad_BH
if(!"Salinidad_BH" %in% colnames(datos)) {
col_Salin <- grep("Salinidad", colnames(datos), ignore.case = TRUE, value = TRUE)[1]
datos$Salinidad_BH <- datos[[col_Salin]]
}
# Extraer variable
Salinidad_BH <- datos$Salinidad_BH
# Convertir a numérico y eliminar NAs
Salinidad_BH <- as.numeric(Salinidad_BH)
ValorNulo <- is.na(Salinidad_BH)
Salinidad_BH <- na.omit(Salinidad_BH)
# 3. CALCULAR INTERVALOS (STURGES - Exacto como tu ejemplo)
n <- length(Salinidad_BH) # Número de datos
k <- 1 + (3.3 * log10(n)) # Fórmula de Sturges
k <- floor(k) # Redondear hacia abajo
min_val <- min(Salinidad_BH) # Mínimo
max_val <- max(Salinidad_BH) # Máximo
R <- max_val - min_val # Rango
A <- R / k # Amplitud
# Redondear valores para mejor presentación
A <- round(A, 2)
min_val <- round(min_val, 2)
max_val <- round(max_val, 2)
# 4. CREAR LÍMITES DE INTERVALOS
Li <- round(seq(from = min_val, to = max_val - A, by = A), 4) # Límites inferiores
Ls <- round(seq(from = min_val + A, to = max_val, by = A), 4) # Límites superiores
MC <- round((Li + Ls) / 2, 2) # Marcas de clase
# 5. CALCULAR FRECUENCIAS ABSOLUTAS (ni)
ni <- numeric(length(Li))
for (i in 1:length(Li)) {
if (i < length(Li)) {
ni[i] <- sum(Salinidad_BH >= Li[i] & Salinidad_BH < Ls[i])
} else {
ni[i] <- sum(Salinidad_BH >= Li[i] & Salinidad_BH <= max_val)
}
}
# 6. CALCULAR TODAS LAS FRECUENCIAS
hi <- ni / sum(ni) * 100 # Frecuencia relativa porcentual
Niasc <- cumsum(ni) # Frecuencia acumulada ascendente
Nidsc <- rev(cumsum(rev(ni))) # Frecuencia acumulada descendente
Hiasc <- cumsum(hi) # Frecuencia relativa acumulada ascendente
Hidsc <- rev(cumsum(rev(hi))) # Frecuencia relativa acumulada descendente
TDF_Salinidad_BH_es <- data.frame(
Li = Li,
Ls = Ls,
MC = MC,
"ni (FA)" = ni,
"hi (FR)" = round(hi, 4),
"Ni (FAAa)" = Niasc,
"Hi (FRAa)" = round(Hiasc, 4),
"Ni (FAAd)" = Nidsc,
"Hi (FRAd)" = round(Hidsc, 4)
)
# 7. MOSTRAR TABLA DE DISTRIBUCIÓN
cat("\n" , rep("=", 80), "\n", sep = "")
##
## ================================================================================
cat("TABLA 1: DISTRIBUCIÓN DE FRECUENCIAS - Salinidad_BH \n")
## TABLA 1: DISTRIBUCIÓN DE FRECUENCIAS - Salinidad_BH
cat(rep("=", 80), "\n\n", sep = "")
## ================================================================================
print(TDF_Salinidad_BH_es)
## Li Ls MC ni..FA. hi..FR. Ni..FAAa. Hi..FRAa. Ni..FAAd. Hi..FRAd.
## 1 0 5080 2540 96 15.0235 96 15.0235 639 100.0000
## 2 5080 10160 7620 162 25.3521 258 40.3756 543 84.9765
## 3 10160 15240 12700 168 26.2911 426 66.6667 381 59.6244
## 4 15240 20320 17780 67 10.4851 493 77.1518 213 33.3333
## 5 20320 25400 22860 25 3.9124 518 81.0642 146 22.8482
## 6 25400 30480 27940 21 3.2864 539 84.3505 121 18.9358
## 7 30480 35560 33020 70 10.9546 609 95.3052 100 15.6495
## 8 35560 40640 38100 13 2.0344 622 97.3396 30 4.6948
## 9 40640 45720 43180 4 0.6260 626 97.9656 17 2.6604
## 10 45720 50800 48260 13 2.0344 639 100.0000 13 2.0344
# =============================================================================
# GRÁFICAS CORTAS: LOCAL vs GLOBAL
# =============================================================================
cat("\n" , rep("=", 60), "\n", sep = "")
##
## ============================================================
cat("GRÁFICAS: LOCAL vs GLOBAL\n")
## GRÁFICAS: LOCAL vs GLOBAL
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# 1. PREPARAR
if("Pozo" %in% colnames(datos)) {
pozos <- unique(datos$Pozo)
local_pozos <- pozos[1:min(9, length(pozos))]
Salinidad_Local <- na.omit(as.numeric(datos$Salinidad_BH[datos$Pozo %in% local_pozos]))
Salinidad_global <- Salinidad_BH
cat("LOCAL:", length(local_pozos), "pozos,", length(Salinidad_Local), "datos\n")
cat("GLOBAL:", length(Salinidad_BH), "datos\n\n")
# 2. CREAR 2 GRÁFICAS
par(mfrow = c(1, 2))
# Gráfica LOCAL
hist(Salinidad_Local, breaks = 10, col = rgb(0.2,0.4,0.8,0.6), border = "white",
main = paste("LOCAL\n", length(local_pozos), " pozos", sep = ""),
xlab = "Salinidad", ylab = "Frecuencia", las = 1)
# Gráfica GLOBAL
hist(Salinidad_global, breaks = k, col = rgb(0.8,0.4,0.2,0.6), border = "white",
main = "GLOBAL\nTodos los pozos",
xlab = "Salinidad", ylab = "Frecuencia", las = 1)
par(mfrow = c(1, 1))
} else {
cat("Solo GLOBAL (no hay columna Pozo)\n\n")
hist(Salinidad_BH, breaks = k, col = "lightblue", border = "white",
main = "Salinidad_BH (GLOBAL)", xlab = "Salinidad", ylab = "Frecuencia", las = 1)
}
## LOCAL: 9 pozos, 71 datos
## GLOBAL: 639 datos
# =============================================================================
# GRÁFICA OJIVA ASCENDENTE
# =============================================================================
cat("\n" , rep("=", 60), "\n", sep = "")
##
## ============================================================
cat("OJIVA ASCENDENTE - Salinidad_BH (GLOBAL)\n")
## OJIVA ASCENDENTE - Salinidad_BH (GLOBAL)
cat(rep("=", 60), "\n\n", sep = "")
## ============================================================
# Crear puntos para la ojiva
puntos_x <- c(Ls) # Los límites superiores
puntos_y <- cumsum(ni) # Frecuencia acumulada
# 1. Configurar márgenes para espacio a la derecha
par(mar = c(5, 4, 4, 8)) # Más espacio a la derecha (8)
# 2. Crear ojiva
plot(puntos_x, puntos_y,
type = "b",
pch = 16,
col = "darkblue",
lwd = 2,
main = "OJIVA ASCENDENTE\ Salinidad - GLOBAL",
xlab = "Límite Superior del Intervalo (Salinidad)",
ylab = "Frecuencia Acumulada (Ni)",
las = 1,
xlim = c(min_val, max_val),
ylim = c(0, n * 1.05))
# 3. Grid
grid(col = "gray80", lty = "dotted")
# 4. Eje derecho con porcentajes
porcentajes <- c(0, 25, 50, 75, 100)
axis(4,
at = porcentajes/100 * n,
labels = paste0(porcentajes, "%"),
las = 1,
col.axis = "darkred",
cex.axis = 0.9)
mtext("Porcentaje Acumulado", side = 4, line = 3, col = "darkred", cex = 0.9)
# 5. Puntos importantes (50% y 75%)
# 50%
indice_50 <- which.max(Hiasc >= 50)
if(length(indice_50) > 0) {
x50 <- puntos_x[indice_50 + 1]
y50 <- puntos_y[indice_50 + 1]
points(x50, y50, col = "red", pch = 17, cex = 1.3) # Triángulo
segments(x50, 0, x50, y50, col = "red", lty = 2, lwd = 1)
segments(min_val, y50, x50, y50, col = "red", lty = 2, lwd = 1)
}
# 75%
indice_75 <- which.max(Hiasc >= 75)
if(length(indice_75) > 0) {
x75 <- puntos_x[indice_75 + 1]
y75 <- puntos_y[indice_75 + 1]
points(x75, y75, col = "darkgreen", pch = 15, cex = 1.3) # Cuadrado
segments(x75, 0, x75, y75, col = "darkgreen", lty = 2, lwd = 1)
segments(min_val, y75, x75, y75, col = "darkgreen", lty = 2, lwd = 1)
}
# =============================================================================
# DIAGRAMA DE CAJA
# =============================================================================
cat("\n" , rep("=", 60), "\n", sep = "")
##
## ============================================================
cat("4.2 Diagrama de caja\n")
## 4.2 Diagrama de caja
cat(rep("=", 60), "\n\n", sep = "")
## ============================================================
# Exactamente como el ejemplo
boxplot(Salinidad_BH,
horizontal = TRUE,
col = "blue",
main = "Gráfica: Distribución de Salinidad (Salinidad_BH)",
xlab = "Salinidad")
cat("ESTADÍSTICAS:\n")
## ESTADÍSTICAS:
s <- boxplot.stats(Salinidad_BH)$stats
cat("• Mín:", round(s[1], 2), "Salinidad\n")
## • Mín: 0 Salinidad
cat("• Q1: ", round(s[2], 2), "Salinidad\n")
## • Q1: 6250 Salinidad
cat("• Med:", round(s[3], 2), "Salinidad\n")
## • Med: 11200 Salinidad
cat("• Q3: ", round(s[4], 2), "Salinidad\n")
## • Q3: 16300 Salinidad
cat("• Máx:", round(s[5], 2), "Salinidad\n")
## • Máx: 28500 Salinidad
# =============================================================================
# 5. INDICADORES ESTADÍSTICOS - Salinidad_BH
# =============================================================================
cat("\n" , rep("=", 70), "\n", sep = "")
##
## ======================================================================
cat("5. INDICADORES ESTADÍSTICOS\n")
## 5. INDICADORES ESTADÍSTICOS
cat(rep("=", 70), "\n\n", sep = "")
## ======================================================================
# Cargar librería para asimetría y curtosis
if(!require(e1071)) {
install.packages("e1071")
library(e1071)
}
## Cargando paquete requerido: e1071
# =============================================================================
# 5.1 MEDIDAS DE POSICIÓN
# =============================================================================
cat("5.1 MEDIDAS DE POSICIÓN:\n")
## 5.1 MEDIDAS DE POSICIÓN:
cat("------------------------\n")
## ------------------------
# Media aritmética
x <- mean(Salinidad_BH)
cat("Media aritmética (x̄):", round(x, 4), "Salinidad\n")
## Media aritmética (x̄): 14378.23 Salinidad
# Mediana (de los datos originales)
Me <- median(Salinidad_BH)
cat("Mediana (Me): ", round(Me, 4), "Salinidad\n")
## Mediana (Me): 11200 Salinidad
# Mediana de las marcas de clase (opcional)
Me_MC <- median(MC)
cat("Mediana de MC: ", round(Me_MC, 4), "Salinidad\n")
## Mediana de MC: 25400 Salinidad
# Mínimo y máximo
ri <- min(Salinidad_BH)
rs <- max(Salinidad_BH)
cat("Mínimo (ri): ", round(ri, 4), "Salinidad\n")
## Mínimo (ri): 0 Salinidad
cat("Máximo (rs): ", round(rs, 4), "Salinidad\n")
## Máximo (rs): 50800 Salinidad
# Moda (aproximada)
calcular_moda <- function(x) {
freq <- table(round(x, 1)) # Redondear a 1 decimal
moda_val <- as.numeric(names(freq)[which.max(freq)])
return(moda_val)
}
Mo <- calcular_moda(Salinidad_BH)
cat("Moda aproximada (Mo): ", round(Mo, 4), "°API\n")
## Moda aproximada (Mo): 6250 °API
# Cuartiles
Q1 <- quantile(Salinidad_BH, 0.25)
Q2 <- quantile(Salinidad_BH, 0.50) # Igual que mediana
Q3 <- quantile(Salinidad_BH, 0.75)
cat("\nCuartiles:\n")
##
## Cuartiles:
cat(" Q1 (25%): ", round(Q1, 4), "Salinidad\n")
## Q1 (25%): 6250 Salinidad
cat(" Q2 (50% - Mediana): ", round(Q2, 4), "Salinidad\n")
## Q2 (50% - Mediana): 11200 Salinidad
cat(" Q3 (75%): ", round(Q3, 4), "Salinidad\n")
## Q3 (75%): 16300 Salinidad
# Percentiles
percentiles <- quantile(Salinidad_BH, probs = c(0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99))
cat("\nPercentiles seleccionados:\n")
##
## Percentiles seleccionados:
for(i in 1:length(percentiles)) {
cat(" P", names(percentiles)[i], ": ",
round(percentiles[i], 4), "Salinidad\n", sep = "")
}
## P10%: 3200Salinidad
## P25%: 6250Salinidad
## P50%: 11200Salinidad
## P75%: 16300Salinidad
## P90%: 31900Salinidad
## P95%: 35000Salinidad
## P99%: 50800Salinidad
# =============================================================================
# 5.2 MEDIDAS DE DISPERSIÓN
# =============================================================================
cat("\n\n5.2 MEDIDAS DE DISPERSIÓN:\n")
##
##
## 5.2 MEDIDAS DE DISPERSIÓN:
cat("---------------------------\n")
## ---------------------------
# Rango
R <- rs - ri
cat("Rango (R): ", round(R, 4), "Salinidad\n")
## Rango (R): 50800 Salinidad
# Rango intercuartílico
RIQ <- Q3 - Q1
cat("Rango intercuartílico:", round(RIQ, 4), "Salinidad\n")
## Rango intercuartílico: 10050 Salinidad
# Varianza
varianza <- var(Salinidad_BH)
cat("Varianza (s²): ", round(varianza, 4), "\n")
## Varianza (s²): 127472091
# Desviación estándar
s <- sd(Salinidad_BH)
cat("Desviación estándar (s):", round(s, 4), "Salinidad\n")
## Desviación estándar (s): 11290.35 Salinidad
# Coeficiente de variación
CV <- (s / x) * 100
cat("Coeficiente de variación (CV):", round(CV, 2), "%\n")
## Coeficiente de variación (CV): 78.52 %
# Error estándar de la media
EE <- s / sqrt(n)
cat("Error estándar de la media: ", round(EE, 4), "Salinidad\n")
## Error estándar de la media: 446.6395 Salinidad
# =============================================================================
# 5.3 MEDIDAS DE FORMA
# =============================================================================
cat("\n\n5.3 MEDIDAS DE FORMA:\n")
##
##
## 5.3 MEDIDAS DE FORMA:
cat("---------------------\n")
## ---------------------
# Coeficiente de asimetría (skewness)
As <- skewness(Salinidad_BH)
cat("Coeficiente de asimetría (As):", round(As, 4), "\n")
## Coeficiente de asimetría (As): 1.3049
# Interpretación de asimetría
if(abs(As) < 0.5) {
cat(" Interpretación: Distribución aproximadamente simétrica\n")
} else if(As > 0) {
cat(" Interpretación: Asimetría positiva (cola a la derecha)\n")
} else {
cat(" Interpretación: Asimetría negativa (cola a la izquierda)\n")
}
## Interpretación: Asimetría positiva (cola a la derecha)
# Coeficiente de curtosis (kurtosis)
k <- kurtosis(Salinidad_BH)
cat("\nCoeficiente de curtosis (k): ", round(k, 4), "\n")
##
## Coeficiente de curtosis (k): 1.0662
# Curtosis exceso (comparada con normal = 3)
k_exceso <- k - 3
cat("Curtosis exceso (k-3): ", round(k_exceso, 4), "\n")
## Curtosis exceso (k-3): -1.9338
# Interpretación de curtosis
if(abs(k_exceso) < 0.5) {
cat(" Interpretación: Mesocúrtica (similar a normal)\n")
} else if(k_exceso > 0) {
cat(" Interpretación: Leptocúrtica (más picuda que normal)\n")
} else {
cat(" Interpretación: Platicúrtica (más aplanada que normal)\n")
}
## Interpretación: Platicúrtica (más aplanada que normal)
# =============================================================================
# 5.4 RESUMEN EN TABLA
# =============================================================================
cat("\n\n" , rep("-", 70), "\n", sep = "")
##
##
## ----------------------------------------------------------------------
cat("RESUMEN DE INDICADORES ESTADÍSTICOS\n")
## RESUMEN DE INDICADORES ESTADÍSTICOS
cat(rep("-", 70), "\n\n", sep = "")
## ----------------------------------------------------------------------
# Crear tabla resumen
resumen_estadisticas <- data.frame(
"Categoría" = c(
rep("POSICIÓN", 8),
rep("DISPERSIÓN", 6),
rep("FORMA", 3)
),
"Indicador" = c(
"n (tamaño muestra)",
"Media aritmética",
"Mediana",
"Moda aproximada",
"Mínimo",
"Máximo",
"Q1 (25%)",
"Q3 (75%)",
"Rango",
"Rango intercuartílico",
"Varianza",
"Desviación estándar",
"Coeficiente variación",
"Error estándar",
"Asimetría",
"Curtosis",
"Curtosis exceso"
),
"Valor" = c(
n,
round(x, 4),
round(Me, 4),
round(Mo, 4),
round(ri, 4),
round(rs, 4),
round(Q1, 4),
round(Q3, 4),
round(R, 4),
round(RIQ, 4),
round(varianza, 4),
round(s, 4),
paste0(round(CV, 2), "%"),
round(EE, 4),
round(As, 4),
round(k, 4),
round(k_exceso, 4)
),
"Unidad" = c(
"observaciones",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"Salinidad",
"(Salinidad)²",
"Salinidad",
"%",
"Salinidad",
"adimensional",
"adimensional",
"adimensional"
)
)
print(resumen_estadisticas)
## Categoría Indicador Valor Unidad
## 1 POSICIÓN n (tamaño muestra) 639 observaciones
## 2 POSICIÓN Media aritmética 14378.2349 Salinidad
## 3 POSICIÓN Mediana 11200 Salinidad
## 4 POSICIÓN Moda aproximada 6250 Salinidad
## 5 POSICIÓN Mínimo 0 Salinidad
## 6 POSICIÓN Máximo 50800 Salinidad
## 7 POSICIÓN Q1 (25%) 6250 Salinidad
## 8 POSICIÓN Q3 (75%) 16300 Salinidad
## 9 DISPERSIÓN Rango 50800 Salinidad
## 10 DISPERSIÓN Rango intercuartílico 10050 Salinidad
## 11 DISPERSIÓN Varianza 127472090.8325 (Salinidad)²
## 12 DISPERSIÓN Desviación estándar 11290.3539 Salinidad
## 13 DISPERSIÓN Coeficiente variación 78.52% %
## 14 DISPERSIÓN Error estándar 446.6395 Salinidad
## 15 FORMA Asimetría 1.3049 adimensional
## 16 FORMA Curtosis 1.0662 adimensional
## 17 FORMA Curtosis exceso -1.9338 adimensional
# =============================================================================
# 6. TABLA RESUMEN FINAL
# =============================================================================
cat("\n" , rep("=", 70), "\n", sep = "")
##
## ======================================================================
cat("6. TABLA RESUMEN\n")
## 6. TABLA RESUMEN
cat(rep("=", 70), "\n\n", sep = "")
## ======================================================================
# Ya tienes resumen_estadisticas, solo decir que es la tabla resumen
cat("La tabla anterior (Resumen de Indicadores Estadísticos) es la tabla resumen.\n")
## La tabla anterior (Resumen de Indicadores Estadísticos) es la tabla resumen.
cat("Contiene todas las medidas calculadas organizadas por categoría.\n")
## Contiene todas las medidas calculadas organizadas por categoría.
# =============================================================================
# 7. CONCLUSIÓN
# =============================================================================
cat("\n" , rep("=", 70), "\n", sep = "")
##
## ======================================================================
cat("7. CONCLUSIÓN\n")
## 7. CONCLUSIÓN
cat(rep("=", 70), "\n\n", sep = "")
## ======================================================================
cat("ANÁLISIS DE SALINIDAD_BH EN EL CAMPO SACHA:\n\n")
## ANÁLISIS DE SALINIDAD_BH EN EL CAMPO SACHA:
cat("1. El análisis de", n, "observaciones de Salinidad_BH muestra:\n")
## 1. El análisis de 639 observaciones de Salinidad_BH muestra:
cat(" • Promedio:", round(x, 1), "Salinidad ")
## • Promedio: 14378.2 Salinidad
# CORRECCIÓN: Todo en una línea o con llaves
if(x < 10) {
cat("Agua Dulce/Baja Salinidad")
} else if(x < 22) {
cat("Agua Salobre/Salinidad Media")
} else if(x < 31) {
cat("Agua Salada/Alta Salinidad")
} else {
cat("Brina/Hipersalina")
}
## Brina/Hipersalina
cat(")\n")
## )
cat(" • Variabilidad:", round(CV, 1), "% (")
## • Variabilidad: 78.5 % (
# Misma corrección aquí
if(CV > 30) {
cat("alta")
} else if(CV > 15) {
cat("moderada")
} else {
cat("baja")
}
## alta
cat(")\n")
## )
cat(" • Distribución: ")
## • Distribución:
# Y aquí
if(abs(As) < 0.5) {
cat("Cercana a simétrica")
} else if(As > 0) {
cat("sesgada a la derecha (valores altos)")
} else {
cat("sesgada a la izquierda (valores bajos)")
}
## sesgada a la derecha (valores altos)
cat("\n\n")
cat("2. El análisis estadístico completo proporciona una caracterización\n")
## 2. El análisis estadístico completo proporciona una caracterización
cat(" detallada de la calidad del petróleo en el campo Sacha.\n")
## detallada de la calidad del petróleo en el campo Sacha.
cat("\n✓ Análisis estadístico completado exitosamente\n")
##
## ✓ Análisis estadístico completado exitosamente