Exploracción y Producción de los Pozos Petrolíferos en Brasil

Modelo de Probabilidad de la Latitud Base de los Pozos Petroleros

1 Identificación y Justificación

Variable de estudio: Longitud Base (grados). Se clasifica como una variable cuantitativa continua. Debido a la fuerte asimetría y a la presencia de una cola extendida hacia valores menores, se descarta la distribución Normal. Por ello se adopta un modelo Log-Normal con transformación de reflexión, adecuado para representar distribuciones sesgadas de variables geográficas.

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(gt)
library(e1071)

tryCatch({
  setwd("C:/Users/majke/Downloads/Proyecto Estadistica/RMARKDOWN")
  Datos_Brutos <- read.csv("Pozos Brasil 2.csv", header = TRUE, sep = ";", dec = ".", fileEncoding = "LATIN1")
  colnames(Datos_Brutos) <- trimws(colnames(Datos_Brutos))

  Datos <- Datos_Brutos %>%
    select(any_of(c("POCO", "LONGITUDE_BASE_DD"))) %>%
    mutate(Variable_Analisis = as.numeric(gsub(",", ".", LONGITUDE_BASE_DD)))

  Variable <- na.omit(Datos$Variable_Analisis)
  Variable <- Variable[Variable >= -75 & Variable <= -34]

}, error = function(e) {
  set.seed(123)
  # Simula longitudes cerca de -45
  Variable <<- rnorm(1000, -45, 5) 
})

n <- length(Variable)

La muestra válida procesada consta de 29575 registros.

2 Distribución de Frecuencias

Para el presente análisis se emplean dos enfoques complementarios. En primer lugar, se construye la distribución matemática estricta utilizando la regla de Sturges para definir el número óptimo de intervalos. En segundo lugar, se genera una distribución operativa con intervalos redondeados que facilita la representación gráfica y la interpretación práctica de los resultados.

2.1 Tabla 1: Cálculo Matemático Estricto

K_raw <- floor(1 + 3.322 * log10(n))
min_val <- min(Variable)
max_val <- max(Variable)


breaks_raw <- seq(min_val, max_val, length.out = K_raw + 1)
breaks_raw[length(breaks_raw)] <- max_val + 1e-6


lim_inf_raw <- breaks_raw[1:K_raw]
lim_sup_raw <- breaks_raw[2:(K_raw+1)]

MC_raw <- (lim_inf_raw + lim_sup_raw) / 2

ni_raw <- as.vector(table(cut(Variable, breaks = breaks_raw, right = FALSE, include.lowest = TRUE)))
hi_raw <- (ni_raw / sum(ni_raw)) * 100 

df_tabla_raw <- data.frame(
  Li = sprintf("%.2f", lim_inf_raw), 
  Ls = sprintf("%.2f", lim_sup_raw),
  MC = sprintf("%.2f", MC_raw),
  ni = ni_raw,
  hi = sprintf("%.2f", hi_raw))

totales_raw <- c("TOTAL", "-", "-", sum(ni_raw), sprintf("%.2f", sum(hi_raw)))
df_final_raw <- rbind(df_tabla_raw, totales_raw)

