ESTADÍSTICA INFERENCIAL - VARIABLES CONTINUAS

0.1 Carga de datos

setwd("/cloud/project")
datos <- read.csv("DERRAMES_GLOBALEST.csv", header = TRUE, sep = ";" , dec = ".")
str(datos)

## 'data.frame':    3550 obs. of  23 variables:
##  $ Id                              : int  6786 6250 8220 6241 6216 6620 6262 6229 6201 6221 ...
##  $ Dia                             : int  19 3 21 16 19 7 10 12 18 29 ...
##  $ Mes                             : int  1 6 4 3 12 10 2 5 3 1 ...
##  $ Año                             : chr  "A" "1979" "2010" "1978" ...
##  $ Nombre                          : chr  "Arabian Gulf Spills; Persian Gulf, Kuwait" "IXTOC I; Bahia de Campeche, Mexico" "Deepwater Horizon; Gulf of Mexico" "Amoco Cadiz; Brittany, France" ...
##  $ Ubicacion                       : chr  "Persian Gulf, Kuwait" "Bahia de Campeche, Mexico" "Gulf of Mexico" "Brittany, France" ...
##  $ Latitud                         : chr  "29,5" "19,4083" "28,7367" "48,5833" ...
##  $ Longuitud                       : chr  "48" "-92,325" "-88,3872" "-4,71667" ...
##  $ Amenaza                         : chr  "Oil" "Oil" "Oil" "Oil" ...
##  $ Etiquetas                       : chr  "" "Collision" "" "Grounding" ...
##  $ Tipo_de_crudo                   : chr  "Kuwait crude oil" "IXTOC I crude oil" "Diesel, crude oil" "Arabian light crude, Iranian light crude, Bunker C" ...
##  $ Cantidad_recuperada_superficie  : int  NA NA 1 NA NA NA NA NA NA NA ...
##  $ Cantidad_recuperada_costas      : int  NA NA 1 NA NA NA NA NA NA NA ...
##  $ Cantidad_tratada_biologicamente : int  1 NA 1 1 NA NA NA NA NA NA ...
##  $ Cantidad_dispersada_quimicamente: int  NA 1 1 1 NA NA NA 1 1 1 ...
##  $ Cantidad_quemada                : int  NA 1 1 NA NA NA NA 1 1 1 ...
##  $ Maximo_liberacion_galones       : int  336000009 NA 205000000 68000017 NA NA NA 9240000 36100000 NA ...
##  $ Barreras_de_contencion_flotantes: int  35 12 182 17 3 3 7 8 5 6 ...
##  $ Causa_principal                 : chr  "Daño del tanque  " "Incendio y explosion " "Incendio y explosion " "Daño del tanque " ...
##  $ Volumen_derramados_galones      : chr  "336.000.000" "365.000.000" "600.000.000" "68.000.000" ...
##  $ Respuesta_actual_galones        : chr  "336000000" "252000000" "168000000" "68700000" ...
##  $ Fuente_respuesta                : chr  "description and posts" "posts" "description" "posts" ...
##  $ etiqueta_actualizacion          : chr  "RA updated" "RA newly acquired" "RA updated" "RA updated" ...

1 Latitud

1.1 Extraer variable

datos$Latitud <- as.numeric(as.character(datos$Latitud))

## Warning: NAs introduced by coercion

latitud <- na.omit(datos$Latitud)

1.2 Tabla de Frecuencia

n <- length(latitud)
k <- ceiling(1 + 3.322 * log10(n))   # Sturges

# Límites
min_lat <- min(latitud)
max_lat <- max(latitud)

# Amplitud
A <- (max_lat - min_lat) / k

# Crear breaks
breaks_Latitud <- seq(min_lat, max_lat + A, by = A)

# Crear intervalos
rango_Latitud <- cut(latitud,
                     breaks = breaks_Latitud,
                     right = FALSE,
                     include.lowest = TRUE)

# TABLA DE FRECUENCIAS
TDF_Latitud <- table(rango_Latitud)
tabla_Latitud <- as.data.frame(TDF_Latitud)

