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

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
retornoV <- na.omit(retorno)

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

## [1] 898
## [1] 11
## [1] 126.7845
# CALCULO ni y hi
# Frecuencia absoluta (ni)
ni <- numeric(k)

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

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

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

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

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

# Mostrar tabla
TDFretorno
##              Intervalo  ni     hi
## 1           [7.23,134] 483  53.79
## 2            (134,261] 220  24.50
## 3            (261,388]  96  10.69
## 4            (388,514]  35   3.90
## 5            (514,641]  24   2.67
## 6            (641,768]  11   1.22
## 7            (768,895]  12   1.34
## 8       (895,1.02e+03]   5   0.56
## 9  (1.02e+03,1.15e+03]   4   0.45
## 10 (1.15e+03,1.28e+03]   2   0.22
## 11  (1.28e+03,1.4e+03]   6   0.67
## 12               TOTAL 898 100.00

TABLA DE FRECUENCIAS

Distribución de Frecuencias del Período de Retorno de los Volcanes
Análisis de frecuencias globales según los intervalos del período de retorno
Intervalo del Período de Retorno (años) Frecuencia Absoluta (ni) Frecuencia Relativa (hi %)
[7.23,134] 483 53.79
(134,261] 220 24.50
(261,388] 96 10.69
(388,514] 35 3.90
(514,641] 24 2.67
(641,768] 11 1.22
(768,895] 12 1.34
(895,1.02e+03] 5 0.56
(1.02e+03,1.15e+03] 4 0.45
(1.15e+03,1.28e+03] 2 0.22
(1.28e+03,1.4e+03] 6 0.67
TOTAL 898 100.00

GRÁFICA DE DISTRIBUCIÓN DE PROBABILIDAD

DIAGRAMA DE DISTRIBUCIÓN DE FRECUENCIA ABSOLUTA 

hist(
  retorno,
  breaks = breaks,
  freq = TRUE,
  xaxt = "n",
  col = "lightblue",
  border = "black",
  main = "Histograma de frecuencias absolutas",
  xlab = "Periodo de retorno (años)",
  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 LOGNORMAL

CALCULO DE PARAMETROS 

# PARÁMETRO μ (media del logaritmo)

medialog <- mean(log(retorno))
medialog
## [1] 4.790609
# PARÁMETRO σ (desviación estándar del logaritmo)

sd_log <- sd(log(retorno))
sd_log
## [1] 0.9695931

MODELO Y REALIDAD 

#1
# GRÁFICA
par(oma = c(1, 1, 1, 1))

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

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

histograma <- hist(
  retorno,
  plot = FALSE,
  breaks = breaks
)

h <- length(histograma$counts)
h
## [1] 11
curve(
  dlnorm(
    x,
    meanlog = medialog,
    sdlog = sd_log
  ),
  from = min(retorno),
  to = max(retorno),
  add = TRUE,
  col = "red",
  lwd = 3
)

FRECUENCIAS OBSERVADAS Y ESPERADAS 

# Tamaño de la muestra completa
n <- length(retorno)

# Frecuencias observadas
Fo <- histograma$counts

# Número de intervalos
h <- length(Fo)


# Probabilidades esperadas según el modelo Lognormal

P <- numeric(h)

for (i in 1:h) {
  
  P[i] <- 
    plnorm(
      histograma$breaks[i+1],
      meanlog = medialog,
      sdlog = sd_log
    ) -
    plnorm(
      histograma$breaks[i],
      meanlog = medialog,
      sdlog = sd_log
    )
}


# Frecuencias esperadas
Fe <- P * n


# Conversión a porcentajes

# Frecuencia observada (%)
Fo_perc <- (Fo / n) * 100

Fo_perc
##  [1] 53.7861915 24.4988864 10.6904232  3.8975501  2.6726058  1.2249443
##  [7]  1.3363029  0.5567929  0.4454343  0.2227171  0.6681514
# Frecuencia esperada (%)
Fe_perc <- (Fe / n) * 100

Fe_perc
##  [1] 54.2204257 24.3319094  9.8705382  4.6828299  2.4829707  1.4265085
##  [7]  0.8706761  0.5572053  0.3704805  0.2542198  0.1791310

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

# Gráfica
plot(
  Fo_perc,
  Fe_perc,
  main = "Gráfica N°3. Correlación de frecuencias en el modelo Lognormal del período de retorno",
  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 

# Correlación entre realidad y modelo Gamma
Correlacion <- cor(Fo, Fe) * 100
Correlacion
## [1] 99.9664
# 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 

# Evitar divisiones por cero
Fe[Fe == 0] <- 1e-9

# Estadístico Chi-cuadrado
x2 <- sum((Fo - Fe)^2 / Fe)

x2
## [1] 16.61715
# Umbral de aceptación
umbral_aceptacion <- qchisq(0.99, df = k - 3)
umbral_aceptacion
## [1] 20.09024
# Evaluación del Chi-cuadrado

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

TABLA DE RESUMEN

Tabla. Resumen del test de bondad de ajuste del modelo Lognormal
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Período de retorno de los volcanes (años) 99.97 16.62 20.09

CALCULO DE PROBABILIDAD

PREGUNTA

CALCULO DE VALORES 

# Aquí trabajamos con todo nuestro dataset para el cálculo de los parámetros

n_completo <- length(retorno)
n_completo
## [1] 898
# PARÁMETRO MEDIA DEL LOGARITMO (μ)

# Calculamos μ con todos los datos globales
medialog_completo <- mean(log(retorno))

medialog_completo
## [1] 4.790609
# PARÁMETRO DESVIACIÓN ESTÁNDAR DEL LOGARITMO (σ)
sd_log_completo <- sd(log(retorno))
sd_log_completo
## [1] 0.9695931

CALCULO DE PROBABILIDAD 

# USAMOS LOS PARÁMETROS COMPLETOS RECALCULADOS

probabilidad_retorno <- plnorm(
  200,
  meanlog = medialog_completo,
  sdlog = sd_log_completo
) -
plnorm(
  50,
  meanlog = medialog_completo,
  sdlog = sd_log_completo
)

# PROBABILIDAD (%)

probabilidad_retorno * 100
## [1] 51.7301

AREA DE PROBABILIDAD

INTERVALOS DE CONFIANZA

## [1] 192.6094
## [1] 240.5895
## [1] 8.028575
## [1] 176.5522
## [1] 208.6665

TABLA RESUMEN

Tabla Nro.3: Intervalo de confianza de la media poblacional del período de retorno de los volcanes activos
Límite inferior Media poblacional Límite superior Desviación estándar poblacional
176.55 Período de retorno de los volcanes activos (años) 208.67 240.59

CONCLUSIÓN