1. Carga de librerías y datos

library(readxl)
library(dplyr)
library(gt)

datos <- read_excel("dataset_deslizamientos.xlsx")

2. Selección y depuración de la variable

Se seleccionó la variable Latitud (°) y se eliminaron los valores faltantes.

latitude <- datos$latitude

latitude <- latitude[!is.na(latitude)]

n_lat <- length(latitude)

n_lat
## [1] 11033

3. Tabla de distribución de frecuencias

3.1 Construcción de intervalos de clase

k_lat <- 12

min_lat <- min(latitude)
max_lat <- max(latitude)

R_lat <- max_lat - min_lat

A_real <- R_lat / k_lat

A_lat <- ifelse(
  A_real <= 2,
  2,
  ifelse(
    A_real <= 5,
    5,
    ifelse(
      A_real <= 10,
      10,
      ceiling(A_real / 10) * 10
    )
  )
)

Li0 <- floor(min_lat / A_lat) * A_lat

Li_lat <- seq(
  Li0,
  by = A_lat,
  length.out = k_lat
)

Ls_lat <- Li_lat + A_lat

MC_lat <- round((Li_lat + Ls_lat) / 2, 2)

3.2 Cálculo de frecuencias

ni_lat <- numeric(k_lat)

for(i in 1:k_lat){

  if(i < k_lat){

    ni_lat[i] <- sum(
      latitude >= Li_lat[i] &
      latitude < Ls_lat[i]
    )

  } else {

    ni_lat[i] <- sum(
      latitude >= Li_lat[i] &
      latitude <= max_lat
    )

  }

}

hi_lat <- (ni_lat / sum(ni_lat)) * 100

Ni_asc <- cumsum(ni_lat)

Ni_dsc <- rev(cumsum(rev(ni_lat)))

Hi_asc <- cumsum(hi_lat)

Hi_dsc <- rev(cumsum(rev(hi_lat)))

3.3 Elaboración de la tabla

TDF_latitude <- data.frame(
  Li = Li_lat,
  Ls = Ls_lat,
  MC = MC_lat,
  ni = ni_lat,
  hi = hi_lat,
  Ni_asc = Ni_asc,
  Ni_dsc = Ni_dsc,
  Hi_asc = Hi_asc,
  Hi_dsc = Hi_dsc
)

TDF_latitude <- rbind(
  TDF_latitude,
  data.frame(
    Li = "TOTAL",
    Ls = "",
    MC = "",
    ni = sum(ni_lat),
    hi = 100,
    Ni_asc = "",
    Ni_dsc = "",
    Hi_asc = "",
    Hi_dsc = ""
  )
)

3.4 Presentación de la tabla

tabla_latitude <- TDF_latitude %>%
  mutate(
    across(
      c(ni, hi, Ni_asc, Ni_dsc, Hi_asc, Hi_dsc, MC),
      as.numeric
    )
  ) %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 1**"),
    subtitle = md(
      "Distribución de frecuencias de la variable Latitud (°) de los eventos a nivel mundial"
    )
  )

tabla_latitude
Tabla N° 1
Distribución de frecuencias de la variable Latitud (°) de los eventos a nivel mundial
Li Ls MC ni hi Ni_asc Ni_dsc Hi_asc Hi_dsc
-50 -40 -45 99 0.8973081 99 11033 0.8973081 100.0000000
-40 -30 -35 147 1.3323665 246 10934 2.2296746 99.1026919
-30 -20 -25 257 2.3293755 503 10787 4.5590501 97.7703254
-20 -10 -15 158 1.4320674 661 10530 5.9911176 95.4409499
-10 0 -5 498 4.5137315 1159 10372 10.5048491 94.0088824
0 10 5 1015 9.1996737 2174 9874 19.7045228 89.4951509
10 20 15 1394 12.6348228 3568 8859 32.3393456 80.2954772
20 30 25 1842 16.6953684 5410 7465 49.0347140 67.6606544
30 40 35 2554 23.1487356 7964 5623 72.1834497 50.9652860
40 50 45 2603 23.5928578 10567 3069 95.7763074 27.8165503
50 60 55 419 3.7976978 10986 466 99.5740053 4.2236926
60 70 65 47 0.4259947 11033 47 100.0000000 0.4259947
TOTAL NA 11033 100.0000000 NA NA NA NA

4. Modelo probabilístico 1: Exponencial

4.1 Conjetura

Se seleccionó el intervalo de latitudes comprendido entre -50° y -20°, debido a que la forma del histograma presenta una disminución progresiva de las frecuencias, sugiriendo un posible comportamiento exponencial.

lat_50_20 <- latitude[latitude >= -50 & latitude <= -20]

Li_50_20 <- seq(-50, -30, by = 10)
LS_50_20 <- Li_50_20 + 10
MC_50_20 <- (Li_50_20 + LS_50_20) / 2

k_50_20 <- length(Li_50_20)

ni_50_20 <- numeric(k_50_20)

