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] 887
## [1] 7
## [1] 1000
# 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 <= Ls[i],
                 na.rm = TRUE)
  }
}

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


# TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

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

# 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    [0 - 1000] 196  22.10
## 2 (1000 - 2000] 329  37.09
## 3 (2000 - 3000] 201  22.66
## 4 (3000 - 4000] 110  12.40
## 5 (4000 - 5000]  17   1.92
## 6 (5000 - 6000]  32   3.61
## 7 (6000 - 7000]   2   0.23
## 8         TOTAL 887 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 %)
[0 - 1000] 196 22.10
(1000 - 2000] 329 37.09
(2000 - 3000] 201 22.66
(3000 - 4000] 110 12.40
(4000 - 5000] 17 1.92
(5000 - 6000] 32 3.61
(6000 - 7000] 2 0.23
TOTAL 887 100.00

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE FRECUENCIA ABSOLUTA 

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


# Eje X con intervalos de 500 m
axis(
  side = 1,
  at = breaks_elevacion,
  labels = round(breaks_elevacion, 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] 1966.742
varianza_observada
## [1] 1533599
shape
## [1] 2.522219
scale
## [1] 779.7664

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.10 37.09 22.66 12.40  1.92  3.61  0.23
Fe
## [1] 22.829499 36.573447 22.852382 10.750304  4.411982  1.670899  0.600152
# Frecuencias observadas (%)
Fo <- TDFelevacion$hi[TDFelevacion$Intervalo != "TOTAL"]

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

GRAFICA 

## [1] 7

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

#GRAFICA CORRELACION
# 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 entre 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.38871
# 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] 4
# Nivel de significancia
nivel_significancia <- 0.9999999

# Estadístico Chi-cuadrado
x2 <- sum((Fo_abs - Fe_abs)^2 / Fe_abs)
x2
## [1] 37.04345
# Valor crítico
umbral_aceptacion <- qchisq(
  nivel_significancia,
  grados_libertad
)

umbral_aceptacion
## [1] 38.2396
# 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.39 37.04 38.24

CALCULO DE PROBABILIDAD

PREGUNTA

LIMITES INFERIOR Y SUPERIOR 

# Límites del intervalo
limite_inferior <- 2000
limite_superior <- 3000

MEDIA ARITMETICA 

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

CALCULO DE PROBABILIDAD 

# Cálculo de probabilidad mediante distribución Gamma
probabilidad_elevacion <- pgamma(
  limite_superior,
  shape = shape,
  scale = scale
) -
  pgamma(
    limite_inferior,
    shape = shape,
    scale = scale
  )

# Resultado en porcentaje
probabilidad_elevacion_porcentaje <- probabilidad_elevacion * 100
probabilidad_elevacion_porcentaje
## [1] 22.85238

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] 68.55715

INTERVALOS DE CONFIANZA

#Media aritmetica
x <- mean(elevacion)
x
## [1] 1985.594
#Desviación estandar
sigma <- sd(elevacion)
sigma
## [1] 1237.792
#Tamaño muestral
n <- length(elevacion)
n
## [1] 887
#P(x-2e<u<x+2e)=95%
e <- sigma/sqrt(n)
e
## [1] 41.56098
#Limite inferior
li <- x - 2*e
li
## [1] 1902.472
#Limite superior
ls <- x + 2*e
ls
## [1] 2068.716

TABLA

Tabla Nro.3: Media poblacional de la elevación de los volcanes activos
Limite inferior Media poblacional Limite superior Desviación estandar
1902.47 Elevación de los volcanes activos 2068.72 41.56

CONCLUSIÓN