**LIBERACIÓN DE CRUDO DE PETRÓLEO*+ ## Cargar los datos

setwd("/cloud/project")
datos <- read.csv("DATOS.csv", header = TRUE, sep = ";" , dec = ".")
str(datos)
## 'data.frame':    10190 obs. of  17 variables:
##  $ Distrito_edit                        : chr  "1" "1" "1" "1" ...
##  $ Year_edit_Fecha_del_derrame          : int  2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
##  $ Mes_edit_Fecha_del_derrame           : int  6 3 4 4 6 6 3 9 10 6 ...
##  $ Categoria_Instalaciones              : chr  "Instalacion fija" "Pozos" "Pozos" "Pozos" ...
##  $ Operacion_general                    : chr  "Produccion" "Otro" "Produccion" "Produccion" ...
##  $ Categoria_Fuente                     : chr  NA "Tanques/Almacenamiento" "Lineas/Tuberias" "Infraestructura Fija" ...
##  $ Grupo_causas_probable                : chr  NA "Afectaciones externas" "Factores humanos" "Problemas tecnicos" ...
##  $ Liberacion_petroleo_crudo_edicion    : num  0 0 0 0 0 ...
##  $ Edicion_recuperacion_petroleo_crudo  : num  NA 0 0 0 0 0 0 0 0 NA ...
##  $ Volumen_liberado_Cond_Final          : num  0 0 0 10 0 0 0 1 0 0 ...
##  $ Liberacion_agua_de_produccion_edicion: num  6720 3780 5040 420 10920 ...
##  $ Liberacion_volumen_gas               : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Volumen_condensado_recuperado        : num  NA 0 0 1 0 0 0 0 0 NA ...
##  $ Edicion_Recuperacion_agua_producida  : num  NA 420 4620 0 10920 ...
##  $ Derrame_sobre_agua_limpio            : chr  "NO" "NO" "NO" "NO" ...
##  $ Estado_general                       : chr  "Observaciones tecnicas" NA NA NA ...
##  $ Codigo_area                          : int  1 1 1 1 1 1 1 1 1 3 ...

1 Extraer variable

# Limpieza de datos
Liberacion_petroleo_crudo_edicion <- datos$Liberacion_petroleo_crudo_edicion
Liberacion_petroleo_crudo_edicion <- as.numeric(Liberacion_petroleo_crudo_edicion)
Liberacion_petroleo_crudo_edicion <- Liberacion_petroleo_crudo_edicion[!is.na(Liberacion_petroleo_crudo_edicion) & Liberacion_petroleo_crudo_edicion > 0]

Gráfica de la variable

hist(Liberacion_petroleo_crudo_edicion, 
     main = "Histograma de Liberacion_petroleo_crudo_edicion",
     xlab = "Liberacion_petroleo_crudo_edicion", 
     ylab = "Frecuencia",
     breaks = 70,  # Ajusta el número de intervalos según el rango de los datos
     col = "cyan4", 
     xaxt = "n")  # Desactiva los ejes X automáticos

Ajuste del modelo al histograma

Modelo Log-Normal 

log_Liberacion_petroleo_crudo_edicion <- log(Liberacion_petroleo_crudo_edicion)
ulog <- mean(log_Liberacion_petroleo_crudo_edicion)
ulog
## [1] 6.401712
sigmalog <- sd(log_Liberacion_petroleo_crudo_edicion)
sigmalog
## [1] 1.413285
# Histograma con curva Log-Normal
HistogramaLogNormal <- hist(Liberacion_petroleo_crudo_edicion, freq = FALSE,
                            main = "Gráfica No.1: Histograma de Liberación de Petróleo Crudo con Curva Log-Normal",
                            xlab = "Liberación de Petróleo Crudo",
                            ylab = "Densidad",
                            breaks = 70,
                            col = "cyan4",
                            cex.lab = 1.5,
                            axes = FALSE)


# Personalizar los ejes
# Eje X con saltos cada 10000 sin notación científica
axis(1, at = seq(0, max(Liberacion_petroleo_crudo_edicion), by = 10000),
     labels = format(seq(0, max(Liberacion_petroleo_crudo_edicion), by = 10000), scientific = FALSE))
eje_y <- axTicks(2)
axis(2,
     at = eje_y,
     labels = format(eje_y, scientific = FALSE),
     las = 1,
     cex.axis = 0.5)
# Ajustar la curva log-normal
x <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), length.out = 1000)
lines(x, dlnorm(x, meanlog = ulog, sdlog = sigmalog), col = "red", lwd = 2)

Frecuencia simple observada

FO <- HistogramaLogNormal$counts
FO
##  [1] 5730  120   28    8    3    2    1    1    0    0    0    0    0    0    0
## [16]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## [31]    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0
## [46]    0    0    0    0    0    0    0    0    0    0    1

2 Frecuencia simple esperada

liminf <- HistogramaLogNormal$breaks[-length(HistogramaLogNormal$breaks)]
liminf
##  [1]      0  10000  20000  30000  40000  50000  60000  70000  80000  90000
## [11] 100000 110000 120000 130000 140000 150000 160000 170000 180000 190000
## [21] 200000 210000 220000 230000 240000 250000 260000 270000 280000 290000
## [31] 300000 310000 320000 330000 340000 350000 360000 370000 380000 390000
## [41] 400000 410000 420000 430000 440000 450000 460000 470000 480000 490000
## [51] 500000 510000 520000 530000 540000 550000
limsup <- HistogramaLogNormal$breaks[-1]
limsup
##  [1]  10000  20000  30000  40000  50000  60000  70000  80000  90000 100000
## [11] 110000 120000 130000 140000 150000 160000 170000 180000 190000 200000
## [21] 210000 220000 230000 240000 250000 260000 270000 280000 290000 300000
## [31] 310000 320000 330000 340000 350000 360000 370000 380000 390000 400000
## [41] 410000 420000 430000 440000 450000 460000 470000 480000 490000 500000
## [51] 510000 520000 530000 540000 550000 560000
Plog <- numeric(length(FO))
for (j in 1:length(FO)) {
  Plog[j] <- plnorm(limsup[j], meanlog = ulog, sdlog = sigmalog) -
    plnorm(liminf[j], meanlog = ulog, sdlog = sigmalog)
}
FElog <- Plog * length(Liberacion_petroleo_crudo_edicion)

3 #Test chi-cuadrado

valores_validos_log <- FElog > 0
X2_log <- sum((FO[valores_validos_log] - FElog[valores_validos_log])^2 / FElog[valores_validos_log])
X2_log
## [1] 3870.945
df_log <- sum(valores_validos_log) - 1 - 2
Vc_log <- qchisq(0.90, df = df_log)
Vc_log
## [1] 66.5482
# ¿Se ajusta a distribución log-normal?
X2_log < Vc_log
## [1] FALSE