for (i in 1:k_50_20) {

  if (i < k_50_20) {

    ni_50_20[i] <- sum(
      lat_50_20 >= Li_50_20[i] &
      lat_50_20 < LS_50_20[i]
    )

  } else {

    ni_50_20[i] <- sum(
      lat_50_20 >= Li_50_20[i] &
      lat_50_20 <= -20
    )

  }

}

hi_50_20 <- (ni_50_20 / sum(ni_50_20)) * 100

NI_50_20 <- cumsum(ni_50_20)

NID_50_20 <- rev(cumsum(rev(ni_50_20)))

HI_50_20 <- cumsum(hi_50_20)

HID_50_20 <- rev(cumsum(rev(hi_50_20)))

Tabla N° 2

Tabla_50_20 <- data.frame(
  Li  = Li_50_20,
  LS  = LS_50_20,
  MC  = MC_50_20,
  ni  = ni_50_20,
  hi  = round(hi_50_20, 2),
  NI  = NI_50_20,
  NID = NID_50_20,
  HI  = round(HI_50_20, 2),
  HID = round(HID_50_20, 2)
)

tabla_gt_50_20 <- Tabla_50_20 %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 2**"),
    subtitle = md(
      "Distribución de frecuencias de la variable Latitud entre (-50° a -20°) de los eventos a nivel mundial"
    )
  ) %>%
  fmt_number(
    columns = c(hi, HI, HID),
    decimals = 2
  ) %>%
  tab_source_note(
    source_note = md(
      "Elaborado por: Grupo 2 – Carrera de Geología"
    )
  )

tabla_gt_50_20
Tabla N° 2
Distribución de frecuencias de la variable Latitud entre (-50° a -20°) de los eventos a nivel mundial
Li LS MC ni hi NI NID HI HID
-50 -40 -45 99 19.68 99 503 19.68 100.00
-40 -30 -35 147 29.22 246 404 48.91 80.32
-30 -20 -25 257 51.09 503 257 100.00 51.09
Elaborado por: Grupo 2 – Carrera de Geología

Gráfica N° 2

HistoLat1 <- hist(
  lat_50_20,
  breaks = seq(-50, -20, by = 10),
  freq = TRUE,
  col = "grey",
  main = "Gráfica 2: Conjetura del modelo exponencial",
  xlab = "Latitud (°)",
  ylab = "Frecuencia absoluta"
)

Fo_abs <- HistoLat1$counts

Fo_abs
## [1]  99 147 257

4.2 Cálculo de parámetros

Para la estimación de los parámetros del modelo exponencial se empleó una linealización mediante logaritmos naturales de las frecuencias observadas.

ni_adj <- ni_50_20

ni_adj[ni_adj == 0] <- 0.5

ajuste_exp <- lm(
  log(ni_adj) ~ MC_50_20
)

summary(ajuste_exp)
## 
## Call:
## lm(formula = log(ni_adj) ~ MC_50_20)
## 
## Residuals:
##        1        2        3 
##  0.02722 -0.05444  0.02722 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 6.714300   0.169454   39.62   0.0161 *
## MC_50_20    0.047698   0.004715   10.12   0.0627 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06668 on 1 degrees of freedom
## Multiple R-squared:  0.9903, Adjusted R-squared:  0.9806 
## F-statistic: 102.3 on 1 and 1 DF,  p-value: 0.06273

Parámetros estimados

a_exp <- coef(ajuste_exp)[1]

b_exp <- coef(ajuste_exp)[2]

lambda_exp <- -b_exp

parametros_exp <- data.frame(
  Parametro = c("a", "b", "Lambda"),
  Valor = c(a_exp, b_exp, lambda_exp)
)

parametros_exp
##   Parametro       Valor
## 1         a  6.71429958
## 2         b  0.04769781
## 3    Lambda -0.04769781

4.3 Modelo exponencial

Una vez estimados los parámetros, se superpone la curva exponencial al histograma de frecuencias observadas para evaluar visualmente el ajuste.

hist(
  lat_50_20,
  breaks = seq(-50, -20, by = 10),
  freq = TRUE,
  col = "grey",
  main = "Gráfica 3: Ajuste del modelo exponencial",
  xlab = "Latitud (°)",
  ylab = "Frecuencia absoluta"
)

x_exp <- seq(
  -50,
  -20,
  length.out = 300
)

y_exp <- exp(
  b_exp * x_exp
)

y_exp <- y_exp /
  max(y_exp) *
  max(ni_50_20)

lines(
  x_exp,
  y_exp,
  col = "red",
  lwd = 2
)

legend(
  "topleft",
  legend = c(
    "Histograma",
    "Modelo exponencial"
  ),
  col = c(
    "grey",
    "red"
  ),
  lwd = c(
    10,
    2
  ),
  bty = "n"
)

5. Validación del modelo exponencial

5.1 Frecuencias esperadas

f_exp <- function(x) {
  exp(b_exp * x)
}

h <- length(Fo_abs)

P <- numeric(h)