# Porcentajes
hi_Latitud <- (tabla_Latitud$Freq / sum(tabla_Latitud$Freq)) * 100
tabla_Latitud$hi <- round(hi_Latitud, 2)

# Frecuencias acumuladas
Niasc_Latitud <- cumsum(tabla_Latitud$Freq)
Hiasc_Latitud <- cumsum(hi_Latitud)
Nidsc_Latitud <- rev(cumsum(rev(tabla_Latitud$Freq)))
Hidsc_Latitud <- rev(cumsum(rev(hi_Latitud)))


# TABLA FINAL
tabla_Latitud_Final <- data.frame(
  Rango_Latitud = levels(rango_Latitud),  
  Frecuencia = tabla_Latitud$Freq,
  Porcentaje = tabla_Latitud$hi,
  Niasc = Niasc_Latitud,
  Hiasc = round(Hiasc_Latitud, 2),
  Nidsc = Nidsc_Latitud,
  Hidsc = round(Hidsc_Latitud, 2)
)

print(tabla_Latitud_Final )

##   Rango_Latitud Frecuencia Porcentaje Niasc  Hiasc Nidsc  Hidsc
## 1   [-78,-54.8)          1       3.57     1   3.57    28 100.00
## 2 [-54.8,-31.7)          2       7.14     3  10.71    27  96.43
## 3  [-31.7,-8.5)          0       0.00     3  10.71    25  89.29
## 4   [-8.5,14.7)          2       7.14     5  17.86    25  89.29
## 5   [14.7,37.8)         13      46.43    18  64.29    23  82.14
## 6     [37.8,61)          9      32.14    27  96.43    10  35.71
## 7     [61,84.2]          1       3.57    28 100.00     1   3.57

1.3 Gráfica

h <- hist(latitud,
          breaks = breaks_Latitud,
          freq = TRUE,
          main = "Gráfica N°1: Histograma de la Latitud de derrames Petroleros a Nivel Mundial",
          cex.main = 0.8, 
          xlab = "Latitud (grados)",
          ylab = "Frecuencia",
          col = terrain.colors(length(breaks_Latitud) - 1),
          border = "black",
          las = 1,
          xlim = c(min(breaks_Latitud), max(breaks_Latitud)))  

1.4 Conjetura de Modelo

1.5 Modelo Normal

# Rango de 0 a 50
latitud_filtrada <- latitud[latitud >= 0 & latitud <= 50]
breaks_intervalo <- seq(0, 50, by = 10)
mu_lat <- mean(latitud_filtrada)
sigma_lat <- sd(latitud_filtrada)

# Valores para la curva normal
x_vals <- seq(min(breaks_intervalo), max(breaks_intervalo), length.out = 100)
y_vals <- dnorm(x_vals, mean = mu_lat, sd = sigma_lat)
y_vals_scaled <- y_vals * length(latitud_filtrada) * (breaks_intervalo[2] - breaks_intervalo[1])

# Histograma con curva normal
hist(latitud_filtrada,
     breaks = breaks_intervalo,
     freq = TRUE,
     main = "Histograma Latitud 30-50 con Ajuste Normal",
     xlab = "Latitud (grados)",
     ylab = "Frecuencia",
     col = terrain.colors(length(breaks_intervalo) - 1),
     border = "black",
     las = 1,
     xlim = c(min(breaks_intervalo), max(breaks_intervalo)))
lines(x_vals, y_vals_scaled, col = "red", lwd = 2)

# Frecuencias observadas y esperadas
n_obs <- length(latitud_filtrada)
u <- mu_lat
sigma <- sigma_lat
limites <- breaks_intervalo
limites[1] <- -Inf
limites[length(limites)] <- Inf
hist_result <- hist(latitud_filtrada, breaks = breaks_intervalo, plot = FALSE)
FO <- hist_result$counts

