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

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
elevacionV <- Volcanes_Globales$elevation_m

elevacion <- elevacionV[elevacionV > 0]

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

## [1] 885
## [1] 6
## [1] 1061
# CALCULO ni y hi
# Frecuencia absoluta ni
ni <- numeric(k)


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

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

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

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

# Unir fila Total con la tabla
TDFelevacion <- rbind(TDF, fila_total)

# Mostrar tabla
TDFelevacion
##             Intervalo  ni     hi
## 1       [11,1.07e+03] 200  22.60
## 2 (1.07e+03,2.13e+03] 338  38.19
## 3 (2.13e+03,3.19e+03] 199  22.49
## 4 (3.19e+03,4.26e+03] 103  11.64
## 5 (4.26e+03,5.32e+03]  30   3.39
## 6 (5.32e+03,6.38e+03]  15   1.69
## 7               TOTAL 885 100.00

TABLA DE FRECUENCIAS

Distribución de Frecuencias de la Elevación de los Volcanes Activos
Análisis de frecuencias globales según los intervalos de elevación registrados
Intervalo de Elevación (m) Frecuencia Absoluta (ni) Frecuencia Relativa (hi %)
[11,1.07e+03] 200 22.60
(1.07e+03,2.13e+03] 338 38.19
(2.13e+03,3.19e+03] 199 22.49
(3.19e+03,4.26e+03] 103 11.64
(4.26e+03,5.32e+03] 30 3.39
(5.32e+03,6.38e+03] 15 1.69
TOTAL 885 100.00

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE FRECUENCIA ABSOLUTA 

hist(
  elevacion,
  breaks = breaks,
  freq = TRUE,
  xaxt = "n",
  col = "lightblue",
  border = "black",
  main = "Histograma de frecuencias absolutas",
  xlab = "Elevación (m.s.n.m)",
  ylab = "Cantidad"
)

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

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL SE CONJETURA UN MODELO GAMMA

CALCULO DE PARAMETROS 

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

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

# Media observada mediante frecuencias
marcas_clase <- (Li + Ls) / 2
media_observada <- sum(marcas_clase * x) / n

# Varianza observada mediante datos agrupados
varianza_observada <- sum(x * (marcas_clase - media_observada)^2) / (n-1)

# Parámetros Gamma método de momentos
shape <- media_observada^2 / varianza_observada
scale <- varianza_observada / media_observada

media_observada
## [1] 2028.099
varianza_observada
## [1] 1500875
shape
## [1] 2.740525
scale
## [1] 740.0404

MODELO Y REALIDAD 

# FRECUENCIAS OBSERVADAS
Fo <- TDFelevacion$hi[TDFelevacion$Intervalo != "TOTAL"]


# FRECUENCIAS ESPERADAS MODELO GAMMA
P_gamma <- numeric(length(Li))

for(i in 1:length(Li)){
  
  P_gamma[i] <- pgamma(
    Ls[i],
    shape = shape,
    scale = scale
  ) -
    pgamma(
      Li[i],
      shape = shape,
      scale = scale
    )
}

Fe <- P_gamma * 100

Fo
## [1] 22.60 38.19 22.49 11.64  3.39  1.69
Fe
## [1] 22.884788 38.377994 23.069256 10.053645  3.750561  1.276586
# GRÁFICA
par(mar = c(5,5,5,2))
barplot(
  rbind(Fo, Fe),
  beside = TRUE,
  col = c("skyblue", "blue"),
  names.arg = TDFelevacion$Intervalo[TDFelevacion$Intervalo != "TOTAL"],
  main = "Gráfica N°2: Comparación de la realidad con el\nmodelo Gamma de la elevación de los volcanes activos",
  ylab = "Probabilidad (%)",
  xlab = "Intervalos de elevación (m)",
  ylim = c(0, max(c(Fo,Fe))*1.3),
  las = 2
)


# Leyenda

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

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

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

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

# Gráfica
plot(
  Fo,
  Fe,
  main = "Gráfica N°3: Correlación de frecuencias\nentre la realidad y el modelo Gamma",
  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 Gamma
Correlacion <- cor(Fo, Fe) * 100
Correlacion
## [1] 99.85691
# 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 <- TDFelevacion$ni[TDFelevacion$Intervalo != "TOTAL"]

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

# Frecuencias absolutas esperadas
Fe_abs <- P_gamma * n

# Grados de libertad
grados_libertad <- length(Fo_abs) - 3
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] 3.901212
# Valor crítico
umbral_aceptacion <- qchisq(
  nivel_significancia,
  grados_libertad
)

umbral_aceptacion
## [1] 7.814728
# 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 99.86 3.9 7.81

CALCULO DE PROBABILIDAD

PREGUNTA

LIMITES INFERIOR Y SUPERIOR 

limite_inferior <- 2133
limite_superior <- 3254

MEDIA ARITMETICA 

media_completa  <- mean(elevacion)
media_completa
## [1] 1990.081

CALCULO DE PROBABILIDAD 

# Probabilidad Gamma
probabilidad_elevacion <- pgamma(
  limite_superior,
  shape = shape,
  scale = scale
) -
  pgamma(
    limite_inferior,
    shape = shape,
    scale = scale
  )

# En porcentaje
probabilidad_elevacion * 100
## [1] 23.91324

AREA DE PROBABILIDAD

CALCULO DE PROBABILIDAD 2

NUEVAS MEDICIONES 

# Número de nuevas mediciones
n_nuevas <- 300

# Probabilidad obtenida anteriormente
p <- probabilidad_elevacion

# Número esperado de volcanes dentro del intervalo
volcanes_esperados <- n_nuevas * p

volcanes_esperados
## [1] 71.73971

INTERVALOS DE CONFIANZA

## [1] 1990.081
## [1] 1235.579
## [1] 885
## [1] 41.53352
## [1] 1907.014
## [1] 2073.148
Tabla Nro.3: Media poblacional de la elevación de los volcanes activos
Limite inferior Media poblacional Limite superior Desviación estandar poblacional
1907.01 Elevación de los volcanes activos 2073.15 41.53

CONCLUSIÓN