for (i in 1:h) {

  P[i] <- integrate(
    f_exp,
    lower = HistoLat1$breaks[i],
    upper = HistoLat1$breaks[i + 1]
  )$value

}

P <- P / sum(P)

Fe_abs <- P * sum(Fo_abs)

Fe_abs
## [1]  96.5978 155.6382 250.7640

5.2 Correlación entre frecuencias observadas y esperadas

Coeficiente de correlación

pearson_r <- cor(Fo_abs, Fe_abs)

pearson_pct <- pearson_r * 100

pearson_pct
## [1] 99.60892

Gráfica de correlación

plot(
  Fo_abs,
  Fe_abs,
  main = "Gráfica 4: Correlación entre frecuencias observadas y esperadas del modelo exponencial",
  xlab = "Frecuencia observada",
  ylab = "Frecuencia esperada",
  pch = 19,
  col = "blue3"
)

abline(
  lm(Fe_abs ~ Fo_abs),
  col = "red",
  lwd = 2
)

Porcentaje de correlación

Correlacion <- cor(Fo_abs, Fe_abs) * 100

Correlacion
## [1] 99.60892

5.3 Prueba de bondad de ajuste Chi-cuadrado

Reagrupación de frecuencias

Fo <- (Fo_abs / sum(Fo_abs)) * 100

Fe <- (Fe_abs / sum(Fe_abs)) * 100

tabla_chi <- data.frame(
  Fo = Fo,
  Fe = Fe
)

Fo_r <- c()
Fe_r <- c()

acum_Fo <- 0
acum_Fe <- 0

for (i in 1:nrow(tabla_chi)) {

  acum_Fo <- acum_Fo + tabla_chi$Fo[i]
  acum_Fe <- acum_Fe + tabla_chi$Fe[i]

  if (acum_Fe >= 5) {

    Fo_r <- c(Fo_r, acum_Fo)
    Fe_r <- c(Fe_r, acum_Fe)

    acum_Fo <- 0
    acum_Fe <- 0

  }

}

if (acum_Fe > 0) {

  Fo_r[length(Fo_r)] <- Fo_r[length(Fo_r)] + acum_Fo

  Fe_r[length(Fe_r)] <- Fe_r[length(Fe_r)] + acum_Fe

}

tabla_chi_final <- data.frame(
  Fo = Fo_r,
  Fe = Fe_r
)

tabla_chi_final
##         Fo       Fe
## 1 19.68191 19.20433
## 2 29.22465 30.94199
## 3 51.09344 49.85368

Estadístico Chi-cuadrado

x2 <- sum(
  (tabla_chi_final$Fo - tabla_chi_final$Fe)^2 /
    tabla_chi_final$Fe
)

x2
## [1] 0.1380221

Grados de libertad

k_validas <- nrow(tabla_chi_final)

m <- 1

gl <- k_validas - 1 - m

gl
## [1] 1

Valor crítico

alpha <- 0.05

chi_crit <- qchisq(
  1 - alpha,
  gl
)

chi_crit
## [1] 3.841459

Decisión

if (x2 < chi_crit) {

  "NO se rechaza H0: el modelo exponencial es adecuado"

} else {

  "SE rechaza H0: el modelo exponencial NO es adecuado"

}
## [1] "NO se rechaza H0: el modelo exponencial es adecuado"
x2 < chi_crit
## [1] TRUE

6. Aplicación del modelo exponencial

6.1 Función de densidad de probabilidad

f_exp <- function(x) {
  exp(b_exp * x)
}

C_exp <- integrate(
  f_exp,
  lower = -50,
  upper = -20
)$value

f_exp_norm <- function(x) {

  f_exp(x) / C_exp

}

Gráfica de densidad

x <- seq(
  -50,
  -20,
  length.out = 500
)

y <- f_exp_norm(x)

plot(
  x,
  y,
  type = "l",
  lwd = 2,
  col = "skyblue3",
  main = "Gráfica 5: Cálculo de probabilidades del modelo exponencial",
  xlab = "Latitud (°)",
  ylab = "Densidad de probabilidad"
)

x_sec <- seq(
  -35,
  -25,
  by = 0.01
)

y_sec <- f_exp_norm(x_sec)

polygon(
  c(x_sec, rev(x_sec)),
  c(y_sec, rep(0, length(y_sec))),
  col = rgb(1, 0, 0, 0.6),
  border = NA
)

lines(
  x_sec,
  y_sec,
  col = "red",
  lwd = 2
)

legend(
  "topleft",
  legend = c(
    "Modelo exponencial",
    "Área de probabilidad"
  ),
  col = c(
    "skyblue3",
    "red"
  ),
  lwd = 4,
  bty = "n"
)

6.2 Cálculo de probabilidades

prob_decimal <- integrate(
  f_exp_norm,
  lower = -35,
  upper = -25
)$value

prob_porcentaje <- round(
  prob_decimal * 100,
  1
)

paste(
  "La probabilidad es de:",
  prob_porcentaje,
  "%"
)
## [1] "La probabilidad es de: 39.3 %"