# Probabilidades según distribución normal
P <- numeric(length(FO))
for (i in 1:length(FO)) {
  P[i] <- pnorm(limites[i+1], mean = u, sd = sigma) - pnorm(limites[i], mean = u, sd = sigma)
}

FE <- P * n_obs
valores_validos <- FE >= 1
FO_filt <- FO[valores_validos]
FE_filt <- FE[valores_validos]
tabla <- data.frame(Intervalo = 1:length(FO_filt), FO = FO_filt, FE = FE_filt)
print(tabla)

##   Intervalo FO       FE
## 1         1  5 3.579548
## 2         2  8 7.588438
## 3         3  6 6.892056
## 4         4  2 3.156104

1.6 Test de Pearson

X2_normal <- sum((FO_filt - FE_filt)^2 / FE_filt)

df_normal <- length(FO_filt) - 1 - 2   # k - 1 - parámetros (μ y σ)

Vc_normal <- qchisq(0.95, df = df_normal)

p_valor_normal <- 1 - pchisq(X2_normal, df = df_normal)

cat("Chi-cuadrado:", X2_normal, "\n")

## Chi-cuadrado: 1.124942

cat("Grados de libertad:", df_normal, "\n")

## Grados de libertad: 1

cat("Valor crítico (95%):", Vc_normal, "\n")

## Valor crítico (95%): 3.841459

cat("p-value:", p_valor_normal, "\n")

## p-value: 0.2888568

cat("¿Se ajusta a la Normal?:", X2_normal < Vc_normal, "\n")

## ¿Se ajusta a la Normal?: TRUE

1.7 Test de Chi-cuadrado

X2_normal <- sum((FO_filt - FE_filt)^2 / FE_filt)
df_normal <- length(FO_filt) - 1 - 2  # Restar 1 y los parámetros estimados (media y sigma)
Vc_normal <- qchisq(0.95, df = df_normal)

cat("Chi-cuadrado:", X2_normal, "\n")

## Chi-cuadrado: 1.124942

cat("Valor crítico (95%):", Vc_normal, "\n")

## Valor crítico (95%): 3.841459

cat("¿Se ajusta a la Normal?:", X2_normal < Vc_normal, "\n")

## ¿Se ajusta a la Normal?: TRUE

1.8 Cálculo de probabilidad

# Parámetros de la distribución normal
u <- mean(latitud_filtrada)   
sigma <- sd(latitud_filtrada)

# Definir el intervalo deseado
lim_inf <- 10   # por ejemplo 10°
lim_sup <- 30   # por ejemplo 30°

# 3Calcular la probabilidad de que una observación esté entre lim_inf y lim_sup
probabilidad_normal <- pnorm(lim_sup, mean = u, sd = sigma) - 
                       pnorm(lim_inf, mean = u, sd = sigma)
probabilidad_normal * 100  # Resultado como porcentaje

## [1] 50.76357

#  Graficar la curva normal y sombrear el área correspondiente
x_vals <- seq(min(latitud), max(latitud), length.out = 1000)
plot(x_vals, dnorm(x_vals, mean = u, sd = sigma),
     type = "l", lwd = 2, col = "skyblue3",
     main = "Gráfica: Probabilidad de Latitud (Normal)",
     xlab = "Latitud (grados)", 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))

# Añadir leyenda
legend("topright", legend = c("Modelo Normal", "Área de Probabilidad"), 
       col = c("skyblue3", "red"), lwd = 2, pch = c(NA, 15))

1.9 Intervalo de Confianza

# Intervalo de confianza de la media al 95%
n <- length(latitud)           # tamaño de muestra
error <- sigma / sqrt(n)       # error estándar

li <- u - 2 * error            # límite inferior
ls <- u + 2 * error            # límite superior

# Mostrar resultados
u      # media muestral

## [1] 28.86364

sigma  # desviación estándar

## [1] 10.4576

error  # error estándar

## [1] 1.976301

li     # límite inferior

## [1] 24.91103

ls     # límite superior

## [1] 32.81624

1.10 Conclusión

