CARGA DE LIBRERIAS Y DATOS

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#CARGA DE DATASET
Volcanes_Globales <- read.csv("global_volcano_eruption_intelligence.csv", header = T, sep = ";", dec = ".")

#SELECCION VARIABLE
risk_vei <- Volcanes_Globales$vei_risk_tier

#LIMPIEZA DE DATOS
sum(is.na(risk_vei))

## [1] 0

risk_vei <- risk_vei[risk_vei != "Unknown"]

ORDEN DE CATEGORIAS

# Orden según la jerarquía (sin Unknown)
niveles_riesgo <- c(
  "Minimal",
  "Low",
  "Moderate",
  "High",
  "Very High",
  "Extreme",
  "Catastrophic"
)

risk_vei <- factor(
  risk_vei,
  levels = niveles_riesgo,
  ordered = TRUE
)

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

CALCULO DE FRECUENCIAS

# Generar tabla de frecuencias
TDFRiesgo <- table(risk_vei)

# Convertir datos a dataframe
TDFRiesgo <- as.data.frame(TDFRiesgo)

# Cambiar nombres
colnames(TDFRiesgo) <- c("Nivel_Riesgo", "Freq")

# Frecuencia absoluta (ni) y frecuencia relativa (hi)
TDFRiesgoFinal <- TDFRiesgo %>%
  group_by(Nivel_Riesgo) %>%
  summarise(
    ni = sum(Freq),
    
    # Frecuencia relativa porcentual
    hi = round((ni / sum(TDFRiesgo$Freq)) * 100, 2),
    
    # Frecuencia relativa decimal
    hi_decimal = round(ni / sum(TDFRiesgo$Freq), 4)
  )

# Agregar numeración en tabla
TDFRiesgoFinal <- TDFRiesgoFinal %>%
  mutate(Nro = row_number())

# Añadir fila de totales
TDFRiesgoFinal <- TDFRiesgoFinal %>%
  add_row(
    Nro = NA,
    Nivel_Riesgo = "Total",
    ni = sum(TDFRiesgoFinal$ni),
    hi = 100,
    hi_decimal = 1
  )

# Ver tabla
TDFRiesgoFinal

## # A tibble: 8 × 5
##   Nivel_Riesgo    ni     hi hi_decimal   Nro
##   <chr>        <int>  <dbl>      <dbl> <int>
## 1 Minimal         68   9.52     0.0952     1
## 2 Low            240  33.6      0.336      2
## 3 Moderate       222  31.1      0.311      3
## 4 High           114  16.0      0.160      4
## 5 Very High       30   4.2      0.042      5
## 6 Extreme         36   5.04     0.0504     6
## 7 Catastrophic     4   0.56     0.0056     7
## 8 Total          714 100        1         NA

TABLA DE FRECUENCIAS Y PROBABILIDAD

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE DENSIDAD DE PROBABILIDAD

par(mar = c(9,4,4,2))
barplot(
  TDFRiesgoFinal$hi[1:(nrow(TDFRiesgoFinal)-1)],
  names.arg = TDFRiesgoFinal$Nro[1:(nrow(TDFRiesgoFinal)-1)],
  col = "colores_riesgo" <- c(
    "#008000",  # Minimal
    "#66CC66",  # Low
    "#FFFF66",  # Moderate
    "#FFCC66",  # High
    "#FF9933",  # Very High
    "#FF6666",  # Extreme
    "#CC0000"   # Catastrophic 
  ),
  ylim = c(0, 100),
  main = "Gráfica N°1: Distribución de probabilidad del nivel de riesgo VEI",
  ylab = "Probabilidad (%)",
  xlab = ""
)

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL SE CONJETURA UN MODELO GAMMA

CALCULO DE PARAMETROS

#CALCULO DE PARAMETROS
n <- sum(TDFRiesgoFinal$ni)

x <- TDFRiesgoFinal$ni

X <- 1:length(x)

media_observada <- sum(X*x)/n

p <- media_observada/length(x)

P_binomial <- dbinom(X, size=length(x), prob=p)

P_binomial

## [1] 0.001843676 0.013754790 0.058638841 0.156241648 0.266433126 0.283961621
## [7] 0.172939032 0.046079150

MODELO Y REALIDAD

n <- sum(TDFRiesgoFinal$ni[1:(nrow(TDFRiesgoFinal)-1)])

x <- TDFRiesgoFinal$ni[1:(nrow(TDFRiesgoFinal)-1)]

X <- 1:length(x)

media_observada <- sum(X*x)/n

p <- media_observada/length(x)

P_binomial <- dbinom(X, size=length(x), prob=p)

Fo <- TDFRiesgoFinal$hi[1:(nrow(TDFRiesgoFinal)-1)]

Fe <- P_binomial*100

# GRÁFICA
par(mar = c(5,5,5,2))
barplot(
  rbind(Fo, Fe),
  beside = TRUE,
  col = c("skyblue", "blue"),
  names.arg = TDFRiesgoFinal$Nro[1:(nrow(TDFRiesgoFinal)-1)],
  main = "Gráfica N°2: Comparación de la realidad con el\nmodelo binomial del nivel de riesgo VEI",
  ylab = "Probabilidad (%)",
  xlab = "Nivel de riesgo (Nro)",
  ylim = c(0,100)
)

# Leyenda
legend(
  "topright",
  legend = c("Realidad", "Modelo binomial"),
  fill = c("skyblue", "blue"),
  border = "black",
  bty = "n"
)

GRÁFICA DE CORRELACIÓN DEL MODELO BINOMIAL Y LA REALIDAD

plot(Fo, Fe,
     main = "Gráfica N°3: Correlación de frecuencias en el modelo binomial del nivel de riesgo VEI",
     xlab = "Frecuencia Observada (%)",
     ylab = "Frecuencia Esperada (%)",
     pch = 19,
     col = "darkblue",
     xlim = c(0,100),
     ylim = c(0,100))

# Línea de ajuste perfecto
abline(a = 0, b = 1, col = "red", lwd = 2)

TESTS DE APROBACIÓN

PEARSON

Correlacion <- cor(Fo, Fe) * 100
Correlacion

## [1] 93.91243

CHI CUADRADO

# Grados de libertad
gl <- length(Fo)-1
# Chi cuadrado calculado
x2 <- sum((Fo-Fe)^2/Fe)
# Valor crítico al 99%
vc <- qchisq(0.99, gl)
# Resultados
x2

## [1] 11.95202

vc

## [1] 16.81189

x2 < vc

## [1] TRUE

#RESUMEN
library(knitr)
tabla_resumen <- data.frame(
  Variable="Nivel de Riesgo VEI",
  Pearson=round(Correlacion,2),
  Chi2=round(x2,2),
  Umbral=round(vc,2)
)

TABLA DE RESUMEN

CALCULO DE PROBABILIDAD

# 1. Calculamos el valor de la probabilidad dinámicamente
# Categoría Extreme (riesgo extremo)

porcentaje <- round(
  (TDFRiesgoFinal$ni[TDFRiesgoFinal$Nivel_Riesgo == "Extreme"] /
     TDFRiesgoFinal$ni[TDFRiesgoFinal$Nivel_Riesgo == "Total"]) * 100,
  1
)

porcentaje

## [1] 5

CONCLUSIÓN