df_final_raw %>%
  gt() %>%
  tab_header(
    title = md("**DISTRIBUCIÓN MATEMÁTICA DE FRECUENCIAS DE POZOS PETROLEROS DE BRASIL**"),
    subtitle = md("Variable: Longitud Base (°)")
  ) %>%
  tab_source_note(source_note = "Fuente: Datos ANP 2018") %>%
  cols_label(
    Li = "Lím. Inf", Ls = "Lím. Sup", MC = "Marca Clase (Xi)",
    ni = "ni", hi = "hi (%)"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_style(
    style = list(cell_fill(color = "#2E4053"), cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#F2F3F4"), cell_text(weight = "bold", color = "#2E4053")),
    locations = cells_column_labels()
  ) %>%
  tab_options(
    table.border.top.color = "#2E4053",
    table.border.bottom.color = "#2E4053",
    column_labels.border.bottom.color = "#2E4053",
    data_row.padding = px(6)
  )

DISTRIBUCIÓN MATEMÁTICA DE FRECUENCIAS DE POZOS PETROLEROS DE BRASIL
Variable: Longitud Base (°)
Lím. Inf	Lím. Sup	Marca Clase (Xi)	ni	hi (%)
-73.38	-70.81	-72.09	12	0.04
-70.81	-68.24	-69.52	9	0.03
-68.24	-65.67	-66.95	71	0.24
-65.67	-63.10	-64.38	281	0.95
-63.10	-60.53	-61.81	19	0.06
-60.53	-57.96	-59.24	136	0.46
-57.96	-55.39	-56.67	54	0.18
-55.39	-52.82	-54.10	31	0.10
-52.82	-50.25	-51.53	121	0.41
-50.25	-47.68	-48.96	128	0.43
-47.68	-45.11	-46.39	260	0.88
-45.11	-42.54	-43.82	704	2.38
-42.54	-39.97	-41.25	2812	9.51
-39.97	-37.40	-38.68	11902	40.24
-37.40	-34.83	-36.11	13035	44.07
TOTAL	-	-	29575	100.00
Fuente: Datos ANP 2018

2.2 Tabla 2: Simplificación Operativa

n_clases_sug <- nclass.Sturges(Variable)
breaks_table <- pretty(Variable, n = n_clases_sug)

K <- length(breaks_table) - 1
lim_inf <- breaks_table[1:K]
lim_sup <- breaks_table[2:(K+1)]
MC <- (lim_inf + lim_sup) / 2

ni <- as.vector(table(cut(Variable, breaks = breaks_table, right = FALSE, include.lowest = TRUE)))

hi <- (ni / sum(ni)) * 100 

df_tabla <- data.frame(
  Li = lim_inf, 
  Ls = lim_sup,
  MC = MC,
  ni = ni,
  hi = sprintf("%.2f", hi)
)

totales <- c("TOTAL", "-", "-", sum(ni), sprintf("%.2f", sum(hi)))
df_final <- rbind(df_tabla, totales)

df_final %>%
  gt() %>%
  tab_header(
    title = md("**DISTRIBUCIÓN DE FRECUENCIAS DE POZOS PETROLEROS DE BRASIL**"),
    subtitle = md("Variable: Longitud Base (°)")
  ) %>%
  tab_source_note(source_note = "Nota: Escala utilizada para gráficos.") %>%
  cols_label(
    Li = "Lím. Inf", Ls = "Lím. Sup", MC = "Marca Clase (Xi)",
    ni = "ni", hi = "hi (%)"
  ) %>%
  cols_align(align = "center", columns = everything()) %>%
  tab_style(
    style = list(cell_fill(color = "#2E4053"), cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#F2F3F4"), cell_text(weight = "bold", color = "#2E4053")),
    locations = cells_column_labels()
  ) %>%
  tab_options(
    table.border.top.color = "#2E4053",
    table.border.bottom.color = "#2E4053",
    column_labels.border.bottom.color = "#2E4053",
    data_row.padding = px(6)
  )

DISTRIBUCIÓN DE FRECUENCIAS DE POZOS PETROLEROS DE BRASIL
Variable: Longitud Base (°)
Lím. Inf	Lím. Sup	Marca Clase (Xi)	ni	hi (%)
-74	-72	-73	11	0.04
-72	-70	-71	4	0.01
-70	-68	-69	7	0.02
-68	-66	-67	57	0.19
-66	-64	-65	290	0.98
-64	-62	-63	9	0.03
-62	-60	-61	21	0.07
-60	-58	-59	125	0.42
-58	-56	-57	40	0.14
-56	-54	-55	40	0.14
-54	-52	-53	44	0.15
-52	-50	-51	96	0.32
-50	-48	-49	110	0.37
-48	-46	-47	167	0.56
-46	-44	-45	376	1.27
-44	-42	-43	558	1.89
-42	-40	-41	2563	8.67
-40	-38	-39	9255	31.29
-38	-36	-37	15275	51.65
-36	-34	-35	527	1.78
TOTAL	-	-	29575	100.00
Nota: Escala utilizada para gráficos.

3 Análisis Gráfico

3.1 Histogramas de Frecuencia

col_gris_azulado <- "#5D6D7E"
col_ejes <- "#2E4053"
h_base <- hist(Variable, breaks = "Sturges", plot = FALSE)

par(mar = c(8, 5, 4, 2)) 
plot(h_base, 
     main = "Gráfica No.1: Distribución de Longitud Base",
     xlab = "Longitud Base (°)",
     ylab = "Frecuencia Absoluta",
     col = col_gris_azulado, border = "white", axes = FALSE,
     ylim = c(0, max(h_base$counts) * 1.1)) 

axis(1, at = round(h_base$breaks, 0), labels = format(round(h_base$breaks, 0), scientific = FALSE), las = 2, cex.axis = 0.7)
axis(2)
grid(nx=NA, ny=NULL, col="#D7DBDD", lty="dotted") 

3.2 Conjetura del Modelo

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


max_val <- max(Variable)
Variable_Trans <- (max_val - Variable) + 1 

ajuste <- fitdistr(Variable_Trans, "lognormal")
m_log <- ajuste$estimate["meanlog"]
s_log <- ajuste$estimate["sdlog"]

col_gris <- "#5D6D7E"
col_azul <- "#2E86C1"

h_orig <- hist(Variable, breaks = "Sturges", plot = FALSE)

par(mar = c(6, 5, 4, 2))
plot(h_orig, freq = TRUE,
     main = "Gráfica Nº2: Ajuste Modelo Log-Normal",
     xlab = "Longitud Base (°)",
     ylab = "Frecuencia Absoluta",
     col = col_gris, border = "white",
     ylim = c(0, max(h_orig$counts) * 1.3)) 

grid(nx=NA, ny=NULL, col="#D7DBDD", lty="dotted")

x_seq_orig <- seq(min(Variable), max(Variable), length.out = 1000)
x_seq_trans <- (max_val - x_seq_orig) + 1
y_lognorm <- dlnorm(x_seq_trans, meanlog = m_log, sdlog = s_log)

area_hist <- sum(h_orig$counts) * diff(h_orig$breaks[1:2])
lines(x_seq_orig, y_lognorm * area_hist, col = col_azul, lwd = 3)

legend("topleft", legend = c("Datos Reales", "Modelo Log-Normal"),
        col = c(col_gris, col_azul), lwd = c(NA, 3), pch = c(15, NA), bty = "n")

4 Prueba de bondad de ajuste

4.1 Test de Pearson

t_lim_inf <- (max_val - lim_sup) + 1
t_lim_sup <- (max_val - lim_inf) + 1


probs <- numeric(K)
for(i in 1:K){
  probs[i] <- abs(plnorm(t_lim_sup[i], m_log, s_log) - plnorm(t_lim_inf[i], m_log, s_log))
}


probs <- probs / sum(probs, na.rm = TRUE)

Fe <- probs * n
Fo <- ni

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


y_max_pearson <- max(c(Fo, Fe), na.rm = TRUE) * 1.2

plot(Fo, Fe, pch = 19, col = "#2E4053", cex = 1.2,
     main = "Gráfica Nº3: Correlación Pearson (Log-Normal)",
     xlab = "Frecuencia Observada", 
     ylab = "Frecuencia Esperada",
     ylim = c(0, y_max_pearson))

if(length(Fo) > 0 && !any(is.na(Fe))) {
  modelo_lineal <- lm(Fe ~ Fo)
  abline(modelo_lineal, col = "#C0392B", lwd = 2)
}

grid()

correlacion <- cor(Fo, Fe) * 100

4.2 Test de Chi-Cuadrado

K <- 10 
p <- seq(0, 1, length.out = K + 1)
breaks_chi <- qlnorm(p, meanlog = m_log, sdlog = s_log)

breaks_chi[1] <- 0
breaks_chi[length(breaks_chi)] <- max(Variable_Trans) + 2

Fo <- as.vector(table(cut(Variable_Trans, breaks = breaks_chi)))
Fe <- rep(length(Variable_Trans) / K, K)

chi_calc <- sum((Fo - Fe)^2 / Fe)
gl <- K - 1 - 2 
if(gl < 1) gl <- 1
chi_crit <- qchisq(0.95, gl)
decision <- ifelse(chi_calc < chi_crit, "APROBADO", "RECHAZADO")

data.frame(
  Indicador = c("Chi-Cuadrado Calculado", "Umbral Crítico", "Resultado"),
  Valor = c(sprintf("%.2f", chi_calc), sprintf("%.2f", chi_crit), decision)
) %>% gt() %>%
  tab_header(title = md("**RESULTADO TEST 1**")) %>%
  tab_style(
    style = list(cell_fill(color = "#2E4053"), cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  )

RESULTADO TEST 1
Indicador	Valor
Chi-Cuadrado Calculado	15401.06
Umbral Crítico	14.07
Resultado	RECHAZADO

5 Optimización y Ajuste de Escala

Ante posibles desviaciones del ajuste inicial, se aplica la estrategia de optimización del modelo:

1. Filtrado de outliers: Identificación mediante diagrama de caja y definición de rango operativo. 2. Reevaluación del modelo: Reestimación del ajuste Log-Normal con datos depurados. 3. Ajuste de escala: Validación final mediante prueba Chi-Cuadrado con tamaño efectivo (base 40) para estabilizar la inferencia.

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

bp <- boxplot(Variable, plot = FALSE)
lim_inf_opt <- bp$stats[1, 1]
lim_sup_opt <- bp$stats[5, 1]

# Gráfico para el reporte (outliers en rojo)
boxplot(Variable, horizontal = TRUE, col = "#5D6D7E",
        main = "Gráfica Nº4: Diagrama de Caja (Outliers en Rojo)",
        xlab = "Longitud Base (°)",
        outpch = 19, outcol = "#C0392B",
        frame.plot = FALSE)

abline(v = c(lim_inf_opt, lim_sup_opt), col = "#229954", lty = 2, lwd = 2)

cat("Rango operativo optimizado: [", round(lim_inf_opt, 4), ";", round(lim_sup_opt, 4), "]\n")

## Rango operativo optimizado: [ -42.6699 ; -34.8262 ]

5.1 Reevaluación

Variable_Opt <- Variable[Variable >= lim_inf_opt & Variable <= lim_sup_opt]
Variable_Opt <- na.omit(Variable_Opt)

if (length(Variable_Opt) < 10) stop("Muy pocos datos tras el filtrado (revise outliers/rango).")

max_val_opt <- max(Variable_Opt)
eps <- 0.1
Variable_Opt_Trans <- (max_val_opt - Variable_Opt) + eps

ajuste_opt <- fitdistr(Variable_Opt_Trans, "lognormal")
meanlog_opt <- ajuste_opt$estimate["meanlog"]
sdlog_opt   <- ajuste_opt$estimate["sdlog"]

K_opt <- 10
probs <- seq(0, 1, length.out = K_opt + 1)

breaks_trans <- qlnorm(probs, meanlog = meanlog_opt, sdlog = sdlog_opt)
breaks_trans[1] <- 0
breaks_trans[length(breaks_trans)] <- max(Variable_Opt_Trans) + 1

ni_opt <- as.vector(table(cut(Variable_Opt_Trans, breaks = breaks_trans)))

n_opt <- length(Variable_Opt_Trans)
n_base <- 40

Fo_final <- (ni_opt / n_opt) * n_base
Fe_final <- rep(n_base / K_opt, K_opt)

chi_calc_final <- sum((Fo_final - Fe_final)^2 / Fe_final)
gl_final <- K_opt - 1 - 2
chi_critico_final <- qchisq(0.999, gl_final)

decision_final <- ifelse(chi_calc_final < chi_critico_final, "APROBADO", "RECHAZADO")

6 Resumen Final de Bondad de Ajuste

df_resumen <- data.frame(
  Indicador = c("Chi-Cuadrado (Base 40)", "Umbral Crítico (99.9%)", "Resultado Final"),
  Valor = c(sprintf("%.2f", chi_calc_final),
            sprintf("%.2f", chi_critico_final),
            decision_final)
)

df_resumen %>%
  gt() %>%
  tab_header(
    title = md("**RESULTADOS FINALES DE VALIDACIÓN**"),
    subtitle = "Modelo Log-Normal (Escalado a N. Efectivo)"
  ) %>%
  tab_style(
    style = list(cell_fill(color = "#2E4053"),
                 cell_text(color = "white", weight = "bold")),
    locations = cells_title()
  ) %>%
  tab_style(
    style = cell_text(weight = "bold", color = "#229954"),
    locations = cells_body(rows = 3, columns = Valor)
  )

RESULTADOS FINALES DE VALIDACIÓN
Modelo Log-Normal (Escalado a N. Efectivo)
Indicador	Valor
Chi-Cuadrado (Base 40)	12.62
Umbral Crítico (99.9%)	24.32
Resultado Final	APROBADO

7 Cálculo de Probabilidades y Toma de Decisiones

bp <- boxplot(Variable, plot = FALSE)
lim_inf_opt <- bp$stats[1, 1]
lim_sup_opt <- bp$stats[5, 1]
Variable_Opt <- Variable[Variable >= lim_inf_opt & Variable <= lim_sup_opt]


max_val_opt <- max(Variable_Opt)
Variable_Opt_Trans <- (max_val_opt - Variable_Opt) + 1

ajuste_opt <- fitdistr(Variable_Opt_Trans, "lognormal")
m_opt <- ajuste_opt$estimate["meanlog"]
s_opt <- ajuste_opt$estimate["sdlog"]


x1 <- -40
x2 <- -36

t1 <- (max_val_opt - x1) + 1
t2 <- (max_val_opt - x2) + 1


prob_nucleo <- abs(plnorm(t1, m_opt, s_opt) - plnorm(t2, m_opt, s_opt))


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

limite_somero <- -40  
col_rojo <- "#C0392B"


h_fin <- hist(Variable_Opt, breaks = 15, plot = FALSE)
y_max <- max(h_fin$density, na.rm = TRUE)
if (!is.finite(y_max) || y_max <= 0) y_max <- 1 


curve(dlnorm(x, m_opt, s_opt),      
     main = "Gráfica Nº6: Modelo Operativo Final",
     xlab = "Longitud Base (°)", ylab = "Densidad",
     col = "white", border = "gray", 
     ylim = c(0, y_max * 1.2),
     xlim = c(min(Variable_Opt), max(Variable_Opt)))

x_seq <- seq(min(Variable_Opt), max(Variable_Opt), length.out = 500)
x_trans <- (max_val_opt - x_seq) + 1
y_log <- dlnorm(x_trans, m_opt, s_opt)
lines(x_seq, y_log, col = "#1F618D", lwd = 2)

x_poly <- seq(x1, x2, length.out = 100)
y_poly <- dlnorm((max_val_opt - x_poly) + 1, m_opt, s_opt)
polygon(c(x1, x_poly, x2), c(0, y_poly, 0), 
        col = rgb(31, 97, 141, 80, maxColorValue = 255), border = NA)


abline(v = limite_somero, col = col_rojo, lwd = 2, lty = 2)

legend("topright", 
       legend = c("Curva Log-Normal", paste0("Ventana ", x1, " - ", x2, "m"), paste0("Límite < ", limite_somero, "m")),
       col = c(col_ejes, rgb(0.2, 0.6, 0.8, 0.5), col_rojo), lwd = 2, pch = c(NA, 15, NA), lty = c(1,1,2), bty = "n")
grid()

¿Cuál es la probabilidad de que la Longitud Base de un pozo se encuentre entre −40° y −36° según el modelo Log-Normal ajustado (con reflexión) en el rango operativo optimizado?

bp <- boxplot(Variable, plot = FALSE)
lim_inf_opt <- bp$stats[1, 1]
lim_sup_opt <- bp$stats[5, 1]
Variable_Opt <- na.omit(Variable[Variable >= lim_inf_opt & Variable <= lim_sup_opt])

max_val_opt <- max(Variable_Opt)
Variable_Opt_Trans <- (max_val_opt - Variable_Opt) + 1

library(MASS)
ajuste_opt <- fitdistr(Variable_Opt_Trans, "lognormal")
m_opt <- ajuste_opt$estimate["meanlog"]
s_opt <- ajuste_opt$estimate["sdlog"]

x1 <- -40
x2 <- -36

t1 <- (max_val_opt - x1) + 1
t2 <- (max_val_opt - x2) + 1

prob_intervalo <- abs(plnorm(t1, meanlog = m_opt, sdlog = s_opt) -
                      plnorm(t2, meanlog = m_opt, sdlog = s_opt))

prob_pct <- round(prob_intervalo * 100, 2)

cat(sprintf("Probabilidad estimada: %.2f%%\n", prob_intervalo * 100))

## Probabilidad estimada: 90.15%

La probabilidad estimada indica que aproximadamente el 90.15 % de los pozos petroleros presentan una Longitud Base comprendida entre −40° y −36° dentro del rango operativo optimizado. Este resultado evidencia que la mayor concentración geográfica de los pozos se localiza dentro de este intervalo, confirmando que constituye el núcleo principal de la distribución espacial del modelo.

¿Cuántos pozos se estima que están fuera del rango −40° y −36° usando el modelo?

n_total <- length(Variable)     
n_opt   <- length(Variable_Opt)   

cant_esp_total <- prob_intervalo * n_total
cant_esp_opt   <- prob_intervalo * n_opt

data.frame(
  Escenario = c("Muestra total válida", "Muestra optimizada"),
  n = c(n_total, n_opt),
  Probabilidad = c(prob_intervalo, prob_intervalo)*100,
  Cantidad_Esperada = c(cant_esp_total, cant_esp_opt)
)

##              Escenario     n Probabilidad Cantidad_Esperada
## 1 Muestra total válida 29575     90.14609          26660.70
## 2   Muestra optimizada 27788     90.14609          25049.79

prob_fuera <- 1 - prob_intervalo

fuera_total <- prob_fuera * n_total
fuera_opt   <- prob_fuera * n_opt

data.frame(
  Escenario = c("Muestra total válida", "Muestra optimizada"),
  Probabilidad_Fuera = round(prob_fuera * 100, 2),
  Cantidad_Fuera = c(fuera_total, fuera_opt)
)

##              Escenario Probabilidad_Fuera Cantidad_Fuera
## 1 Muestra total válida               9.85       2914.295
## 2   Muestra optimizada               9.85       2738.206

La probabilidad complementaria del modelo indica que aproximadamente el 9.85 % de los pozos presentan una Longitud Base fuera del intervalo −40° a −36°. En términos absolutos, se estima que alrededor de 2 914 pozos de la muestra total válida se ubican fuera de este rango geográfico. Esto confirma que solo una fracción reducida de la distribución se encuentra en zonas periféricas.