La variable Latitud ajusta al modelo normal, con un parámetro estimado de 28.86 y desviación estándar de 10.46. La probabilidad de que derrames hayan ocurrido entre las latitudes de entre 10 y 30 unidades es de: 43.42%. Además, mediante el teorema del límite central, sabemos que la media poblacional se encuentra en el intervalo entre 24.91 y 32.82 (sin depurar y depurada) con un 95% de confianza.

2 Maxima Liberación

2.1 Extraer variable

datos$Maximo_liberacion_galones <- as.numeric(as.character(datos$Maximo_liberacion_galones))
Maximo_liberacion_galones <- na.omit(datos$Maximo_liberacion_galones)

2.2 Tabla de Frecuencia

# Definir rango fijo
rango_min <- 0
rango_max <- 100

# Filtrar datos dentro del rango
x_rango1 <- Maximo_liberacion_galones[
  Maximo_liberacion_galones >= rango_min &
    Maximo_liberacion_galones <= rango_max
]
##Tabla de distribución de frecuencias mediante el teorema de Sturges

n <- length(x_rango1)
k <- round(1 + 3.322 * log10(n))
amplitud <- (rango_max - rango_min) / k

# Intervalos según Sturges en 0-40000
breaks_1 <- seq(rango_min, rango_max, by = amplitud)
breaks_1[length(breaks_1)] <- rango_max

# TABLA DE FRECUENCIAS
cut_1 <- cut(x_rango1, breaks = breaks_1,
             right = TRUE, include.lowest = TRUE)

TDF_1 <- table(cut_1)
Tabla_1 <- as.data.frame(TDF_1)

hi_1 <- (Tabla_1$Freq / sum(Tabla_1$Freq)) * 100
Ni_asc_1 <- cumsum(Tabla_1$Freq)
Hi_asc_1 <- cumsum(hi_1)
Ni_dsc_1 <- rev(cumsum(rev(Tabla_1$Freq)))
Hi_dsc_1 <- rev(cumsum(rev(hi_1)))

# ETIQUETAS LEGIBLES
intervalos_num <- gsub("\\[|\\]|\\(|\\)", "", as.character(Tabla_1$cut_1))
intervalos_partes <- strsplit(intervalos_num, ",")

etiquetas_legibles <- sapply(intervalos_partes, function(x) {
  paste0(
    format(round(as.numeric(x[1]), 0), big.mark = ","),
    " - ",
    format(round(as.numeric(x[2]), 0), big.mark = ",")
  )
})

# TABLA FINAL
Tabla_Final_1 <- data.frame(
  Intervalo = etiquetas_legibles,
  ni = Tabla_1$Freq,
  hi = round(hi_1, 2),
  Ni_asc = Ni_asc_1,
  Hi_asc = round(Hi_asc_1, 2),
  Ni_dsc = Ni_dsc_1,
  Hi_dsc = round(Hi_dsc_1, 2)
)

print(Tabla_Final_1)

##   Intervalo ni    hi Ni_asc Hi_asc Ni_dsc Hi_dsc
## 1    0 - 11 32 16.67     32  16.67    192 100.00
## 2   11 - 22 10  5.21     42  21.88    160  83.33
## 3   22 - 33 16  8.33     58  30.21    150  78.12
## 4   33 - 44 23 11.98     81  42.19    134  69.79
## 5   44 - 56 24 12.50    105  54.69    111  57.81
## 6   56 - 67 13  6.77    118  61.46     87  45.31
## 7   67 - 78  9  4.69    127  66.15     74  38.54
## 8   78 - 89 19  9.90    146  76.04     65  33.85
## 9  89 - 100 46 23.96    192 100.00     46  23.96

2.3 Gráfica

h <- hist(
  x_rango1,
  breaks = breaks_1,
  freq = TRUE,
  main = "Histograma del Máximo de Liberación Petrolera (Sturges 0-40000)",
  xlab = "Máximo de liberación (galones)",
  ylab = "Frecuencia",
  col = terrain.colors(length(breaks_1) - 1),
  border = "black"
)

