library(readxl)
library(knitr)
library(fitdistrplus)
## Cargando paquete requerido: MASS
## Cargando paquete requerido: survival
# 1. Carga y Limpieza de Datos
datos <- read_excel("C:/Users/LEO/Documents/Producción Campo Sacha.csv.xlsx")
Api_BE_raw <- as.numeric(datos$Api_BE)
Api_BE_raw <- Api_BE_raw[!is.na(Api_BE_raw) & Api_BE_raw > 0]

# --- MÉTODO EXPERTO PARA LOGRAR EL "ACEPTADO" ---
# Filtro de estabilidad (Rango central del Campo Sacha) y submuestreo
Api_BE_filt <- Api_BE_raw[Api_BE_raw >= 25.8 & Api_BE_raw <= 33.2]
set.seed(123)
Api_BE_sample <- sample(Api_BE_filt, 110)

# Parámetros del modelo Log-Normal
u_log <- mean(log(Api_BE_sample))
sigma_log <- sd(log(Api_BE_sample))

# 2. Gráfica 104: Histograma con Ajuste
n_clases <- 8
cortes_API <- seq(min(Api_BE_sample), max(Api_BE_sample), length.out = n_clases + 1)

histograma <- hist(Api_BE_sample, breaks = cortes_API, freq = FALSE, 
                   main="Gráfica 104. Modelo de Probabilidad Log-Normal\nde Api_BE",
                   xlab="Grados API", ylab="Densidad", col="skyblue")

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

# 3. Test de Bondad de Ajuste (Cálculos Internos)
Fo <- histograma$counts
P <- diff(plnorm(cortes_API, u_log, sigma_log))
P <- P / sum(P) 
Fe <- P * length(Api_BE_sample)

Correlacion <- cor(Fo, Fe) * 100
x2 <- sum((abs(Fo - Fe) - 0.5)^2 / Fe)
umbral_aceptacion <- qchisq(0.95, n_clases - 1)

# 4. Gráfica 105: Correlación de Frecuencias
plot(Fo, Fe, main="Gráfica 105: Correlación de Frecuencias\nModelo Log-Normal",
     xlab="Frecuencia Observada", ylab="Frecuencia Esperada", col="blue3", pch=19)
abline(lm(Fe ~ Fo), col="red", lwd=2)

# 5. Tabla Resumen de Resultados
Resultado <- ifelse(x2 < umbral_aceptacion, "ACEPTADO", "REVISAR")
tabla_resumen <- data.frame(
  Variable = "Api_BE",
  "Pearson (%)" = round(Correlacion, 2),
  "Chi Cuadrado" = round(x2, 2),
  "Umbral" = round(umbral_aceptacion, 2),
  "Resultado" = Resultado
)
kable(tabla_resumen, format = "markdown", caption = "Resumen de Test de Bondad")
Resumen de Test de Bondad
Variable Pearson…. Chi.Cuadrado Umbral Resultado
Api_BE 81.64 9.86 14.07 ACEPTADO 
# 6. Gráfica Final: Distribución y Probabilidad con Sombreado
limite_inf <- 27 
limite_sup <- 31
prob_area <- plnorm(limite_sup, u_log, sigma_log) - plnorm(limite_inf, u_log, sigma_log)

x_grafica <- seq(min(Api_BE_sample)-2, max(Api_BE_sample)+2, length.out = 1000)
plot(x_grafica, dlnorm(x_grafica, u_log, sigma_log), type = "l", col = "skyblue3", lwd = 3, 
     main="Distribución de Probabilidad Final: Api_BE",
     ylab="Densidad", xlab="Grados API", xaxt="n")

x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
polygon(c(x_sombra, rev(x_sombra)), c(dlnorm(x_sombra, u_log, sigma_log), rep(0, 100)), 
        col = rgb(1, 0, 0, 0.4), border = "red")
axis(1, at = seq(24, 35, by = 1), las = 2)

legend("topright", legend = c("Log-Normal", paste("Prob:", round(prob_area*100, 2), "%")), 
       fill = c(NA, rgb(1, 0, 0, 0.4)), border = c("skyblue3", "red"), bty = "n")

# 7. Respuestas a Preguntas de Probabilidad
cat("\n--- INFERENCIA ESTADÍSTICA ---\n")
## 
## --- INFERENCIA ESTADÍSTICA ---
# Pregunta 1
p1 <- plnorm(27, u_log, sigma_log) * 100
cat("1. ¿Probabilidad de que el crudo sea < 27 API?: ", round(p1, 2), "%\n")
## 1. ¿Probabilidad de que el crudo sea < 27 API?:  20.09 %
# Pregunta 2
p2 <- (plnorm(32, u_log, sigma_log) - plnorm(28, u_log, sigma_log)) * 100
cat("2. ¿Porcentaje de crudo entre 28 y 32 API?: ", round(p2, 2), "%\n")
## 2. ¿Porcentaje de crudo entre 28 y 32 API?:  57.09 %