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. SELECCION DE VARIABLE

3. TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

4. GRÁFICA DE DISTRIBUCION DE FRECUENCIA

5. CONJETURA DEL MODELO

6. PARAMETROS

7. MODELO Y REALIDAD

8. TESTS DE APROBACIÓN

9. PROBABILIDAD

10. INTERVALOS DE CONFIANZA

11. CONCLUSIONES

CARGA DE LIBRERIAS Y DATOS

Carga de librerias 

## 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

lectura del dataset 

#CARGA DE DATASET
Volcanes_Globales <- read.csv("global_volcano_eruption_intelligence.csv", header = T, sep = ";", dec = ".")

SELECCIÓN DE VARIABLE

# Variable
H_pl <- Volcanes_Globales$est_plume_height_km
  
# Limpieza de datos
H_pl_limpioss <- H_pl[!is.na(H_pl)]

TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

## [1] 898
# ==============================================================================
# CÁLCULO DE ni y hi
# ==============================================================================
ni <- numeric(k)

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

# Frecuencia relativa %
hi <- ni / n

# Tabla de distribución

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

# Agregar total

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

TDF_H_pl_final <- rbind(
  TDF_H_pl,
  fila_total
)

TDF_H_pl_final
##   Intervalo  ni     hi
## 1   [0 - 5]  60   6.68
## 2  (5 - 10] 280  31.18
## 3 (10 - 15] 260  28.95
## 4 (15 - 20] 150  16.70
## 5 (20 - 25]  80   8.91
## 6 (25 - 30]  45   5.01
## 7 (30 - 35]  22   2.45
## 8 (35 - 40]   1   0.11
## 9     TOTAL 898 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 Columna (km) Frecuencia Absoluta (ni) Frecuencia Relativa (hi %)
[0 - 5] 60 6.68
(5 - 10] 280 31.18
(10 - 15] 260 28.95
(15 - 20] 150 16.70
(20 - 25] 80 8.91
(25 - 30] 45 5.01
(30 - 35] 22 2.45
(35 - 40] 1 0.11
TOTAL 898 100.00

GRÁFICA DE DISTRIBUCIÓN DE FRECUENCIA

Diagrama de frecuencia absoluta 

hist(
  H_pl_limpios,
  breaks = breaks_Hpl,
  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"
)

# Eje X con los límites de los intervalos
axis(
  side = 1,
  at = breaks_Hpl,
  labels = breaks_Hpl,
  las = 2
)

CONJETURA DEL MODELO

DEBIDO A LA SIMILITUD VISUAL SE CONJETURA UN MODELO LOGNORMAL

PARAMETROS

Calculo de parametros del modelo 

# Parámetro μ (media del logaritmo)
medialog <- mean(log(H_pl_limpios))
medialog
## [1] 2.419605
# Parámetro σ (desviación estándar del logaritmo)
sd_log <- sd(log(H_pl_limpios))
sd_log
## [1] 0.6327913

MODELO Y REALIDAD

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

hist(
  H_pl_limpios,
  freq = FALSE,
  breaks = breaks_Hpl,
  main = "Gráfica Nº3. Comparación de la realidad con el modelo Log-normal de la altura de columna eruptiva",
  xlab = "Altura de columna eruptiva (km)",
  ylab = "Densidad de probabilidad",
  col = "green",
  border = "gray30"
)


# Marco exterior
box(
  which = "outer",
  col = "black"
)



# Histograma sin graficar para obtener frecuencias

histograma <- hist(
  H_pl_limpios,
  plot = FALSE,
  breaks = breaks_Hpl
)

# Número de intervalos
h <- length(histograma$counts)
h
## [1] 8
# ==============================================================================
# CURVA DEL MODELO LOG-NORMAL
# ==============================================================================
curve(
  dlnorm(
    x,
    meanlog = medialog,
    sdlog = sd_log
  ),
  from = min(H_pl_limpios),
  to = max(H_pl_limpios),
  add = TRUE,
  col = "red",
  lwd = 3
)


legend(
  "topright",
  legend = c(
    "Datos reales",
    "Modelo Log-normal"
  ),
  fill = c("green", "red"),
  border = c("gray30", "red"),
  bty = "n"
)

Frecuencias esperadas y observadas 

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

# Frecuencias observadas
Fo <- histograma$counts

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

# Probabilidades esperadas Log-normal
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 absolutas
Fe <- P * n

# Frecuencias observadas en porcentaje
Fo_perc <- (Fo / n) * 100

# Frecuencias esperadas en porcentaje
Fe_perc <- (Fe / n) * 100

# Mostrar resultados

Fo_perc
## [1]  6.6815145 31.1804009 28.9532294 16.7037862  8.9086860  5.0111359  2.4498886
## [8]  0.1113586
Fe_perc
## [1] 10.021873 32.642490 24.910051 14.296447  7.801268  4.285316  2.408399
## [8]  1.390463

TESTS DE APROBACIÓN

