ANÁLISIS ESTADÍSTICO DE VOLCANES ACTIVOS A NIVEL GLOBAL
MODELO EXPONENCIAL

CARRERA DE GEOLOGÍA

GRUPO N°2

ANÁLISIS ESTADÍSTICO INFERENCIAL

1. CARGA DE LIBRERIAS Y DATOS

2. TABLA DE DISTRIBUCIÓN DE PROBABILIDAD

3. GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

4. CONJETURA DEL MODELO

5. TEST DE APROBACIÓN

6. CALCULO DE PROBABILIDAD

7. INTERVALOS DE CONFIANZA

8. CONCLUSIÓN

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
retornoV<- Volcanes_Globales$avg_eruption_return_period_years

# Eliminar NA
retorno <- na.omit(retornoV)

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

## [1] 898
## [1] 8
## [1] 200
# Frecuencia absoluta
ni <- numeric(k)

for (i in 1:k) {
  
  if (i == 1) {
    
    # Primer intervalo: [0, 200]
    ni[i] <- sum(
      retorno >= Li[i] &
        retorno <= Ls[i],
      na.rm = TRUE
    )
    
  } else {
    
    # Resto de intervalos: (Li, Ls]
    ni[i] <- sum(
      retorno > Li[i] &
        retorno <= Ls[i],
      na.rm = TRUE
    )
    
  }
}

# Frecuencia relativa
hi <- ni / n

# ==============================================================================
# TABLA DE DISTRIBUCIÓN DE FRECUENCIAS
# ==============================================================================

TDFretorno <- data.frame(
  Intervalo = etiquetas,
  ni = ni,
  hi = round(hi * 100, 2)
)

# Agregar fila TOTAL
TDFretorno <- rbind(
  TDFretorno,
  data.frame(
    Intervalo = "TOTAL",
    ni = sum(ni),
    hi = 100
  )
)

# Mostrar tabla
TDFretorno
##       Intervalo  ni     hi
## 1     [0 - 200] 626  69.71
## 2   (200 - 400] 176  19.60
## 3   (400 - 600]  46   5.12
## 4   (600 - 800]  24   2.67
## 5  (800 - 1000]  13   1.45
## 6 (1000 - 1200]   6   0.67
## 7 (1200 - 1400]   6   0.67
## 8 (1400 - 1600]   1   0.11
## 9         TOTAL 898 100.00

TABLA DE FRECUENCIAS

Distribución de Frecuencias del Período de Retorno de los Volcanes
Análisis de frecuencias mediante ocho intervalos de amplitud de 200 años
Intervalo del Período de Retorno (años) Frecuencia Absoluta (ni) Frecuencia Relativa (hi) (%)
[0 - 200] 626 69.71
(200 - 400] 176 19.60
(400 - 600] 46 5.12
(600 - 800] 24 2.67
(800 - 1000] 13 1.45
(1000 - 1200] 6 0.67
(1200 - 1400] 6 0.67
(1400 - 1600] 1 0.11
TOTAL 898 100.00

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE FRECUENCIA ABSOLUTA 

hist(
  retorno,
  breaks = breaks_retorno,
  freq = TRUE,
  xaxt = "n",
  col = "lightblue",
  border = "black",
  main = "Histograma de Frecuencias Absolutas",
  xlab = "Período de Retorno (años)",
  ylab = "Cantidad"
)

axis(
  side = 1,
  at = breaks_retorno,
  labels = round(breaks_retorno, 0),
  las = 2
)

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL SE CONJETURA UN MODELO EXPONENCIAL

CALCULO DE PARAMETROS 

#PARAMETROS
# Número total de volcanes
n <- sum(TDFretorno$ni[TDFretorno$Intervalo != "TOTAL"])

# Datos agrupados (frecuencias absolutas)
x <- TDFretorno$ni[TDFretorno$Intervalo != "TOTAL"]

# Marcas de clase
marcas_clase <- (Li + Ls) / 2

# Media observada mediante datos agrupados
media_observada <- sum(marcas_clase * x) / n

# Parámetro de la distribución Exponencial LAMBDA
lambda <- 1 / media_observada

# Mostrar resultados
media_observada
## [1] 203.5635
lambda
## [1] 0.004912473

MODELO Y REALIDAD 

#==============================================================================
# FRECUENCIAS OBSERVADAS
# ==============================================================================

Fo <- TDFretorno$hi[TDFretorno$Intervalo != "TOTAL"]


# ==============================================================================
# FRECUENCIAS ESPERADAS - MODELO EXPONENCIAL
# ==============================================================================

P_exp <- numeric(length(Li))

for(i in 1:length(Li)){
  
  P_exp[i] <- pexp(
    Ls[i],
    rate = lambda
  ) -
    pexp(
      Li[i],
      rate = lambda
    )
}

# Frecuencia esperada porcentual
Fe <- P_exp * 100


# Mostrar resultados
Fo
## [1] 69.71 19.60  5.12  2.67  1.45  0.67  0.67  0.11
Fe
## [1] 62.56239590 23.42186209  8.76858400  3.28274776  1.22898211  0.46010146
## [7]  0.17225096  0.06448663
# GRÁFICA
par(oma = c(1, 1, 1, 1))

