Variables discretas aleatorias
1 Cargar librerías y datos
library(knitr)
setwd("D:/Data")
datos <- read.csv("derrames_globales_.csv", header = TRUE, sep = ";", dec =".")2 Variable Año [1987-2004]
años <- as.numeric(datos$ANo)## Warning: NAs introducidos por coerciónaños <- na.omit(años)
# Filtrar para quitar el primer intervalo bajo (Outlier)
años_filtrados <- años[años >= 1987 & años <= 2004]
# Definición de intervalos (Sexenios)
min_val <- 1987
max_val <- 2004
breaks <- seq(min_val, max_val, by = 6)
if (max(breaks) <= max_val) {
breaks <- c(breaks, max_val + 1)
}
labels_intervalos <- paste(breaks[-length(breaks)], breaks[-1] - 1, sep = "-")
años_intervalos <- cut(años_filtrados,
breaks = breaks,
labels = labels_intervalos,
include.lowest = TRUE,
right = FALSE) 2.1 Tabla de frecuencias
TDF_años <- table(años_intervalos)
Tabla_años <- as.data.frame(TDF_años)
colnames(Tabla_años) <- c("años_intervalos", "ni_FA")
# Cálculos estadísticos
hi_FR <- Tabla_años$ni_FA / sum(Tabla_años$ni_FA)
Niasc <- cumsum(Tabla_años$ni_FA)
Hiasc <- cumsum(hi_FR)
Nidsc <- rev(cumsum(rev(Tabla_años$ni_FA)))
Hidsc <- rev(cumsum(rev(hi_FR)))
# Dataframe Final
Tabla_años_final <- data.frame(
Años = as.character(Tabla_años$años_intervalos),
ni_FA = Tabla_años$ni_FA,
hi_FR = round(hi_FR, 4),
Ni_FAAa = Niasc,
Hi_FRAa = round(Hiasc, 4),
Ni_FAAd = Nidsc,
Hi_FRAd = round(Hidsc, 4)
)
print("--- Tabla de Frecuencias (1987-2004) ---")## [1] "--- Tabla de Frecuencias (1987-2004) ---"print(Tabla_años_final)## Años ni_FA hi_FR Ni_FAAa Hi_FRAa Ni_FAAd Hi_FRAd
## 1 1987-1992 353 0.3378 353 0.3378 1045 1.0000
## 2 1993-1998 349 0.3340 702 0.6718 692 0.6622
## 3 1999-2004 343 0.3282 1045 1.0000 343 0.32822.2 Gráfica 1: Distribución aleatoria de años
barplot(Tabla_años_final$ni_FA,
names.arg = Tabla_años_final$Años,
col = "salmon3",
main = "Distribución aleatoria de años",
xlab = "Intervalos de años",
ylab = "Cantidad de Derrames",
ylim = c(0, max(Tabla_años_final$ni_FA) * 1.2),
cex.names= 0.7)
2.3 Conjetura de Modelo Uniforme
FO_absoluta <- Tabla_años_final$ni_FA
k <- nrow(Tabla_años_final)
n_total <- sum(FO_absoluta)
FE_absoluta <- rep(n_total / k, k)
hi_observada <- FO_absoluta / n_total
Prob_uniforme <- rep(1/k, k)2.4 Gráfica 2: Distribución comparativa aleatoria de años (modelo)
barplot(
rbind(hi_observada, Prob_uniforme),
beside = TRUE,
col = c("salmon3", "yellow3"),
legend = c("Real (FO)", "Uniforme (FE)"),
names.arg = Tabla_años_final$Años,
main = "Gráfica N°2: Ajuste del Modelo Uniforme",
xlab = "Intervalos",
ylab = "Probabilidad",
cex.names = 0.7,
args.legend = list(x = "bottomright", bty = "n")
)
2.5 Gráfica 3: Boxplot
boxplot(años_filtrados, horizontal = TRUE, col = "lightblue",
main = "Gráfica N°2.1: Diagrama de Caja - Años",
xlab = "Años")
2.6 Gráfica 4: Correlación FO vs FE
Ni_observada <- cumsum(FO_absoluta)
Ni_esperada <- cumsum(FE_absoluta)
cor_acumulada <- cor(Ni_observada, Ni_esperada)
plot(Ni_observada, Ni_esperada,
main = "Gráfica N°3: FO vs FE",
xlab = "Frecuencia Acumulada Observada",
ylab = "Frecuencia Acumulada Esperada",
pch = 19, col = "darkblue")
abline(lm(Ni_esperada ~ Ni_observada), col = "red", lwd = 2)
text(min(Ni_observada), max(Ni_esperada),
labels = paste("R =", round(cor_acumulada, 4)),
pos = 4, col = "red", cex = 1.2)
2.7 Test de correlación
x2_real <- sum((FO_absoluta - FE_absoluta)^2 / FE_absoluta)
gl <- k - 1
vc <- qchisq(0.95, df = gl)2.8 Tabla de resumen
Variable <- c("Años (Uniforme)")
tabla_resumen <- data.frame(
Variable,
Correlacion_R = round(cor_acumulada, 4),
Chi_Cuadrado = round(x2_real, 2),
Umbral_Aceptacion = round(vc, 2),
Decision = ifelse(x2_real < vc, "Se acepta Modelo", "Se rechaza Modelo")
)
kable(tabla_resumen, format = "markdown", caption = "Resumen Estadístico")| Variable | Correlacion_R | Chi_Cuadrado | Umbral_Aceptacion | Decision |
|---|---|---|---|---|
| Años (Uniforme) | 1 | 0.15 | 5.99 | Se acepta Modelo |
2.9 Segundo Intervalo [2004-2020]->[2004-2011]
años_p1 <- años[años >= 2004 & años <= 2011]
breaks_p1 <- seq(2004, 2012, by = 2)
cut_p1 <- cut(años_p1, breaks = breaks_p1, include.lowest = TRUE, right = FALSE)
Tabla_p1 <- as.data.frame(table(cut_p1))
colnames(Tabla_p1) <- c("Intervalo", "ni_FA")
# Frecuencias Relativas
n_p1 <- sum(Tabla_p1$ni_FA)
hi_p1 <- Tabla_p1$ni_FA / n_p12.9.1 Conjetura de Modelo Binomial
k_p1 <- nrow(Tabla_p1)
x_p1 <- 0:(k_p1 - 1)
n_trials_p1 <- k_p1 - 1
media_x_p1 <- sum(x_p1 * hi_p1)
p_est_p1 <- media_x_p1 / n_trials_p1
# Probabilidades de binomial
prob_p1 <- dbinom(x_p1, size = n_trials_p1, prob = p_est_p1)2.9.2 Tests de correlación
x2_p1 <- sum((hi_p1 - prob_p1)^2 / prob_p1)
pearson_p1 <- cor(hi_p1, prob_p1)
vc_p1 <- qchisq(0.95, df = max(1, k_p1 - 2))2.9.3 Gráfica 5: Distribucion comparativa aleatoria de años [2005-2011]
barplot(rbind(hi_p1, prob_p1), beside = TRUE, col = c("steelblue", "pink"),
names.arg = Tabla_p1$Intervalo, main = "Distribución comparativa [2005-2011]",
legend = c("Real", "Binomial"), args.legend = list(x = "topright", bty = "n"),
ylab = "Probabilidad")
2.9.4 Intervalo de Confianza
media_p1 <- mean(años_p1)
desv_p1 <- sd(años_p1)
n_p1 <- length(años_p1)
# Calculamos el error estándar con la distribución normal (Teorema del Límite Central)
error_p1 <- qnorm(0.975) * (desv_p1 / sqrt(n_p1))
LI_p1 <- media_p1 - error_p1
LS_p1 <- media_p1 + error_p1
# Desviación poblacional necesaria para las conclusiones
desv_pob_p1 <- sqrt(sum((años_p1 - media_p1)^2) / n_p1)
cat("Intervalo de confianza (95%): [", round(LI_p1, 2), ",", round(LS_p1, 2), "]\n")## Intervalo de confianza (95%): [ 2007.3 , 2007.58 ]2.10 Segundo Intervalo [2004-2020]->[2011-2020]
años_p2 <- años[años >= 2012 & años <= 2020]
breaks_p2 <- seq(2012, 2022, by = 2)
cut_p2 <- cut(años_p2, breaks = breaks_p2, include.lowest = TRUE, right = FALSE)
Tabla_p2 <- as.data.frame(table(cut_p2))
colnames(Tabla_p2) <- c("Intervalo", "ni_FA")
# Frecuencias Relativas
n_p2 <- sum(Tabla_p2$ni_FA)
hi_p2 <- Tabla_p2$ni_FA / n_p22.10.1 Conjetura de modelo binomial
k_p2 <- nrow(Tabla_p2)
x_p2 <- 0:(k_p2 - 1)
n_trials_p2 <- k_p2 - 1
media_x_p2 <- sum(x_p2 * hi_p2)
p_est_p2 <- media_x_p2 / n_trials_p2
prob_p2 <- dbinom(x_p2, size = n_trials_p2, prob = p_est_p2)
# Tests de correlación
x2_p2 <- sum((hi_p2 - prob_p2)^2 / prob_p2)
pearson_p2 <- cor(hi_p2, prob_p2)
vc_p2 <- qchisq(0.95, df = max(1, k_p2 - 2))2.10.2 Gráfica 6: Distribucion comparativa aleatoria de años [2012-2020]
barplot(rbind(hi_p2, prob_p2), beside = TRUE, col = c("darkblue", "lightgreen"),
names.arg = Tabla_p2$Intervalo, main = "Distribución comparativa [2012-2020]",
legend = c("Real", "Binomial"), args.legend = list(x = "topright", bty = "n"),
ylab = "Probabilidad")
### Intervalo de confianza
media_p2 <- mean(años_p2)
desv_p2 <- sd(años_p2)
n_p2 <- length(años_p2)
# Calculamos el error estándar con la distribución normal
error_p2 <- qnorm(0.975) * (desv_p2 / sqrt(n_p2))
LI_p2 <- media_p2 - error_p2
LS_p2 <- media_p2 + error_p2
# Desviación poblacional necesaria para las conclusiones
desv_pob_p2 <- sqrt(sum((años_p2 - media_p2)^2) / n_p2)
cat("Intervalo de confianza (95%): [", round(LI_p2, 2), ",", round(LS_p2, 2), "]\n")## Intervalo de confianza (95%): [ 2016.43 , 2016.71 ]2.11 Tabla de resumen
tabla_final <- data.frame(
Periodo = c("2004-2011", "2012-2020"),
Modelo = c("Binomial", "Binomial"),
Pearson_R = c(round(pearson_p1, 4), round(pearson_p2, 4)),
Chi_Cuadrado_Rel = c(round(x2_p1, 4), round(x2_p2, 4)),
Valor_Critico = c(round(vc_p1, 2), round(vc_p2, 2)),
Decision = c(ifelse(x2_p1 < vc_p1, "Se Acepta", "Se Rechaza"),
ifelse(x2_p2 < vc_p2, "Se Acepta", "Se Rechaza"))
)
kable(tabla_final, format = "markdown", caption = "Resumen de Ajuste (Cálculo Relativo)")| Periodo | Modelo | Pearson_R | Chi_Cuadrado_Rel | Valor_Critico | Decision |
|---|---|---|---|---|---|
| 2004-2011 | Binomial | 0.8141 | 0.0766 | 5.99 | Se Acepta |
| 2012-2020 | Binomial | 0.9745 | 0.1877 | 7.81 | Se Acepta |
2.12 Conclusiones
##
##
## === CONCLUSIONES FINALES ===##
## --- Para el Periodo 1 (2004-2011) ---## La variable Años se explica mediante una distribución Binomial, con un intervalo de confianza que se encuentra entre 2007.3 y 2007.58 lo que afirmamos con un 95% de confianza. Además la desviación estándar poblacional es de 1.93 lo que representa la variabilidad en los años observados, indicando que los valores tienden a estar muy cerca de la media calculada.##
##
## --- Para el Periodo 2 (2012-2020) ---## La variable Años se explica mediante una distribución Binomial, con un intervalo de confianza que se encuentra entre 2016.43 y 2016.71 lo que afirmamos con un 95% de confianza. Además la desviación estándar poblacional es de 2.39 lo que representa la variabilidad en los años observados, indicando que los valores tienden a estar muy cerca de la media calculada.3 Variable Barreras de contención
barreras_total <- as.numeric(datos$Barreras_de_contencion_flotantes)
barreras_total <- na.omit(barreras_total)3.1 Primer intervalo [0-30]
barreras_p1 <- barreras_total[barreras_total >= 0 & barreras_total <= 30]
breaks_p1 <- seq(0, 31, by = 5)
labels_p1 <- paste(breaks_p1[-length(breaks_p1)], breaks_p1[-1] - 1, sep = "-")
if(length(labels_p1) < (length(breaks_p1)-1)) {
labels_p1 <- c(labels_p1, paste(tail(breaks_p1, 2)[1], tail(breaks_p1, 1), sep="-"))
}
cut_p1 <- cut(barreras_p1,
breaks = breaks_p1,
labels = labels_p1,
include.lowest = TRUE,
right = FALSE)3.1.1 Tabla de frecuencias
TDF_p1 <- table(cut_p1)
Tabla_p1 <- as.data.frame(TDF_p1)
colnames(Tabla_p1) <- c("Intervalo", "ni_FA")
# Cálculos estadísticos (Acumuladas)
hi_FR_p1 <- Tabla_p1$ni_FA / sum(Tabla_p1$ni_FA)
Niasc_p1 <- cumsum(Tabla_p1$ni_FA)
Hiasc_p1 <- cumsum(hi_FR_p1)
Nidsc_p1 <- rev(cumsum(rev(Tabla_p1$ni_FA)))
Hidsc_p1 <- rev(cumsum(rev(hi_FR_p1)))
# Dataframe Final P1
Tabla_p1_final <- data.frame(
Intervalo = as.character(Tabla_p1$Intervalo),
ni_FA = Tabla_p1$ni_FA,
hi_FR = round(hi_FR_p1, 4),
Ni_FAAa = Niasc_p1,
Hi_FRAa = round(Hiasc_p1, 4),
Ni_FAAd = Nidsc_p1,
Hi_FRAd = round(Hidsc_p1, 4)
)
print("--- Tabla de Frecuencias [0-30] ---")## [1] "--- Tabla de Frecuencias [0-30] ---"print(Tabla_p1_final)## Intervalo ni_FA hi_FR Ni_FAAa Hi_FRAa Ni_FAAd Hi_FRAd
## 1 0-4 2456 0.7150 2456 0.7150 3435 1.0000
## 2 5-9 609 0.1773 3065 0.8923 979 0.2850
## 3 10-14 186 0.0541 3251 0.9464 370 0.1077
## 4 15-19 76 0.0221 3327 0.9686 184 0.0536
## 5 20-24 56 0.0163 3383 0.9849 108 0.0314
## 6 25-29 52 0.0151 3435 1.0000 52 0.01513.1.2 Conjetura de Modelo Geométrico (P1)
n_p1 <- sum(Tabla_p1$ni_FA)
media_p1 <- mean(barreras_p1)
p_hat_p1 <- 1 / (media_p1 + 1) # Estimador de p
puntos_medios_p1 <- breaks_p1[-length(breaks_p1)] + 2.5
prob_geo_p1 <- dgeom(round(puntos_medios_p1), prob = p_hat_p1)
prob_geo_p1 <- prob_geo_p1 / sum(prob_geo_p1) 3.1.3 Tests de correlación P1
x2_p1 <- sum((hi_FR_p1 - prob_geo_p1)^2 / prob_geo_p1)
pearson_p1 <- cor(hi_FR_p1, prob_geo_p1)
vc_p1 <- qchisq(0.95, df = max(1, nrow(Tabla_p1) - 2)) 3.1.4 Gráfica 1: Distribucion comparativa [0-30]
barplot(rbind(hi_FR_p1, prob_geo_p1),
beside = TRUE,
col = c("steelblue", "pink"),
names.arg = Tabla_p1_final$Intervalo,
main = "Distribución aleatoria barrreras de contención",
legend = c("Real", "Geométrico"),
args.legend = list(x = "topright", bty = "n"),
ylab = "Probabilidad")
## Segundo intervalo [31-380]
barreras_p2 <- barreras_total[barreras_total >= 31 & barreras_total <= 380]
breaks_p2 <- seq(31, 381, by = 35) # Ajustado para cubrir hasta 380+
labels_p2 <- paste(breaks_p2[-length(breaks_p2)], breaks_p2[-1] - 1, sep = "-")
cut_p2 <- cut(barreras_p2,
breaks = breaks_p2,
labels = labels_p2,
include.lowest = TRUE,
right = FALSE)3.1.5 Tabla de frecuencias
TDF_p2 <- table(cut_p2)
Tabla_p2 <- as.data.frame(TDF_p2)
colnames(Tabla_p2) <- c("Intervalo", "ni_FA")
# Cálculos estadísticos (Acumuladas)
hi_FR_p2 <- Tabla_p2$ni_FA / sum(Tabla_p2$ni_FA)
Niasc_p2 <- cumsum(Tabla_p2$ni_FA)
Hiasc_p2 <- cumsum(hi_FR_p2)
Nidsc_p2 <- rev(cumsum(rev(Tabla_p2$ni_FA)))
Hidsc_p2 <- rev(cumsum(rev(hi_FR_p2)))
# Dataframe Final P2
Tabla_p2_final <- data.frame(
Intervalo = as.character(Tabla_p2$Intervalo),
ni_FA = Tabla_p2$ni_FA,
hi_FR = round(hi_FR_p2, 4),
Ni_FAAa = Niasc_p2,
Hi_FRAa = round(Hiasc_p2, 4),
Ni_FAAd = Nidsc_p2,
Hi_FRAd = round(Hidsc_p2, 4)
)
print("--- Tabla de Frecuencias [31-380] ---")## [1] "--- Tabla de Frecuencias [31-380] ---"print(Tabla_p2_final)## Intervalo ni_FA hi_FR Ni_FAAa Hi_FRAa Ni_FAAd Hi_FRAd
## 1 31-65 67 0.5826 67 0.5826 115 1.0000
## 2 66-100 20 0.1739 87 0.7565 48 0.4174
## 3 101-135 11 0.0957 98 0.8522 28 0.2435
## 4 136-170 7 0.0609 105 0.9130 17 0.1478
## 5 171-205 2 0.0174 107 0.9304 10 0.0870
## 6 206-240 2 0.0174 109 0.9478 8 0.0696
## 7 241-275 1 0.0087 110 0.9565 6 0.0522
## 8 276-310 2 0.0174 112 0.9739 5 0.0435
## 9 311-345 1 0.0087 113 0.9826 3 0.0261
## 10 346-380 2 0.0174 115 1.0000 2 0.01743.1.6 Conjetura de Modelo Geométrico (P2)
n_p2 <- sum(Tabla_p2$ni_FA)
media_p2 <- mean(barreras_p2)
p_hat_p2 <- 1 / (media_p2 + 1)
puntos_medios_p2 <- breaks_p2[-length(breaks_p2)] + 17.5 # Mitad de 35
prob_geo_p2 <- dgeom(round(puntos_medios_p2), prob = p_hat_p2)
prob_geo_p2 <- prob_geo_p2 / sum(prob_geo_p2) 3.1.7 Tests de correlación P2
x2_p2 <- sum((hi_FR_p2 - prob_geo_p2)^2 / prob_geo_p2)
pearson_p2 <- cor(hi_FR_p2, prob_geo_p2)
vc_p2 <- qchisq(0.95, df = max(1, nrow(Tabla_p2) - 2))3.1.8 Gráfica 2: Distribucion comparativa [31-380]
barplot(rbind(hi_FR_p2, prob_geo_p2),
beside = TRUE,
col = c("darkblue", "lightgreen"),
names.arg = Tabla_p2_final$Intervalo,
main = "Distribución aleatoria barrreras de contención",
legend = c("Real", "Geométrico"),
args.legend = list(x = "topright", bty = "n"),
ylab = "Probabilidad",
cex.names = 0.7)
### Intervalos de Confianza
# P1
error_p1 <- qnorm(0.975) * (sd(barreras_p1) / sqrt(n_p1))
LI_p1 <- media_p1 - error_p1
LS_p1 <- media_p1 + error_p1
desv_pob_p1 <- sqrt(sum((barreras_p1 - mean(barreras_p1))^2) / length(barreras_p1))
# P2
error_p2 <- qnorm(0.975) * (sd(barreras_p2) / sqrt(n_p2))
LI_p2 <- media_p2 - error_p2
LS_p2 <- media_p2 + error_p2
desv_pob_p2 <- sqrt(sum((barreras_p2 - mean(barreras_p2))^2) / length(barreras_p2))3.2 Tabla de resumen
tabla_final <- data.frame(
Periodo = c("Intervalo 0-30", "Intervalo 31-380"),
Pearson_R = c(round(pearson_p1, 4), round(pearson_p2, 4)),
Chi_Cuadrado_Rel = c(round(x2_p1, 4), round(x2_p2, 4)),
Valor_Critico = c(round(vc_p1, 2), round(vc_p2, 2)),
Decision = c(ifelse(x2_p1 < vc_p1, "Se Acepta", "Se Rechaza"),
ifelse(x2_p2 < vc_p2, "Se Acepta", "Se Rechaza"))
)
kable(tabla_final, format = "markdown", caption = "Resumen de Ajuste (Cálculo Relativo)")| Periodo | Pearson_R | Chi_Cuadrado_Rel | Valor_Critico | Decision |
|---|---|---|---|---|
| Intervalo 0-30 | 0.9996 | 0.1774 | 9.49 | Se Acepta |
| Intervalo 31-380 | 0.9256 | 0.2787 | 15.51 | Se Acepta |
Conclusiones
##
##
## === CONCLUSIONES FINALES ===##
## --- Para el Intervalo [0 - 30] ---## La variable Barreras se explica mediante una distribución Geométrica, con un intervalo de confianza que se encuentra entre 3.42 y 3.79 lo que afirmamos con un 95% de confianza. Además la desviación estándar poblacional es de 5.51 lo que representa la variabilidad en las barreras observadas.##
##
## --- Para el Intervalo [31 - 380] ---## La variable Barreras se explica mediante una distribución Geométrica, con un intervalo de confianza que se encuentra entre 70.32 y 95.75 lo que afirmamos con un 95% de confianza. Además la desviación estándar poblacional es de 69.25 lo que representa la variabilidad en las barreras observadas.