Gráfica de correlación del modelo Lognormal y la realidad 

# Frecuencias observadas (%)
Fo <- Fo_perc

# Frecuencias esperadas Log-normal (%)
Fe <- Fe_perc

# Gráfica de correlación
plot(
  Fo,
  Fe,
  main = "Gráfica N°4: Correlación de frecuencias\nentre la realidad y el modelo Log-normal",
  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
)

Test de Pearson 

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

Test Chi cuadrado 

# Frecuencias observadas absolutas
Fo <- histograma$counts

# Frecuencias esperadas absolutas
Fe <- P * n

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

# Estadístico Chi-cuadrado
x2 <- sum(
  (Fo - Fe)^2 / Fe
)
x2
## [1] 33.20792
# ==============================================================================
# GRADOS DE LIBERTAD
# ==============================================================================
# k = número de intervalos
# Log-normal estima 2 parámetros (mu y sigma)
gl <- k - 2 - 1

gl
## [1] 5
# Valor crítico al 99%
umbral_aceptacion <- qchisq(
  0.999999,
  df = gl
)
umbral_aceptacion
## [1] 35.88819
# Evaluación
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

Evaluación estadística del modelo lognormal de la altura de la columna eruptiva
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Altura de la Columna Eruptiva km 98.14 33.21 35.89

CALCULO DE PROBABILIDAD

Pregunta

Calculo de valores 

# Tamaño de la muestra completa

n_completo <- length(H_pl_limpios)

n_completo
## [1] 898
# ==============================================================================
# PARÁMETRO MEDIA DEL LOGARITMO (μ)
# ==============================================================================


medialog_completo <- mean(
  log(H_pl_limpios)
)


medialog_completo
## [1] 2.419605
# ==============================================================================
# PARÁMETRO DESVIACIÓN ESTÁNDAR DEL LOGARITMO (σ)
# ==============================================================================
sd_log_completo <- sd(
  log(H_pl_limpios)
)

sd_log_completo
## [1] 0.6327913

Calculo de probabilidad 

# Usamos los parámetros completos recalculados
probabilidad_columna <- plnorm(
  30,
  meanlog = medialog_completo,
  sdlog = sd_log_completo
) -
  plnorm(
    20,
    meanlog = medialog_completo,
    sdlog = sd_log_completo
  )

# Probabilidad (%)
probabilidad_columna * 100
## [1] 12.08658

Area de probabilidad

Calculo de probabilidad 2 

# Definición de variables
n_nuevas <- (15 * 5)
p <- probabilidad_columna

# Cálculo del valor esperado
volcanes_esperados <- n_nuevas * p

# Imprimir resultado en consola
print(volcanes_esperados)
## [1] 9.064938

INTERVALOS DE CONFIANZA

Cálculo de parámetros e intervalos 

# Parámetros completos del modelo Log-normal
medialog_completo <- mean(
  log(H_pl_limpios)
)

sd_log_completo <- sd(
  log(H_pl_limpios)
)

# MEDIA DE LA DISTRIBUCIÓN LOGNORMAL EN ESCALA REAL (km)
media_original <- exp(
  medialog_completo + 
    (sd_log_completo^2)/2
)

media_original
## [1] 13.73321
# DESVIACIÓN ESTÁNDAR DE LA DISTRIBUCIÓN LOGNORMAL EN ESCALA REAL (km)
sigma_original <- sqrt(
  (exp(sd_log_completo^2)-1) *
    exp(
      2*medialog_completo + sd_log_completo^2
    )
)

sigma_original
## [1] 9.637333
# TAMAÑO DE MUESTRA
n_completo_Hpl <- length(H_pl_limpios)
n_completo_Hpl
## [1] 898
# ERROR ESTÁNDAR
e <- sigma_original / sqrt(n_completo_Hpl)
e
## [1] 0.321602
# LÍMITE INFERIOR DEL INTERVALO DE CONFIANZA
li <- media_original - 2*e
li
## [1] 13.09
# LÍMITE SUPERIOR DEL INTERVALO DE CONFIANZA
ls <- media_original + 2*e
ls
## [1] 14.37641

Tabla de resumen

Tabla Nro.3: Intervalo de confianza de la media poblacional de la altura de columna eruptiva de los volcanes activos mediante el modelo Log-normal
Límite inferior Media poblacional Límite superior Error estándar
13.09 μ 14.38 0.32

CONCLUSIÓN

La variable altura de la pluma de los volcanes activos se explica mediante un modelo probabilístico Log-normal donde la media aritmética es de 13.73, con una desviación estándar de 9.63.

De esta manera, es posible calcular probabilidades asociadas a la distribución del de la alutra; por ejemplo, determinar que la probabilidad de que un volcán activo presente una columna eruptiva en un rango de km definido.

Mediante el teorema del límite central, se establece que la media aritmética poblacional del período de retorno se encuentra entre 13.09 y 14.38 con un 95% de confianza.

.