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
H_pl <- Volcanes_Globales$est_plume_height_km

# Limpieza de datos
H_pl_limpios <- H_pl[!is.na(H_pl)]

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

## [1] 714
## [1] 7.98
# CALCULO ni y hi
# Frecuencia absoluta ni
ni <- numeric(k)

for (i in 1:k) {
  
  if (i < k) {
    
    ni[i] <- sum(
      H_pl_limpios >= Li[i] &
        H_pl_limpios < Ls[i],
      na.rm = TRUE
    )
    
  } else {
    
    ni[i] <- sum(
      H_pl_limpios >= Li[i] &
        H_pl_limpios <= maximo,
      na.rm = TRUE
    )
  }
}

# Frecuencia relativa hi
n <- length(H_pl_limpios)
hi <- ni / n

# TABLA DE DISTRIBUCION DE FRECUENCIAS
TDF_H_pl <- data.frame(
  Intervalo = levels(Intervalo),
  ni = ni,
  hi = round(hi,4)*100
)

# Fila de Total
fila_total <- data.frame(
  Intervalo = "TOTAL",
  ni = sum(TDF_H_pl$ni),
  hi = 100
)

# Unir fila Total con la tabla
TDF_H_pl_final <- rbind(
  TDF_H_pl,
  fila_total
)

# Mostrar tabla
TDF_H_pl_final
##     Intervalo  ni     hi
## 1  [0.1,8.08] 530  74.23
## 2 (8.08,16.1] 114  15.97
## 3   (16.1,24]  30   4.20
## 4     (24,32]  36   5.04
## 5     (32,40]   4   0.56
## 6       TOTAL 714 100.00

TABLA DE FRECUENCIAS

Distribución de Frecuencias de la Altura de Columna Eruptiva de los Volcanes Activos
Análisis de frecuencias globales según los intervalos de altura de columna registrados
Intervalo de Altura de Pluma (km) Frecuencia Absoluta (ni) Frecuencia Relativa (hi %)
[0.1,8.08] 530 74.23
(8.08,16.1] 114 15.97
(16.1,24] 30 4.20
(24,32] 36 5.04
(32,40] 4 0.56
TOTAL 714 100.00

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE FRECUENCIA ABSOLUTA 

hist(
  H_pl_limpios,
  breaks = breaks,
  freq = TRUE,
  xaxt = "n",
  col = "lightblue",
  border = "black",
  main = "Histograma de frecuencias absolutas de la altura de columna eruptiva",
  xlab = "Altura de columna eruptiva (km)",
  ylab = "Cantidad"
)

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

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL SE CONJETURA UN MODELO EXPONENCIAL

CALCULO DE PARAMETROS 

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

# Datos agrupados
x <- TDF_H_pl_final$ni[TDF_H_pl_final$Intervalo != "TOTAL"]

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

# Media observada mediante frecuencias agrupadas
media_observada <- sum(marcas_clase * x) / n

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

media_observada
## [1] 7.420588
lambda
## [1] 0.1347602

MODELO Y REALIDAD 

# FRECUENCIAS OBSERVADAS
Fo <- TDF_H_pl_final$hi[TDF_H_pl_final$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
    )
}

Fe <- P_exp * 100

Fo
## [1] 74.23 15.97  4.20  5.04  0.56
Fe
## [1] 65.0015205 22.1763015  7.5657976  2.5811921  0.8806147
# GRÁFICA
par(mar = c(5,5,5,2))

barplot(
  rbind(Fo, Fe),
  beside = TRUE,
  col = c("skyblue", "blue"),
  names.arg = TDF_H_pl_final$Intervalo[TDF_H_pl_final$Intervalo != "TOTAL"],
  main = "Gráfica N°2: Comparación de la realidad con el modelo exponencial de la altura de columna eruptiva",
  ylab = "Probabilidad (%)",
  xlab = "Intervalos de altura de columna (km)",
  ylim = c(0, max(c(Fo, Fe)) * 1.3),
  las = 2
)

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

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

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

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

# Gráfica
plot(
  Fo,
  Fe,
  main = "Gráfica N°3: Correlación de frecuencias entre la realidad y el modelo exponencial",
  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] 98.91383
# 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 <- TDF_H_pl_final$ni[TDF_H_pl_final$Intervalo != "TOTAL"]

# Tamaño de la muestra
n <- sum(Fo_abs)

# Frecuencias absolutas esperadas
Fe_abs <- P_exp * n

# Grados de libertad
grados_libertad <- length(Fo_abs) - 2
grados_libertad
## [1] 3
# Nivel de significancia
nivel_significancia <- 0.95

# Estadístico Chi-cuadrado
x2 <- sum((Fo_abs - Fe_abs)^2 / Fe_abs)
x2
## [1] 50.03373
## [1] 54.23451
# Evaluación
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

Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Elevación de los volcanes activos 98.91 50.03 54.23

CALCULO DE PROBABILIDAD

PREGUNTA

CALCULO DE VALORES 

#Aqui ya trabajamos con todo nuestro dataset para el calculo de los parámetros
H_pl_completo <- as.numeric(Volcanes_Globales$est_plume_height_km)
H_pl_completo <- na.omit(H_pl_completo)

MEDIA ARITMETICA 

media_completa <- mean(H_pl_completo)
media_completa
## [1] 5.466947
# Recalculamos los parámetros reales de toda la población
n_completo <- length(H_pl_completo)
n_completo
## [1] 714
# Lambda Completo
lambda_completo <- 1 / media_completa
lambda_completo
## [1] 0.1829175

CALCULO DE PROBABILIDAD 

# USAMOS LAMBDA COMPLETO
probabilidad_H_pl <- pexp(32, lambda_completo) - pexp(24, lambda_completo)

# Probabilidad:
probabilidad_H_pl * 100
## [1] 0.9530251

AREA DE PROBABILIDAD

INTERVALOS DE CONFIANZA

## [1] 5.466947
## [1] 7.775175
## [1] 714
## [1] 0.2909786
## [1] 4.88499
## [1] 6.048904
Tabla Nro.3: Media poblacional de la altura de columna eruptiva de los volcanes activos
Limite inferior Media poblacional Limite superior Desviación estandar poblacional
4.88 Altura de columna eruptiva de los volcanes activos 6.05 0.29

CONCLUSIÓN