2.4 Conjetura del modelo

2.5 Modelo Normal

# Seleccionar los datos en el rango
x_11_78 <- x_rango1[x_rango1 >= 11 & x_rango1 <= 78]

# Histograma simple con frecuencia absoluta
hist(
  x_11_78,
  breaks = seq(11, 78, by = 11.16),
  freq = TRUE,
  main = "Histograma del Máximo de Liberación (10 a 80)",
  xlab = "Máximo de liberación (galones)",
  ylab = "Frecuencia",
  col = "lightblue",
  border = "black"
)

# Calcular media y desviación estándar
u <- mean(x_11_78)
sigma <- sd(x_11_78)

u      # media

## [1] 43.68041

sigma  # desviación estándar

## [1] 16.33793

# Histograma con densidad
Hist_x_11_78 <- hist(
  x_11_78,
  breaks = seq(11, 78, by = 11.16),
  freq = FALSE,
  main = "Distribución Normal del Máximo de Liberación (10 a 80)",
  xlab = "Máximo de liberación (galones)",
  ylim = c(0, 0.030),
  ylab = "Densidad",
  col = "lightblue",
  border = "black"
)

# Secuencia de valores para la curva normal
x <- seq(min(x_11_78), max(x_11_78), 0.01)

# Curva normal ajustada
curve(
  dnorm(x, mean = u, sd = sigma),
  col = "red",
  lwd = 2,
  add = TRUE
)

# Verificación del área bajo la curva ≈ 1
Altura <- Hist_x_11_78$density
Area <- sum(Altura * diff(Hist_x_11_78$breaks))  # CORREGIDO
Area  # debe dar ≈ 1

## [1] 1

# Frecuencia simple observada
FO <- Hist_x_11_78$counts
sum(FO)

## [1] 97

# Frecuencia simple esperada
liminf <- Hist_x_11_78$breaks[-length(Hist_x_11_78$breaks)]
limsup <- Hist_x_11_78$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(x_11_78)
FE

## [1]  6.901530 16.403894 24.882387 24.095108 14.895183  5.875876

2.6 Test de Pearson

# Frecuencias observadas
FO <- Hist_x_11_78$counts

# Límites de clase
liminf <- Hist_x_11_78$breaks[-length(Hist_x_11_78$breaks)]
limsup <- Hist_x_11_78$breaks[-1]

# Probabilidades teóricas normales
P <- numeric(length(FO))
for (i in 1:length(FO)) {
  P[i] <- pnorm(limsup[i], u, sigma) - 
    pnorm(liminf[i], u, sigma)
}

# Frecuencias esperadas
FE <- P * length(x_11_78)

# Solo clases válidas (FE > 0)
valores_validos <- FE > 0

# Estadístico de Pearson
X2_normal <- sum((FO[valores_validos] - FE[valores_validos])^2 /
                   FE[valores_validos])

X2_normal

## [1] 5.821382

df_normal <- sum(valores_validos) - 1 - 2
df_normal

## [1] 3

Vc_normal <- qchisq(0.90, df = df_normal)
Vc_normal

## [1] 6.251389

X2_normal < Vc_normal

## [1] TRUE

2.7 Test de Chi-cuadrado

# Estadístico Chi-cuadrado
X2 <- sum(((FO - FE)^2) / FE)

# Solo clases válidas
valores_validos <- FE > 0
X2_normal <- sum((FO[valores_validos] - FE[valores_validos])^2 / FE[valores_validos])
X2_normal

## [1] 5.821382

# Grados de libertad: clases válidas - 1 - 2 parámetros estimados
df_normal <- sum(valores_validos) - 1 - 2

# Valor crítico 90%
Vc_normal <- qchisq(0.90, df = df_normal)
Vc_normal

## [1] 6.251389

# ¿Se ajusta a una distribución normal?
ajuste_normal <- ifelse(X2_normal < Vc_normal,
                        "Se ajusta a la distribución normal",
                        "No se ajusta a la distribución normal")
