#UNIVRSIDAD CENTRAL DEL ECUADOR
#Tema:Estadística Inferencial-Variable Continuas
#AUTOR:Grupo 6
#CARGA DE DATOS
setwd("~/")
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 ...
#VARIABLE CONTINUA 1:LIBERACION DE PETRÓLEO CRUDO--------------------------------
# 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]
# Crear histograma con un número específico de intervalos (breaks)
hist(Liberacion_petroleo_crudo_edicion,
main = "Histograma de Liberacion_petroleo_crudo_edicion",
xlab = "Liberacion_petroleo_crudo_edicion",
ylab = "Frecuencia",
breaks = 90,
col = "cyan4",
xaxt = "n")
# Personalizar las etiquetas del eje X
axis(1, at = seq(0, max(Liberacion_petroleo_crudo_edicion), by = 50000),
labels = format(seq(0, max(Liberacion_petroleo_crudo_edicion), by = 50000), scientific = FALSE))
#-------------------
##Modelo Normal
u <- mean(Liberacion_petroleo_crudo_edicion)
u
## [1] 1706.834
sigma <- sd(Liberacion_petroleo_crudo_edicion)
sigma
## [1] 7996.945
# Histograma con frecuencia (freq = TRUE)
HistogramaNormal <- hist(Liberacion_petroleo_crudo_edicion, freq = TRUE,
main = "Gráfica No.2: Histograma de Liberación de Petróleo Crudo con Curva Normal",
cex.main = 0.95,
xlab = "Liberación de Petróleo Crudo",
ylab = "Frecuencia",
breaks = 75,
col = "aquamarine4",
cex.lab = 1.5,
xlim = c(0, 30000),
ylim = c(0, 6000),
xaxt = "n")
# Personalizar el eje X para evitar notación científica
axis(1, at = seq(0, max(Liberacion_petroleo_crudo_edicion), by = 1000),
labels = format(seq(0, max(Liberacion_petroleo_crudo_edicion), by = 1000), scientific = FALSE))
# Generar la curva normal ajustada a las frecuencias
x_vals <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), length.out = 1000)
y_norm <- dnorm(x_vals, mean = u, sd = sigma)
scaling_factor <- max(HistogramaNormal$counts) / max(y_norm)
lines(x_vals, y_norm * scaling_factor, col = "red", lwd = 2)
#Test
# Frecuencia observada
FO <- HistogramaNormal$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
# Frecuencia esperada
liminf <- HistogramaNormal$breaks[-length(HistogramaNormal$breaks)]
limsup <- HistogramaNormal$breaks[-1]
P <- numeric(length(FO))
for (i in 1:length(FO)) {
P[i] <- pnorm(limsup[i], u, sigma) - pnorm(liminf[i], u, sigma)
}
FE <- P * length(Liberacion_petroleo_crudo_edicion)
# Test de chi-cuadrado
valores_validos <- FE > 0
X2_normal <- sum((FO[valores_validos] - FE[valores_validos])^2 / FE[valores_validos])
X2_normal
## [1] 1090894441
df_normal <- sum(valores_validos) - 1 - 2
df_normal
## [1] 4
Vc_normal <- qchisq(0.90, df = df_normal)
Vc_normal
## [1] 7.77944
# ¿Se ajusta a distribución normal?
X2_normal < Vc_normal
## [1] FALSE
##--------------------
#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 = TRUE,
main = "Gráfica No.1: Histograma de Liberación de Petróleo Crudo con Curva Log-Normal",
cex.main=0.8,
xlab = "Liberación de Petróleo Crudo",
ylab = "Frecuencia",
breaks = 75,
col = "pink3",
cex.lab = 1.5,
xlim = c(0, 30000),
ylim = c(0, 6000),
xaxt = "n")
# Personalizar las etiquetas del eje X
axis(1, at = seq(0, max(Liberacion_petroleo_crudo_edicion), by = 5000),
labels = format(seq(0, max(Liberacion_petroleo_crudo_edicion), by = 5000), scientific = FALSE))
# Ajustar la curva log-normal y escalarla para que se ajuste a las frecuencias
x <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), length.out = 1000)
y_lognorm <- dlnorm(x, meanlog = ulog, sdlog = sigmalog)
scaling_factor <- max(HistogramaLogNormal$counts) / max(y_lognorm)
lines(x, y_lognorm * scaling_factor, col = "red", lwd = 2)
#Test
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
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)
#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
df_log
## [1] 53
Vc_log <- qchisq(0.95, df = df_log)
Vc_log
## [1] 70.99345
# ¿Se ajusta a distribución log-normal?
X2_log < Vc_log
## [1] FALSE
#-------------------------------------------------
#Exponencial
lambda <- 1 / mean(Liberacion_petroleo_crudo_edicion)
lambda
## [1] 0.0005858801
# Histograma Exponencial
HistExp <- hist(Liberacion_petroleo_crudo_edicion, freq=TRUE,
main="Gráfica No.10: Histograma de la Liberación de Petróleo Crudo con Curva Exponencial",
cex.main=0.96,
xlab="Liberación de Petróleo Crudo",
ylab="Frecuencia",
breaks = 75,
col = "blue4",
cex.lab = 1.5,
xlim = c(0, 30000),
ylim = c(0, 6000),
xaxt = "n",
yaxt = "n")
# Personalizar el eje X
axis(1, at = seq(0, max(Liberacion_petroleo_crudo_edicion), by = 5000),
labels = format(seq(0, max(Liberacion_petroleo_crudo_edicion), by = 5000), scientific = FALSE))
axis(2,
at = axTicks(2),
labels = format(axTicks(2), scientific = FALSE, digits = 3),
cex.axis = 0.8)
# Ajustar la curva exponencial a la frecuencia
lambda <- 1 / mean(Liberacion_petroleo_crudo_edicion)
x_vals <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), length.out = 1000)
y_exp <- dexp(x_vals, rate = lambda)
scaling_factor <- max(HistExp$counts) / max(y_exp)
lines(x_vals, y_exp * scaling_factor, col = "red", lwd = 2)
#Test:
FO <- HistExp$counts
liminf <- HistExp$breaks[-length(HistExp$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 <- HistExp$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
Pexp <- numeric(length(FO))
for (k in 1:length(FO)) {
Pexp[k] <- pexp(limsup[k], rate = lambda) - pexp(liminf[k], rate = lambda)
}
FEexp <- Pexp * length(Liberacion_petroleo_crudo_edicion)
valores_validos_exp <- FEexp > 0
#Test chi-cuadrado
X2_exp <- sum(((FO[valores_validos_exp] - FEexp[valores_validos_exp])^2) / FEexp[valores_validos_exp])
X2_exp
## [1] 309252578839
df_exp <- sum(valores_validos_exp) - 1 - 1
df_exp
## [1] 5
Vc_exp <- qchisq(0.95, df = df_exp)
Vc_exp
## [1] 11.0705
# ¿Se ajusta a distribución exponencial?
X2_exp < Vc_exp
## [1] FALSE