hist(
  retorno,
  freq = FALSE,
  breaks = breaks_retorno,
  main = "Gráfica Nº3. Comparación de la realidad con el modelo Exponencial del período de retorno volcánico",
  xlab = "Período de retorno (años)",
  ylab = "Densidad de probabilidad",
  col = "cyan",
  border = "gray30"
)

box(which = "outer", col = "black")


# Histograma sin graficar para obtener frecuencias
histograma <- hist(
  retorno,
  plot = FALSE,
  breaks = breaks_retorno
)

h <- length(histograma$counts)
h
## [1] 8
# Curva del modelo exponencial
curve(
  dexp(
    x,
    rate = lambda
  ),
  from = min(retorno),
  to = max(retorno),
  add = TRUE,
  col = "red",
  lwd = 3
)

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

# Frecuencias observadas (%)
Fo <- TDFretorno$hi[TDFretorno$Intervalo != "TOTAL"]

# Frecuencias esperadas (%)
Fe <- P_exp * 100


# Gráfica
plot(
  Fo,
  Fe,
  main = "Gráfica N°4: Correlación de frecuencias\nentre la realidad y el modelo exponencial del período de retorno",
  xlab = "Frecuencia observada (%)",
  ylab = "Frecuencia esperada (%)",
  pch = 19,
  col = "darkblue",
  xlim = c(0, max(c(Fo, Fe)) * 1.1),
  ylim = c(0, max(c(Fo, Fe)) * 1.1)
)

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

TESTS DE APROBACIÓN

PEARSON 

# Correlación entre realidad y modelo Exponencial
Correlacion <- cor(Fo, Fe) * 100
Correlacion
## [1] 99.40248
# Evaluación del Test de Pearson (>80%)
if(Correlacion > 80){
  
  print("APRUEBA EL TEST DE PEARSON")
  
}else{
  
  print("NO APRUEBA EL TEST DE PEARSON")
  
}
## [1] "APRUEBA EL TEST DE PEARSON"

CHI CUADRADO 

# Frecuencias absolutas observadas
Fo_abs <- TDFretorno$ni[TDFretorno$Intervalo != "TOTAL"]


# Tamaño de la muestra
n <- sum(Fo_abs)
n
## [1] 898
# Frecuencias absolutas esperadas
Fe_abs <- P_exp * n


# Grados de libertad
# gl = número de intervalos - parámetros estimados - 1
# La distribución exponencial tiene 1 parámetro: lambda

grados_libertad <- length(Fo_abs) - 2
grados_libertad
## [1] 6
# Nivel de confianza
nivel_confianza <- 0.9999999


# Estadístico Chi-cuadrado

x2 <- sum(
  (Fo_abs - Fe_abs)^2 / Fe_abs
)

x2
## [1] 41.89031
# Valor crítico Chi-cuadrado

umbral_aceptacion <- qchisq(
  nivel_confianza,
  df = grados_libertad
)

umbral_aceptacion
## [1] 43.33776
# Evaluación del modelo

if (x2 < umbral_aceptacion) {
  
  print("APRUEBA EL TEST DE CHI-CUADRADO")
  
} else {
  
  print("NO APRUEBA EL TEST DE CHI-CUADRADO")
  
}
## [1] "APRUEBA EL TEST DE CHI-CUADRADO"

TABLA DE RESUMEN

Evaluación estadística del modelo exponencial del período de retorno volcánico
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Período de retorno (años) 99.4 41.89 43.34

CALCULO DE PROBABILIDAD

PREGUNTA

CALCULO DE VALORES 

## Trabajamos con todo nuestro dataset para el cálculo de parámetros

retorno_completo <- as.numeric(retorno)

retorno_completo <- na.omit(retorno_completo)

# Tamaño de la población

n_completo <- length(retorno_completo)

n_completo
## [1] 898

MEDIA ARITMETICA 

media_completa <- mean(retorno_completo)

media_completa
## [1] 189.8882
# ==============================================================================
# Parámetro lambda del modelo exponencial
# ==============================================================================

lambda_completo <- 1 / media_completa

lambda_completo
## [1] 0.005266256

CALCULO DE PROBABILIDAD 

limite_inferior <- 400
limite_superior <- 800


probabilidad_retorno <- pexp(
  limite_superior,
  rate = lambda_completo
) -
  pexp(
    limite_inferior,
    rate = lambda_completo
  )


# Probabilidad en porcentaje

probabilidad_retorno * 100
## [1] 10.68609

AREA DE PROBABILIDAD

INTERVALOS DE CONFIANZA

# Datos completos
retorno_completo <- as.numeric(retorno_completo)
retorno_completo <- na.omit(retorno_completo)

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

# Media muestral (solo usada para calcular el intervalo)
media_muestral <- mean(retorno_completo)
media_muestral
## [1] 189.8882
# Desviación estándar muestral
desviacion <- sd(retorno_completo)
desviacion
## [1] 209.0615
# Error estándar
e <- desviacion / sqrt(n)
e
## [1] 6.976474
# Límites del intervalo de confianza aproximado 95%
li <- media_muestral - 2 * e
ls <- media_muestral + 2 * e
Tabla Nro.3: Intervalo de confianza de la media poblacional del período de retorno volcánico
Límite inferior Media poblacional Límite superior Error estándar
175.94 μ 203.84 6.98

CONCLUSIÓN