ajuste_normal

## [1] "Se ajusta a la distribución normal"

2.8 Cálculo de probabilidad

# Definir el intervalo deseado
lim_inf <- 20
lim_sup <- 50

# Calcular la probabilidad de que una observación esté entre 20 y 50
probabilidad_normal <- pnorm(lim_sup, mean = u, sd = sigma) - 
  pnorm(lim_inf, mean = u, sd = sigma)
probabilidad_normal * 100 # Resultado como porcentaje

## [1] 57.69384

#Graficar la curva normal y sombrear el área correspondiente
x_vals <- seq(min(x_11_78), max(x_11_78), 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 = "Máximo de liberación (galones)",
     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))

# Añadir leyenda
legend("topright", 
       legend = c("Modelo Normal", "Área de Probabilidad"), 
       col = c("skyblue3", "red"), 
       lwd = 2, pch = c(NA, 15))

2.9 Intervalo de Confianza

# Tamaño de la muestra
n <- length(x_11_78)
# Error estándar de la media
error <- sigma / sqrt(n)
# Intervalo de confianza al 95%: media ± 2 * error estándar
li <- u - 2 * error
ls <- u + 2 * error
# Mostrar los valores
u      # media muestral

## [1] 43.68041

sigma  # desviación estándar muestral

## [1] 16.33793

error  # margen de error

## [1] 1.658866

li     # límite inferior del intervalo

## [1] 40.36268

ls     # límite superior del intervalo

## [1] 46.99814

2.10 Conclusión

La variable Máximo de liberación (galones) es el que ajusta al modelo normal, con un parámetro estimado de 43.6804, y desviación estándar de 16.3379. La probabilidad de que una observación esté entre 20 y 50 unidades es de: 43.42%. Además, mediante el teorema del límite central, sabemos que la media poblacional se encuentra en el intervalo entre 40.36 y 46.99 con un 95% de confianza.

3 Respuesta Actual

3.1 Extraer variable

datos$Respuesta_actual_galones <- as.numeric(as.character(datos$Respuesta_actual_galones))

## Warning: NAs introduced by coercion

Respuesta_actual_galones <- na.omit(datos$Respuesta_actual_galones)

3.2 Tabla de Frecuencia

## Rango de mayor frecuencia y presencia de pocos valores atípicos

# Definir rango fijo
rango_min <- 0
rango_max <- 40000

# Filtrar datos dentro del rango
x_rango1 <- Respuesta_actual_galones[
  Respuesta_actual_galones >= rango_min &
    Respuesta_actual_galones <= rango_max
]

## Tabla de distribución de frecuencias mediante el teorema de Sturges

n <- length(x_rango1)
k <- round(1 + 3.322 * log10(n))
amplitud <- (rango_max - rango_min) / k

# Intervalos según Sturges en 0-100
breaks_1 <- seq(rango_min, rango_max, by = amplitud)
breaks_1[length(breaks_1)] <- rango_max

# TABLA DE FRECUENCIAS
cut_1 <- cut(x_rango1, breaks = breaks_1,
             right = TRUE, include.lowest = TRUE)

TDF_1 <- table(cut_1)
Tabla_1 <- as.data.frame(TDF_1)

hi_1 <- (Tabla_1$Freq / sum(Tabla_1$Freq)) * 100
Ni_asc_1 <- cumsum(Tabla_1$Freq)
Hi_asc_1 <- cumsum(hi_1)
Ni_dsc_1 <- rev(cumsum(rev(Tabla_1$Freq)))
Hi_dsc_1 <- rev(cumsum(rev(hi_1)))

# ETIQUETAS LEGIBLES
intervalos_num <- gsub("\\[|\\]|\\(|\\)", "", as.character(Tabla_1$cut_1))
intervalos_partes <- strsplit(intervalos_num, ",")

etiquetas_legibles <- sapply(intervalos_partes, function(x) {
  paste0(
    format(round(as.numeric(x[1]), 0), big.mark = ","),
    " - ",
    format(round(as.numeric(x[2]), 0), big.mark = ",")
  )
})

