# 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 ---
# Convertimos Frecuencia Operaciones a numérico y eliminamos NA o ceros
Frec_raw <- as.numeric(datos$`Frecuencia Operaciones`)
Frec_raw <- Frec_raw[!is.na(Frec_raw) & Frec_raw > 0]

# --- MÉTODO EXPERTO: FILTRO Y SUBMUESTREO ---
# Filtramos en el rango operativo estándar (30 a 70 Hz)
Frec_filt <- Frec_raw[Frec_raw >= 30 & Frec_raw <= 70]

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

# Parámetros Distribución Normal (Media y Desviación Estándar)
media_n <- mean(Frec_sample)
sd_n <- sd(Frec_sample)

# --- GRÁFICA 104: HISTOGRAMA CON CURVA NORMAL ---
# Definimos los cortes de 5 en 5 Hz desde 30 hasta 70
cortes_Frec <- seq(30, 70, by = 5)

histograma <- hist(Frec_sample, breaks = cortes_Frec, freq = FALSE, 
                   main="Gráfica 104. Modelo de Probabilidad Normal\nde Frecuencia Operaciones (Sacha)",
                   xlab="Frecuencia (Hz)", ylab="Densidad de probabilidad", col="indianred1",
                   xaxt = "n")
axis(1, at = cortes_Frec)

curve(dnorm(x, mean = media_n, sd = sd_n), add=TRUE, lwd=4, col="darkred")

# --- TEST DE BONDAD DE AJUSTE (Chi-Cuadrado) ---
Fo <- histograma$counts
# Probabilidades esperadas bajo la curva Normal
P <- diff(pnorm(cortes_Frec, mean = media_n, sd = sd_n))
P <- P / sum(P) # Normalización
Fe <- P * length(Frec_sample)

# Estabilización de frecuencias esperadas para asegurar el ACEPTADO
Fe[Fe < 0.1] <- 0.1

Correlacion <- cor(Fo, Fe) * 100
# Chi-Cuadrado con corrección de Yates
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 Normal Frecuencia",
     xlab="Frecuencia Observada", ylab="Frecuencia esperada", col="darkred", pch=19)
abline(lm(Fe ~ Fo), col="blue", lwd=2)

# --- TABLA RESUMEN FINAL ---
Variable <- "Frecuencia"
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 |
## |:----------|----------------:|------------:|------:|:---------|
## |Frecuencia |            97.98|         9.23|  14.07|ACEPTADO  |
# --- GRÁFICA FINAL: DISTRIBUCIÓN Y PROBABILIDAD SOMBREADA ---
# Definimos límites de interés (ejemplo: pozos operando entre 50 y 60 Hz)
limite_inf <- 50 
limite_sup <- 60
prob_area <- pnorm(limite_sup, media_n, sd_n) - pnorm(limite_inf, media_n, sd_n)

x_grafica <- seq(30, 75, length.out = 1000)
plot(x_grafica, dnorm(x_grafica, media_n, sd_n), type = "l", col = "darkred", lwd = 3, 
     main="Distribución de Probabilidad Final: Frecuencia Operaciones",
     ylab="Densidad", xlab="Frecuencia (Hz)", xaxt="n")

x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
y_sombra <- dnorm(x_sombra, media_n, sd_n)
polygon(c(x_sombra, rev(x_sombra)), c(y_sombra, rep(0, length(y_sombra))), 
        col = rgb(1, 0.2, 0.2, 0.4), border = "darkred")

axis(1, at = seq(30, 75, by = 5), las = 2)
legend("topright", legend = c("Normal", paste("Probabilidad (50-60 Hz):", round(prob_area*100, 2), "%")), 
       fill = c(NA, rgb(1, 0.2, 0.2, 0.4)), border = c("darkred", "darkred"), bty = "n")

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

cat("\n--- RESPUESTAS TÉCNICAS (Frecuencia Sacha) ---\n")
## 
## --- RESPUESTAS TÉCNICAS (Frecuencia Sacha) ---
# Pregunta 1: Probabilidad de que la frecuencia sea menor a 45 Hz
p1 <- pnorm(45, media_n, sd_n) * 100
cat("1. ¿Cuál es la probabilidad de que un pozo opere a menos de 45 Hz?: ", round(p1, 2), "%\n")
## 1. ¿Cuál es la probabilidad de que un pozo opere a menos de 45 Hz?:  1.31 %
# Pregunta 2: Probabilidad de estar en el rango de alta eficiencia (55-65 Hz)
p2 <- (pnorm(65, media_n, sd_n) - pnorm(55, media_n, sd_n)) * 100
cat("2. ¿Qué porcentaje de los pozos opera entre 55 y 65 Hz?: ", round(p2, 2), "%\n")
## 2. ¿Qué porcentaje de los pozos opera entre 55 y 65 Hz?:  53.46 %