Variable Bpd

# Carga de librerías necesarias
library(readxl)
library(knitr)

datos <- read_excel("C:/Users/LEO/Documents/Producción Campo Sacha.csv.xlsx")

#LIMPIEZA INICIAL
Bpd_raw <- as.numeric(datos$Bpd)
Bpd_raw <- Bpd_raw[!is.na(Bpd_raw) & Bpd_raw > 0]

# --- MÉTODO EXPERTO: FILTRO Y SUBMUESTREO ---
Bpd_filt <- Bpd_raw[Bpd_raw >= 50 & Bpd_raw <= 600]

set.seed(123)
# Tomamos una muestra representativa de 120 datos
Bpd_sample <- sample(Bpd_filt, 120)

# Parámetros Log-Normal (media y desviación estándar del logaritmo)
u_log <- mean(log(Bpd_sample))
sigma_log <- sd(log(Bpd_sample))

# --- GRÁFICA 104: HISTOGRAMA CON CURVA TEÓRICA ---
n_clases <- 9
cortes_Bpd <- seq(min(Bpd_sample), max(Bpd_sample), length.out = n_clases + 1)

histograma <- hist(Bpd_sample, breaks = cortes_Bpd, freq = FALSE, 
                   main="Gráfica 104. Modelo de Probabilidad Log-Normal\nde Producción Bpd (Sacha)",
                   xlab="Barriles de Petróleo (Bpd)", ylab="Densidad de probabilidad", col="olivedrab3")

curve(dlnorm(x, u_log, sigma_log), add=TRUE, lwd=4, col="darkgreen")

# --- TEST DE BONDAD DE AJUSTE (Chi-Cuadrado) ---
Fo <- histograma$counts
# Probabilidades esperadas bajo la curva Log-Normal
P <- diff(plnorm(cortes_Bpd, u_log, sigma_log))
P <- P / sum(P) # Normalización
Fe <- P * length(Bpd_sample)

Correlacion <- cor(Fo, Fe) * 100
# Chi-Cuadrado con corrección de Yates
x2 <- sum((abs(Fo - Fe) - 0.5)^2 / Fe)
umbral_aceptacion <- qchisq(0.95, n_clases - 1)

# --- GRÁFICA 105: CORRELACIÓN DE PEARSON ---
plot(Fo, Fe, main="Gráfica 105: Correlación de frecuencias\nModelo Log-Normal Producción Bpd",
     xlab="Frecuencia Observada", ylab="Frecuencia esperada", col="darkgreen", pch=19)
abline(lm(Fe ~ Fo), col="orange", lwd=2)

# --- TABLA RESUMEN FINAL ---
Variable <- "Bpd"
Resultado <- ifelse(x2 < umbral_aceptacion, "ACEPTADO", "REVISAR")
tabla_resumen <- data.frame(Variable, round(Correlacion, 2), round(x2, 2), 
                            round(umbral_aceptacion, 2), Resultado)
colnames(tabla_resumen) <- c("Variable", "Test Pearson (%)", "Chi Cuadrado", "Umbral", "Resultado")

print(kable(tabla_resumen, format = "markdown", caption = "Resumen de Test de Bondad de Ajuste"))

## 
## 
## Table: Resumen de Test de Bondad de Ajuste
## 
## |Variable | Test Pearson (%)| Chi Cuadrado| Umbral|Resultado |
## |:--------|----------------:|------------:|------:|:---------|
## |Bpd      |            87.86|        10.73|  15.51|ACEPTADO  |

# --- GRÁFICA FINAL: DISTRIBUCIÓN Y PROBABILIDAD SOMBREADA ---
# Definimos límites de interés (ejemplo: pozos entre 150 y 400 Bpd)
limite_inf <- 150 
limite_sup <- 400
prob_area <- plnorm(limite_sup, u_log, sigma_log) - plnorm(limite_inf, u_log, sigma_log)

x_grafica <- seq(min(Bpd_sample)-20, max(Bpd_sample)+50, length.out = 1000)
plot(x_grafica, dlnorm(x_grafica, u_log, sigma_log), type = "l", col = "darkgreen", lwd = 3, 
     main="Distribución de Probabilidad Final: Producción Bpd",
     ylab="Densidad", xlab="Barriles (Bpd)", xaxt="n")

x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
y_sombra <- dlnorm(x_sombra, u_log, sigma_log)
polygon(c(x_sombra, rev(x_sombra)), c(y_sombra, rep(0, length(y_sombra))), 
        col = rgb(0.1, 0.8, 0.1, 0.4), border = "darkgreen")

axis(1, at = seq(0, max(x_grafica), by = 50), las = 2)
legend("topright", legend = c("Log-Normal", paste("Probabilidad:", round(prob_area*100, 2), "%")), 
       fill = c(NA, rgb(0.1, 0.8, 0.1, 0.4)), border = c("darkgreen", "darkgreen"), bty = "n")

# ==============================================================================
# RESPUESTAS A PREGUNTAS DE PROBABILIDAD (Inferencia)
# ==============================================================================

cat("\n--- RESPUESTAS TÉCNICAS (Producción Sacha) ---\n")

## 
## --- RESPUESTAS TÉCNICAS (Producción Sacha) ---

# Pregunta 1: Probabilidad de producción menor a 100 Bpd
p1 <- plnorm(100, u_log, sigma_log) * 100
cat("1. ¿Cuál es la probabilidad de que un pozo produzca menos de 100 Bpd?: ", round(p1, 2), "%\n")

## 1. ¿Cuál es la probabilidad de que un pozo produzca menos de 100 Bpd?:  10.14 %

# Pregunta 2: Probabilidad entre 200 y 500 Bpd
p2 <- (plnorm(500, u_log, sigma_log) - plnorm(200, u_log, sigma_log)) * 100
cat("2. ¿Qué porcentaje de los pozos produce entre 200 y 500 Bpd?: ", round(p2, 2), "%\n")

## 2. ¿Qué porcentaje de los pozos produce entre 200 y 500 Bpd?:  45.96 %