# TABLA FINAL
Tabla_Final_1 <- data.frame(
  Intervalo = etiquetas_legibles,
  ni = Tabla_1$Freq,
  hi = round(hi_1, 2),
  Ni_asc = Ni_asc_1,
  Hi_asc = round(Hi_asc_1, 2),
  Ni_dsc = Ni_dsc_1,
  Hi_dsc = round(Hi_dsc_1, 2)
)


print(Tabla_Final_1 )

##          Intervalo   ni    hi Ni_asc Hi_asc Ni_dsc Hi_dsc
## 1        0 - 3,330 1134 74.90   1134  74.90   1514 100.00
## 2    3,330 - 6,670  138  9.11   1272  84.02    380  25.10
## 3   6,670 - 10,000   81  5.35   1353  89.37    242  15.98
## 4  10,000 - 13,300   37  2.44   1390  91.81    161  10.63
## 5  13,300 - 16,700   22  1.45   1412  93.26    124   8.19
## 6  16,700 - 20,000   23  1.52   1435  94.78    102   6.74
## 7  20,000 - 23,300   11  0.73   1446  95.51     79   5.22
## 8  23,300 - 26,700   21  1.39   1467  96.90     68   4.49
## 9  26,700 - 30,000   20  1.32   1487  98.22     47   3.10
## 10 30,000 - 33,300    5  0.33   1492  98.55     27   1.78
## 11 33,300 - 36,700   12  0.79   1504  99.34     22   1.45
## 12 36,700 - 40,000   10  0.66   1514 100.00     10   0.66

3.3 Gráfica

h <- hist(
  x_rango1,
  breaks = breaks_1,
  freq = TRUE,
  main = "Histograma de la Respuesta Actual (Sturges 0–100)",
  xlab = "Respuesta actual (galones)",
  ylab = "Frecuencia",
  col = terrain.colors(length(breaks_1) - 1),
  border = "black"
)

3.4 Conjetura del modelo

3.5 Modelo Gamma

#Estimación de parámetros del modelo Gamma
media <- mean(x_rango1)
varianza <- var(x_rango1)

shape <- (media^2) / varianza
scale <- varianza / media

shape

## [1] 0.250189

scale

## [1] 14386.11

HistGamma <- hist(x_rango1, freq = FALSE,
                  breaks = breaks_1,
                  main = "Histograma de Respuesta Actual (0–100) con Curva Gamma",
                  xlab = "Respuesta actual (galones)",
                  ylab = "Densidad",
                  ylim = c(0, 0.0003),
                  col = "lightgreen")

x <- seq(min(x_rango1), max(x_rango1), 0.01)
curve(dgamma(x, shape = shape, scale = scale),
      col = "red", lwd = 2, add = TRUE)

#Frecuencias observadas
FO <- HistGamma$counts

#Frecuencia esperada 
liminf <- HistGamma$breaks[-length(HistGamma$breaks)]
limsup <- HistGamma$breaks[-1]
Pgamma <- numeric(length(FO))
for (i in 1:length(FO)) {
  Pgamma[i] <- pgamma(limsup[i], shape = shape, scale = scale) -
    pgamma(liminf[i], shape = shape, scale = scale)
}
FEgamma <- Pgamma * length(x_rango1)

3.6 Test de Pearson

#Frecuencias observadas
FO <- HistGamma$counts
#Frecuencias esperadas (Gamma)
liminf <- HistGamma$breaks[-length(HistGamma$breaks)]
limsup <- HistGamma$breaks[-1]

Pgamma <- numeric(length(FO))

for (i in 1:length(FO)) {
  Pgamma[i] <- pgamma(limsup[i], shape = shape, scale = scale) -
    pgamma(liminf[i], shape = shape, scale = scale)
}

FEgamma <- Pgamma * length(x_rango1)
#Estadístico de Pearson (Chi-cuadrado)
valores_validos_gamma <- FEgamma > 0

