0. Carga de librerías
library(gt)
library(dplyr)
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
1. Carga de datos
df <- read.csv(
"waterPollution.csv",
header = TRUE,
sep = ",",
stringsAsFactors = FALSE
)
2. Selección de la variable
# Extraer variable
Venue <- df$VenueCount
# Eliminar valores faltantes
Venue <- Venue[!is.na(Venue)]
3. Tabla de distribución de frecuencias
# ==========================================
# TABLA DE DISTRIBUCIÓN DE FRECUENCIAS
# ==========================================
tabla_freq <- as.data.frame(
table(Venue)
)
tabla_freq <- tabla_freq[
order(as.numeric(as.character(tabla_freq$Venue))),
]
colnames(tabla_freq) <- c(
"Venue",
"ni"
)
N <- sum(tabla_freq$ni)
tabla_freq$hi <- round(
tabla_freq$ni / N * 100,
2
)
tabla_freq$p_s <- round(
tabla_freq$ni / N,
4
)
fila_total <- data.frame(
Venue = "TOTAL",
ni = N,
hi = 100,
p_s = round(sum(tabla_freq$p_s),4)
)
tabla_final <- rbind(
tabla_freq,
fila_total
)
tabla_final %>%
gt() %>%
tab_header(
title = md("**Tabla N°1**"),
subtitle = md(
"**Distribución de frecuencias original de Venue Count**"
)
) %>%
cols_label(
Venue = "Venue Count",
ni = "ni",
hi = "hi (%)",
p_s = "p(s)"
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
rows = Venue == "TOTAL"
)
) %>%
opt_row_striping()
| Tabla N°1 |
| Distribución de frecuencias original de Venue Count |
| Venue Count |
ni |
hi (%) |
p(s) |
| 0 |
19410 |
97.05 |
0.9705 |
| 1 |
129 |
0.64 |
0.0064 |
| 2 |
113 |
0.56 |
0.0056 |
| 3 |
77 |
0.38 |
0.0039 |
| 4 |
114 |
0.57 |
0.0057 |
| 5 |
65 |
0.32 |
0.0032 |
| 6 |
38 |
0.19 |
0.0019 |
| 7 |
17 |
0.08 |
0.0008 |
| 8 |
5 |
0.03 |
0.0003 |
| 9 |
6 |
0.03 |
0.0003 |
| 11 |
1 |
0.00 |
0.0000 |
| 14 |
4 |
0.02 |
0.0002 |
| 15 |
1 |
0.00 |
0.0000 |
| 16 |
3 |
0.01 |
0.0001 |
| 18 |
2 |
0.01 |
0.0001 |
| 19 |
2 |
0.01 |
0.0001 |
| 22 |
1 |
0.00 |
0.0000 |
| 24 |
2 |
0.01 |
0.0001 |
| 31 |
1 |
0.00 |
0.0000 |
| 33 |
1 |
0.00 |
0.0000 |
| 65 |
2 |
0.01 |
0.0001 |
| 84 |
3 |
0.01 |
0.0001 |
| 100 |
3 |
0.01 |
0.0001 |
| TOTAL |
20000 |
100.00 |
0.9995 |
4. Gráfica de distribución de frecuencias
# ==========================================
# GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIAS
# ==========================================
# Número de barras
n <- length(tabla_freq$Venue)
# Paleta de colores
colores <- colorRampPalette(
c("lightgreen","forestgreen")
)(n)
# Límite del eje Y
ylim_local <- c(
0,
max(tabla_freq$hi) * 1.10
)
# Gráfico de barras
barplot(
height = tabla_freq$hi,
names.arg = tabla_freq$Venue,
main = "Gráfica N°1: Distribución porcentual de Venue Count\nEstudio de contaminación del agua en Europa",
xlab = "Venue Count",
ylab = "Porcentaje",
col = colores,
ylim = ylim_local,
las = 1,
cex.names = 0.9,
cex.main = 0.9
)
abline(
h = pretty(c(0, max(tabla_freq$hi))),
col = "gray80",
lty = 2
)
box()
5. Tratamiento de datos
# ==================================================
# TRATAMIENTO DE DATOS
# ==================================================
# La distribución original presenta una elevada
# concentración de observaciones con valor 0.
# Para el ajuste del modelo de Poisson se trabajará
# únicamente con los valores comprendidos entre 1 y 10,
# ya que representan las categorías de interés para
# el análisis.
# Tratamiento de la variable
Venue_tratado <- Venue[
Venue >= 1 &
Venue <= 10
]
# ==========================================
# TABLA DE FRECUENCIAS TRATADA
# ==========================================
tabla_freq_tratada <- as.data.frame(
table(Venue_tratado)
)
tabla_freq_tratada <- tabla_freq_tratada[
order(as.numeric(as.character(tabla_freq_tratada$Venue_tratado))),
]
colnames(tabla_freq_tratada) <- c(
"Venue",
"ni"
)
# Total
N_tratado <- sum(tabla_freq_tratada$ni)
# Frecuencia relativa
tabla_freq_tratada$hi <- round(
tabla_freq_tratada$ni / N_tratado * 100,
2
)
# Probabilidad observada
tabla_freq_tratada$p_s <- round(
tabla_freq_tratada$ni / N_tratado,
4
)
# ------------------------------------------
# Fila TOTAL
# ------------------------------------------
fila_total <- data.frame(
Venue = "TOTAL",
ni = N_tratado,
hi = 100,
p_s = round(sum(tabla_freq_tratada$p_s),4)
)
tabla_final_tratada <- rbind(
tabla_freq_tratada,
fila_total
)
# ------------------------------------------
# Tabla GT
# ------------------------------------------
tabla_final_tratada %>%
gt() %>%
tab_header(
title = md("**Tabla N°2**"),
subtitle = md(
"**Distribución de frecuencias tratada para el ajuste del modelo de Poisson**"
)
) %>%
cols_label(
Venue = "Venue Count",
ni = "ni",
hi = "hi (%)",
p_s = "p(s)"
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
rows = Venue == "TOTAL"
)
) %>%
opt_row_striping()
| Tabla N°2 |
| Distribución de frecuencias tratada para el ajuste del modelo de Poisson |
| Venue Count |
ni |
hi (%) |
p(s) |
| 1 |
129 |
22.87 |
0.2287 |
| 2 |
113 |
20.04 |
0.2004 |
| 3 |
77 |
13.65 |
0.1365 |
| 4 |
114 |
20.21 |
0.2021 |
| 5 |
65 |
11.52 |
0.1152 |
| 6 |
38 |
6.74 |
0.0674 |
| 7 |
17 |
3.01 |
0.0301 |
| 8 |
5 |
0.89 |
0.0089 |
| 9 |
6 |
1.06 |
0.0106 |
| TOTAL |
564 |
100.00 |
0.9999 |
# ==========================================
# GRÁFICA DE LA DISTRIBUCIÓN TRATADA
# ==========================================
barplot(
tabla_freq_tratada$hi,
names.arg = tabla_freq_tratada$Venue,
main = "Gráfica N°2: Distribución porcentual de Venue Count\nTratamiento de datos para el modelo de Poisson",
xlab = "Venue Count",
ylab = "Frecuencia relativa (%)",
col = "lightgreen",
border = "black",
ylim = c(0, max(tabla_freq_tratada$hi) * 1.10),
las = 1
)
grid(
nx = NA,
ny = NULL,
lty = 2,
col = "gray80"
)
box()
5. Conjetura
# Se conjetura que la variable Venue Count presenta un comportamiento
# compatible con una distribución de Poisson, debido a que corresponde
# a una variable discreta de conteo y la mayor concentración de
# observaciones se encuentra en valores bajos, disminuyendo su
# frecuencia conforme aumenta el número de establecimientos. Esta
# conjetura será verificada mediante las pruebas de bondad de ajuste.
6. Parámetro
# ==========================================
# PARÁMETRO λ
# ==========================================
# Valores de Venue Count
x_mapped <- as.numeric(
as.character(tabla_freq_tratada$Venue)
)
# Probabilidades observadas
p_s_data <- tabla_freq_tratada$p_s
# Esperanza matemática (λ)
lambda <- sum(
x_mapped * p_s_data
)
cat(
"Lambda =",
round(lambda,4),
"\n"
)
## Lambda = 3.2051
7. Sobreposición de la realidad con el modelo de Poisson
# ==========================================
# PROBABILIDADES TEÓRICAS
# ==========================================
# Valores tratados
x_mapped <- as.numeric(as.character(tabla_freq_tratada$Venue))
# Probabilidades observadas
p_s_data <- tabla_freq_tratada$p_s
# Parámetro lambda
lambda <- sum(x_mapped * p_s_data)
# Probabilidades teóricas de Poisson
prob_poisson <- dpois(
x = x_mapped,
lambda = lambda
)
# Convertir a porcentaje
hi_poisson <- prob_poisson * 100
# Mostrar probabilidades
comparacion_poisson <- data.frame(
Venue = tabla_freq_tratada$Venue,
Realidad = round(tabla_freq_tratada$hi, 2),
Poisson = round(hi_poisson, 2)
)
comparacion_poisson
## Venue Realidad Poisson
## 1 1 22.87 13.00
## 2 2 20.04 20.83
## 3 3 13.65 22.25
## 4 4 20.21 17.83
## 5 5 11.52 11.43
## 6 6 6.74 6.11
## 7 7 3.01 2.80
## 8 8 0.89 1.12
## 9 9 1.06 0.40
# ==========================================
# COMPARACIÓN
# ==========================================
datos_poisson <- rbind(
Realidad = tabla_freq_tratada$hi,
`Modelo Poisson` = hi_poisson
)
# ==========================================
# GRÁFICA
# ==========================================
par(mar = c(6.5,4.5,4.5,2))
barplot(
datos_poisson,
beside = TRUE,
names.arg = tabla_freq_tratada$Venue,
col = c("steelblue","firebrick"),
main = "Gráfica N°3: Comparación de la realidad con el modelo de Poisson\nVenue Count",
xlab = "Venue Count",
ylab = "Porcentaje",
ylim = c(0, max(datos_poisson) * 1.25),
las = 1,
cex.names = 0.9
)
abline(
h = pretty(c(0, max(datos_poisson))),
col = "gray80",
lty = 2
)
legend(
"topright",
legend = c("Realidad", "Modelo Poisson"),
fill = c("steelblue", "firebrick"),
bty = "n"
)
box()
9.1 Test de Pearson
# ==========================================
# TEST DE PEARSON
# ==========================================
# Probabilidades teóricas
P_teorica <- hi_poisson / 100
# Frecuencias observadas
fo_pearson <- tabla_freq_tratada$ni
# Total de observaciones
N <- sum(fo_pearson)
# Frecuencias esperadas
fe_pearson <- N * P_teorica
# Coeficiente de Pearson
Coef_Pearson <- cor(
fo_pearson,
fe_pearson
) * 100
cat(
"Coeficiente de Pearson (%):",
round(Coef_Pearson,2),
"\n"
)
## Coeficiente de Pearson (%): 84.93
9.2 Test Chi-cuadrado
# ==========================================
# TEST CHI-CUADRADO
# ==========================================
# Total de observaciones
N_chi <- sum(tabla_freq_tratada$ni)
# Frecuencias relativas observadas
fo <- tabla_freq_tratada$ni / N_chi
# Número de categorías
k <- length(fo)
# Probabilidades teóricas
fe <- dpois(
x_mapped,
lambda
)
# Estadístico Chi-cuadrado
Chi_Calculado <- sum(
(fo - fe)^2 / fe
)
# Grados de libertad
gl <- k - 1
# Valor crítico
Chi_Critico <- qchisq(
0.95,
df = gl
)
cat(
"Chi Calculado:",
round(Chi_Calculado,4),
"\n"
)
## Chi Calculado: 0.1241
cat(
"Chi Crítico:",
round(Chi_Critico,4),
"\n"
)
## Chi Crítico: 15.5073
if(Chi_Calculado < Chi_Critico){
cat(
"\nNo se rechaza H0: El modelo de Poisson es adecuado.\n"
)
}else{
cat(
"\nSe rechaza H0: El modelo de Poisson no es adecuado.\n"
)
}
##
## No se rechaza H0: El modelo de Poisson es adecuado.
10. Cálculo de probabilidades
# ==============================================================================
# 10. CÁLCULO DE PROBABILIDADES
# ==============================================================================
# --------------------------------------
# PREGUNTA DE PORCENTAJE
# ¿Cuál es la probabilidad de que Venue Count
# sea menor o igual a 3?
# --------------------------------------
prob_3 <- ppois(
q = 3,
lambda = lambda
)
prob_3_porcentaje <- prob_3 * 100
cat(
"Probabilidad de que Venue Count sea menor o igual a 3:",
round(prob_3_porcentaje, 2),
"%\n"
)
## Probabilidad de que Venue Count sea menor o igual a 3: 60.14 %
# --------------------------------------
# PREGUNTA DE CANTIDAD
# ¿Cuántas observaciones se espera que tengan
# Venue Count menor o igual a 3?
# --------------------------------------
cantidad_esperada <- prob_3 * length(Venue_tratado)
cat(
"Cantidad esperada de observaciones:",
round(cantidad_esperada),
"observaciones\n"
)
## Cantidad esperada de observaciones: 339 observaciones
# --------------------------------------
# PREGUNTA DE INTERVALO
# ¿Cuál es la probabilidad de que Venue Count
# se encuentre entre 2 y 4?
# --------------------------------------
prob_2_4 <- (
ppois(4, lambda = lambda) -
ppois(1, lambda = lambda)
) * 100
cat(
"Probabilidad entre 2 y 4:",
round(prob_2_4, 2),
"%\n"
)
## Probabilidad entre 2 y 4: 60.92 %
# ==============================================================================
# DEMOSTRACIÓN GRÁFICA
# ==============================================================================
# Valores posibles
x <- 0:10
# Probabilidades del modelo de Poisson
y <- dpois(
x,
lambda = lambda
) * 100
# Colores de las barras
colores <- rep("orange", length(x))
# Verde para P(X ≤ 3)
colores[x <= 3] <- "darkgreen"
# Azul para P(2 ≤ X ≤ 4)
colores[x >= 2 & x <= 4] <- "blue"
# Gráfico
barplot(
height = y,
names.arg = x,
col = colores,
border = "black",
ylim = c(0, max(y) * 1.25),
main = "Gráfica N°4: Cálculo de Probabilidades\nModelo de Poisson",
xlab = "Venue Count",
ylab = "Probabilidad (%)"
)
grid(
nx = NA,
ny = NULL,
lty = 2,
col = "gray80"
)
legend(
"topright",
legend = c(
"Modelo de Poisson",
"P(Venue ≤ 3)",
"P(2 ≤ Venue ≤ 4)"
),
fill = c(
"orange",
"darkgreen",
"blue"
),
bty = "o",
bg = "white"
)
box()
11. intervalo confianza
# ==============================================================================
# 11. INTERVALO DE CONFIANZA
# ==============================================================================
# Variable tratada
Venue_num <- Venue_tratado
# Media
media_venue <- mean(Venue_num)
# Desviación estándar
desviacion_venue <- sd(Venue_num)
# Tamaño de muestra
n_venue <- length(Venue_num)
# Nivel de confianza del 95%
error_venue <- 1.96 * (
desviacion_venue /
sqrt(n_venue)
)
# Límites del intervalo
limite_inferior <- round(
media_venue - error_venue,
2
)
limite_superior <- round(
media_venue + error_venue,
2
)
# Crear tabla
tabla_intervalo <- data.frame(
Intervalo = paste0(
"P [",
limite_inferior,
" < μ < ",
limite_superior,
"] = 95%"
)
)
# Mostrar tabla
tabla_intervalo %>%
gt() %>%
tab_header(
title = md("**Tabla N°3**"),
subtitle = md(
"**Intervalo de confianza para Venue Count**"
)
) %>%
opt_row_striping()
| Tabla N°3 |
| Intervalo de confianza para Venue Count |
| Intervalo |
| P [3.05 < μ < 3.36] = 95% |