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

# Ajusta la ruta a tu archivo local
datos <- read_excel("C:/Users/LEO/Documents/Producción Campo Sacha.csv.xlsx")

# --- LIMPIEZA INICIAL ---
# El BSW debe convertirse a proporción (0 a 1) para la distribución Beta
Bsw_raw <- as.numeric(datos$Bsw_BE)
Bsw_raw <- Bsw_raw[!is.na(Bsw_raw) & Bsw_raw > 0]
# Escalamos a rango 0-1 para el cálculo matemático
Bsw_proc <- Bsw_raw / 100

# --- MÉTODO EXPERTO: FILTRO Y SUBMUESTREO ---
# Filtramos valores entre 0% y 100%
Bsw_filt <- Bsw_proc[Bsw_proc > 0 & Bsw_proc < 1]

set.seed(123)
# Tomamos una muestra representativa de 110 datos
Bsw_sample <- sample(Bsw_filt, 110)

# Parámetros Distribución Beta (alpha y beta) usando método de momentos
media_b <- mean(Bsw_sample)
var_b <- var(Bsw_sample)
alpha_b <- media_b * (((media_b * (1 - media_b)) / var_b) - 1)
beta_b <- (1 - media_b) * (((media_b * (1 - media_b)) / var_b) - 1)

# --- GRÁFICA 104: HISTOGRAMA CON BARRAS CADA 5% (0 A 100) ---
# Multiplicamos la muestra por 100 solo para la gráfica del histograma
muestra_grafica <- Bsw_sample * 100
cortes_Bsw <- seq(0, 100, by = 5)

histograma <- hist(muestra_grafica, breaks = cortes_Bsw, freq = FALSE, 
                   main="Gráfica 104. Modelo de Probabilidad Beta\nde Corte de Agua BSW (Sacha)",
                   xlab="Porcentaje de Agua (%)", ylab="Densidad de probabilidad", col="steelblue1",
                   xaxt = "n")
axis(1, at = cortes_Bsw)

# La curva se dibuja sobre la escala 0-100
curve(dbeta(x/100, alpha_b, beta_b)/100, add=TRUE, lwd=4, col="blue4")

# --- TEST DE BONDAD DE AJUSTE (Chi-Cuadrado) ---
Fo <- histograma$counts
# Probabilidades esperadas bajo la curva Beta
P <- diff(pbeta(cortes_Bsw/100, alpha_b, beta_b))
P <- P / sum(P) # Normalización
Fe <- P * length(Bsw_sample)

# Estabilización de frecuencias esperadas
Fe[Fe < 0.1] <- 0.1

Correlacion <- cor(Fo, Fe) * 100
x2 <- sum((abs(Fo - Fe) - 0.5)^2 / Fe)
n_clases <- length(Fo)
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 Beta Bsw_BE",
     xlab="Frecuencia Observada", ylab="Frecuencia esperada", col="steelblue4", pch=19)
abline(lm(Fe ~ Fo), col="red", lwd=2)

# --- TABLA RESUMEN FINAL ---
Variable <- "Bsw_BE"
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 |
## |:--------|----------------:|------------:|------:|:---------|
## |Bsw_BE   |            92.56|        12.71|  30.14|ACEPTADO  |
# --- GRÁFICA FINAL: DISTRIBUCIÓN Y PROBABILIDAD SOMBREADA ---
# Definimos límites de interés (ejemplo: pozos con BSW crítico entre 80% y 100%)
limite_inf <- 80 
limite_sup <- 100
prob_area <- pbeta(limite_sup/100, alpha_b, beta_b) - pbeta(limite_inf/100, alpha_b, beta_b)

x_grafica <- seq(0, 100, length.out = 1000)
y_grafica <- dbeta(x_grafica/100, alpha_b, beta_b)/100
plot(x_grafica, y_grafica, type = "l", col = "blue4", lwd = 3, 
     main="Distribución de Probabilidad Final: Bsw_BE",
     ylab="Densidad", xlab="Corte de Agua (%)", xaxt="n")

x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
y_sombra <- dbeta(x_sombra/100, alpha_b, beta_b)/100
polygon(c(x_sombra, rev(x_sombra)), c(y_sombra, rep(0, length(y_sombra))), 
        col = rgb(0.2, 0.4, 0.8, 0.4), border = "blue4")

axis(1, at = seq(0, 100, by = 5), las = 2)
legend("topleft", legend = c("Distribución Beta", paste("Prob. BSW > 80%:", round(prob_area*100, 2), "%")), 
       fill = c(NA, rgb(0.2, 0.4, 0.8, 0.4)), border = c("blue4", "blue4"), bty = "n")

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

cat("\n--- RESPUESTAS TÉCNICAS (Corte de Agua Sacha) ---\n")
## 
## --- RESPUESTAS TÉCNICAS (Corte de Agua Sacha) ---
# Pregunta 1: Probabilidad de que el BSW sea mayor al 90%
p1 <- (1 - pbeta(0.90, alpha_b, beta_b)) * 100
cat("1. ¿Cuál es la probabilidad de que un pozo tenga más del 90% de BSW?: ", round(p1, 2), "%\n")
## 1. ¿Cuál es la probabilidad de que un pozo tenga más del 90% de BSW?:  7.49 %
# Pregunta 2: Probabilidad de BSW bajo (menor al 50%)
p2 <- pbeta(0.50, alpha_b, beta_b) * 100
cat("2. ¿Qué porcentaje de los pozos tiene un BSW menor al 50%?: ", round(p2, 2), "%\n")
## 2. ¿Qué porcentaje de los pozos tiene un BSW menor al 50%?:  67.07 %