X2_gamma <- sum((FO[valores_validos_gamma] - 
                   FEgamma[valores_validos_gamma])^2 / 
                  FEgamma[valores_validos_gamma])

X2_gamma

## [1] 64.99155

#Grados de libertad
df_gamma <- sum(valores_validos_gamma) - 1 - 2
df_gamma

## [1] 9

#Valor crítico
Vc_gamma <- qchisq(0.95, df = df_gamma)
Vc_gamma

## [1] 16.91898

X2_gamma < Vc_gamma

## [1] FALSE

3.7 Test de Chi-cuadrado

valores_validos_gamma <- FEgamma > 0
X2_gamma <- sum((FO[valores_validos_gamma] - 
                   FEgamma[valores_validos_gamma])^2 / 
                  FEgamma[valores_validos_gamma])
X2_gamma

## [1] 64.99155

#Grados de libertad
df_gamma <- sum(valores_validos_gamma) - 1 - 2
df_gamma

## [1] 9

#Valor crítico
Vc_gamma <- qchisq(0.95, df = df_gamma)
Vc_gamma

## [1] 16.91898

# ¿Se ajusta a distribución Gamma?
X2_gamma < Vc_gamma

## [1] FALSE

3.8 Cálculo de probabilidad

## Definir el intervalo deseado
lim_inf <- 1000
lim_sup <- 4000
#Calcular la probabilidad de que una observación esté entre 1000 y 4000
prob_gamma <- pgamma(lim_sup, shape = shape, scale = scale) -
  pgamma(lim_inf, shape = shape, scale = scale)
prob_gamma * 100# Resultado como porcentaje

## [1] 20.1141

#Graficar la curva gamma y sombrear el área correspondiente
x_vals <- seq(min(x_rango1), max(x_rango1), length.out = 1000)

plot(x_vals, dgamma(x_vals, shape = shape, scale = scale),
     type = "l", lwd = 2, col = "darkgreen",
     main = "Cálculo de Probabilidad – Modelo Gamma",
     xlab = "Respuesta actual (galones)",
     ylab = "Densidad")

#Sombrear el área entre los límites definidos
x_somb <- seq(lim_inf, lim_sup, length.out = 500)
y_somb <- dgamma(x_somb, shape = shape, scale = scale)
polygon(c(x_somb, rev(x_somb)),
        c(y_somb, rep(0, length(y_somb))),
        col = rgb(1, 0, 0, 0.4))

# Añadir leyenda
legend("topright", 
       legend = c("Modelo Gamma", "Área de Probabilidad"), 
       col = c("skyblue3", "red"), 
       lwd = 2, pch = c(NA, 15))

3.9 Intervalo de Confianza

# Tamaño de la muestra
n <- length(x_rango1)

# Media y desviación estándar muestral
u <- mean(x_rango1)
sigma <- sd(x_rango1)

# Error estándar de la media
error <- sigma / sqrt(n)

# Intervalo de confianza al 95%
li <- u - 2 * error
ls <- u + 2 * error

# Mostrar resultados
u      # media muestral

## [1] 3599.247

sigma  # desviación estándar muestral

## [1] 7195.775

error  # margen de error

## [1] 184.9331

li     # límite inferior

## [1] 3229.381

ls     # límite superior

## [1] 3969.113

3.10 Conclusión

La variable Respuesta_actual_galones se intentó modelar mediante una distribución Gamma, obteniéndose una media estimada de 3582.727 y una desviación estándar de 7176.363. La probabilidad de que una observación esté entre 10 y 80 galones es de 68.17% bajo este modelo teórico. Sin embargo, de acuerdo con la prueba Chi-cuadrado, el estadístico obtenido es mayor que el valor crítico, por lo que se rechaza la hipótesis de ajuste, concluyéndose que el modelo Gamma no describe adecuadamente el comportamiento real de los datos. Aun así, mediante el teorema del límite central, se estima que la media poblacional se encuentra entre 3215.07 y 3950.38 con un 95% de confianza.