#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 = 70, # Ajusta el número de intervalos según el rango de los datos
col = "cyan4",
xaxt = "n") # Desactiva los ejes X automáticos
# 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 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)
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
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
#Modelo Normal
u <- mean(Liberacion_petroleo_crudo_edicion)
u
## [1] 1706.834
sigma <- sd(Liberacion_petroleo_crudo_edicion)
sigma
## [1] 7996.945
# Histograma con curva Normal
HistogramaNormal <- hist(Liberacion_petroleo_crudo_edicion, freq = FALSE,
main = "Gráfica No.2: Histograma de Liberación de Petróleo Crudo con Curva Normal",
xlab = "Liberación de Petróleo Crudo",
ylab = "Densidad",
breaks = 70, # Ajusta el número de intervalos según el rango de los datos
col = "cyan4",
cex.lab = 1.5,
axes = FALSE) # Desactivamos los ejes automáticos
# 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 sin notación científica
eje_y <- axTicks(2)
axis(2,
at = eje_y,
labels = format(eje_y, scientific = FALSE),
las = 1, # Dirección de las etiquetas en el eje Y
cex.axis = 0.5) # Ajusta el tamaño de las etiquetas en el eje Y
# Cálculo de la media y desviación estándar
u <- mean(Liberacion_petroleo_crudo_edicion, na.rm = TRUE) # Media
sigma <- sd(Liberacion_petroleo_crudo_edicion, na.rm = TRUE) # Desviación estándar
# Ajustar la curva normal
x <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), 0.01)
lines(x, dnorm(x, mean = u, sd = sigma), col = "red", lwd = 2)
# 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], mean = u, sd = sigma) - pnorm(liminf[i], mean = u, sd = 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
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
# Cálculo de probabilidad entre 200 y 400
lim_inf <- 200
lim_sup <- 400
probabilidad_normal <- pnorm(lim_sup, mean = u, sd = sigma) -
pnorm(lim_inf, mean = u, sd = sigma)
probabilidad_normal * 100
## [1] 0.9823916
# Graficar la curva normal y sombrear el área correspondiente
x_vals <- seq(min(Liberacion_petroleo_crudo_edicion), max(Liberacion_petroleo_crudo_edicion), length.out = 1000)
plot(x_vals, dnorm(x_vals, mean = u, sd = sigma),
type = "l", lwd = 2, col = "skyblue3",
main = "Gráfica. Cálculo de Probabilidades (Normal)",
xlab = "Liberacion de Petróleo Crudo", ylab = "Densidad de probabilidad")
# Sombrear el área entre los límites definidos
x_sombreado <- seq(lim_inf, lim_sup, length.out = 1000)
y_sombreado <- dnorm(x_sombreado, mean = u, sd = sigma)
polygon(c(x_sombreado, rev(x_sombreado)),
c(y_sombreado, rep(0, length(y_sombreado))),
col = rgb(1, 0, 0, 0.5)) # Área sombreada con color rojo semi-transparente
n <- length(Liberacion_petroleo_crudo_edicion)
# Error estándar de la media (SE = sigma / sqrt(n))
error <- sigma / sqrt(n)
# Intervalo de confianza al 95%: media ± 2 * error estándar
li <- u - 2 * error # Límite inferior
ls <- u + 2 * error # Límite superior
# Mostrar los valores
cat("Media muestral (u):", u, "\n")
## Media muestral (u): 1706.834
cat("Desviación estándar (sigma):", sigma, "\n")
## Desviación estándar (sigma): 7996.945
cat("Error estándar de la media:", error, "\n")
## Error estándar de la media: 104.1643
cat("Límite inferior del intervalo de confianza:", li, "\n")
## Límite inferior del intervalo de confianza: 1498.505
cat("Límite superior del intervalo de confianza:", ls, "\n")
## Límite superior del intervalo de confianza: 1915.162