Volcán Antisana

1. Introducción

Desempeña un papel crucial en el equilibrio ecológico de la zona, siendo fuente de agua para ríos y lagunas que abastecen a comunidades cercanas, incluyendo Quito

Factores cruciales para el equilibrio ecológico de la región:

Temperatura Precipitaciones Vientos Altitud Latitud

2. Planteamiento de Problema

El análisis del clima en el volcán Antisana es fundamental para comprender su impacto en los ecosistemas locales, optimizar la gestión de recursos hídricos y evaluar los riesgos asociados a la actividad volcánica. Mediante un análisis estadístico, se pueden examinar de manera cuantitativa los datos climáticos, transformándolos en información útil para sugerir y obtener conclusiones que faciliten la toma de decisiones en la conservación y gestión sostenible del área. Consideramos como caso de estudio la Reserva Ecológica Antisana, utilizando datos recopilados de fuentes científicas y plataformas de investigación que abordan el clima y la biodiversidad en esta región.

3.Objetivo General

Aplicar la estadística y el análisis climático al estudio del clima del volcán Antisana para evaluar sus patrones meteorológicos y su influencia en el ecosistema circundante.

3.1.Objetivo Específico

Conocer la situación actual de los datos climáticos del volcán Antisana con el fin de evaluar sus características meteorológicas más importantes y sus medidas estadísticas.

3.2. Objetivo Específico

Emplear modelos de probabilidad para establecer conclusiones sobre el clima del volcán Antisana, basándose en los resultados obtenidos de la recolección de datos climáticos.

3.3. Objetivo Específico

Establecer relaciones entre variables climáticas relevantes del volcán Antisana con el fin de realizar estimaciones significativas sobre su impacto en el ecosistema y la biodiversidad de la zona.

4. Población

Todos los registros climáticos (activos o inactivos) en el volcán Antisana durante el periodo de estudio.

4.1.Simbólico:

P={d/d∈todos los registros del clima del volcán Antisana∧lugar(x)=+∞}

5. Muestra

Todos los registros climáticos (activos o inactivos) en el volcán Antisana durante el periodo de estudio.

5.1. Simbólico:

P={d/d∈Registro Climático∧ubicación(d)=“Volcán Antisana”∧fecha(d)∈[01/01/2012, 30/06/2012]}

6. Individuo

Un subconjunto de los registros climáticos (activos o inactivos) en el volcán Antisana durante el periodo de estudio.

6.1. Simbólico:

S={di/di∈P∧i=1,2,…,m}

Donde:

m<n , siendo n el número total de registros en la población P . di : un registro específico en la muestra, donde i indica el número o la posición de ese registro.

1. ESTADÍSTICA DESCRIPTIVA

1.1. Max Temperature

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)

## 
## Adjuntando el paquete: 'e1071'

## The following object is masked from 'package:ggplot2':
## 
##     element

library(dplyr)

## 
## Adjuntando el paquete: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(readr)
weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

weatherdataANTISANA

Temperatura_Máxima <- weatherdataANTISANA$`Max Temperature`

Tabla de Distribución de Frecuencias

#  Aplicar Regla de Sturges para el número de intervalos (k)
n <- length(Temperatura_Máxima)
k <- ceiling(1 + 3.322 * log10(n)) # Regla de Sturges

#  Calcular cortes y límites
cortes <- seq(min(Temperatura_Máxima), max(Temperatura_Máxima), length.out = k + 1)

#  Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Temperatura_Máxima, breaks = cortes, include.lowest = TRUE)))

#  Calcular cada columna solicitada
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],                                 # Límite Inferior
    Ls = cortes[2:(k+1)],                             # Límite Superior
    MC = (Li + Ls) / 2,                               # Marca de Clase
    ni = Freq,                                        # Frecuencia absoluta
    hi_porc = (ni / n) * 100,                         # Frecuencia relativa %
    Ni_asc = cumsum(ni),                              # Frecuencia absoluta acumulada asc.
    Ni_desc = rev(cumsum(rev(ni))),                   # Frecuencia absoluta acumulada desc.
    Hi_asc_porc = (Ni_asc / n) * 100,                 # Frecuencia relativa acumulada asc. %
    Hi_desc_porc = (Ni_desc / n) * 100                # Frecuencia relativa acumulada desc. %
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

# Ver el resultado
print(tabla_frecuencia)

##        Li     Ls      MC ni   hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1  10.320 11.667 10.9935 21  5.737705     21     366    5.737705   100.000000
## 2  11.667 13.014 12.3405 50 13.661202     71     345   19.398907    94.262295
## 3  13.014 14.361 13.6875 70 19.125683    141     295   38.524590    80.601093
## 4  14.361 15.708 15.0345 51 13.934426    192     225   52.459016    61.475410
## 5  15.708 17.055 16.3815 51 13.934426    243     174   66.393443    47.540984
## 6  17.055 18.402 17.7285 53 14.480874    296     123   80.874317    33.606557
## 7  18.402 19.749 19.0755 34  9.289617    330      70   90.163934    19.125683
## 8  19.749 21.096 20.4225 22  6.010929    352      36   96.174863     9.836066
## 9  21.096 22.443 21.7695  9  2.459016    361      14   98.633880     3.825137
## 10 22.443 23.790 23.1165  5  1.366120    366       5  100.000000     1.366120

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Temperatura_Máxima), aes(x = Temperatura_Máxima)) +
  geom_histogram(breaks = cortes, 
                 fill = "steelblue", 
                 color = "black") + # El color blanco crea la división fina pero las barras están unidas
  labs(title = "Gráfica No. 1: Distribución de Temperatura Máxima",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Temperatura (°C)",
       y = "cantidad") +
  scale_x_continuous(breaks = round(cortes, 1)) + # Muestra los límites exactos en el eje X
  theme_minimal()

Histograma Frecuencia Absoluta —————–

# 4. Crear el Histograma
ggplot(data = data.frame(Temperatura_Máxima), aes(x = Temperatura_Máxima)) +
  geom_histogram(breaks = cortes, 
                 fill = "steelblue", 
                 color = "black") + 
  # Ajustar el límite del eje Y al tamaño muestral (n)
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  # Ajustar eje X a los límites de los intervalos
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 1: Distribución de Temperatura Máxima",
       x = "Temperatura Máxima (°C)",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Temperatura_Máxima), aes(x = Temperatura_Máxima)) +
  # El eje Y mostrará el porcentaje real de cada barra sin forzar el límite a 100
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, 
                 fill = "steelblue", 
                 color = "black") + 
  # Ajustar eje X a los límites de los intervalos para mayor precisión
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3:Distribución Porcentual de Temperatura Máxima",
       x = "Temperatura Máxima (°C)",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

# 4. Crear el Histograma de Porcentajes
ggplot(data = data.frame(Temperatura_Máxima), aes(x = Temperatura_Máxima)) +
  # aes(y = after_stat(count)/sum(after_stat(count)) * 100) calcula el % automáticamente
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, 
                 fill = "steelblue", 
                 color = "black") + 
  # Ajustar el límite del eje Y al 100%
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  # Ajustar eje X a los límites de los intervalos
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4:Distribución Porcentual de Temperatura Máxima )",
       x = "Temperatura Máxima (°C)",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

# 2. Asegurarnos de que los datos estén listos
# (Asumiendo que ya tienes tabla_frecuencia del paso anterior)

# 3. Transformar los datos para ggplot
df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

# 4. Crear la Gráfica de Ojivas
ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + # Usamos linewidth en lugar de size (versiones modernas)
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Temperatura Máxima",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Temperatura (°C)",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

# Usando la función básica de R
boxplot(Temperatura_Máxima, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de CajasNo. 6: Temperatura Máxima",
       
        xlab = "Temperatura (°C)",
        border = "darkblue",
        pch = 16,        # Forma de los puntos atípicos
        col.alpha = 0.5) # Transparencia

Cálculo de Indicadores Estadístocos

# 3. Función para calcular la Moda (ya que R base no tiene una directa)
get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

# 4. Cálculo de Indicadores
minimo  <- min(Temperatura_Máxima, na.rm = TRUE)
maximo  <- max(Temperatura_Máxima, na.rm = TRUE)
media   <- mean(Temperatura_Máxima, na.rm = TRUE)
mediana <- median(Temperatura_Máxima, na.rm = TRUE)
moda    <- get_moda(Temperatura_Máxima)
sd_val  <- sd(Temperatura_Máxima, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Temperatura_Máxima, na.rm = TRUE)
k_val   <- kurtosis(Temperatura_Máxima, na.rm = TRUE)

# 5. Crear Tabla de Resultados
indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

# Mostrar tabla con formato
knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de la Temperatura Máxima - Antisana")

Indicadores Estadísticos de la Temperatura Máxima - Antisana

Indicador	Valor
Mínimo	10.32
Máximo	23.79
Media	15.74
Mediana	15.51
Moda	14.77
SD (Desv. Est.)	2.87
CV (Coef. Var. %)	18.22
As (Asimetría)	0.39
K (Curtosis)	-0.56

Conclusión

Temperatura Máxima en el Antisana presenta una alta estabilidad térmica, con una Media y Mediana muy cercanas que indican un comportamiento predecible. La Curtosis (K) positiva refleja una fuerte concentración de datos alrededor del promedio (distribución leptocúrtica), mientras que el bajo Coeficiente de Variación (CV) corrobora la escasa variabilidad anual. No obstante, la Asimetría (As) y la presencia de valores atípicos en el Diagrama de Cajas revelan picos de calor extremos que, aunque inusuales, son los principales responsables del riesgo de deshielo en los glaciares del volcán.

1.2. Temperatura Mínima

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)
library(dplyr)

weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Definir la variable de interés: Temperatura Mínima
Temperatura_Mínima <- weatherdataANTISANA$`Min Temperature`

Tabla de Distribución de Frecuencias

# Aplicar Regla de Sturges
n <- length(Temperatura_Mínima)
k <- ceiling(1 + 3.322 * log10(n)) 

# Calcular cortes y límites
cortes <- seq(min(Temperatura_Mínima), max(Temperatura_Mínima), length.out = k + 1)

# Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Temperatura_Mínima, breaks = cortes, include.lowest = TRUE)))

# Calcular columnas
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],
    Ls = cortes[2:(k+1)],
    MC = (Li + Ls) / 2,
    ni = Freq,
    hi_porc = (ni / n) * 100,
    Ni_asc = cumsum(ni),
    Ni_desc = rev(cumsum(rev(ni))),
    Hi_asc_porc = (Ni_asc / n) * 100,
    Hi_desc_porc = (Ni_desc / n) * 100
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

print(tabla_frecuencia)

##       Li    Ls    MC ni    hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1   2.65  3.47  3.06  1  0.2732240      1     366    0.273224   100.000000
## 2   3.47  4.29  3.88  3  0.8196721      4     365    1.092896    99.726776
## 3   4.29  5.11  4.70  6  1.6393443     10     362    2.732240    98.907104
## 4   5.11  5.93  5.52 13  3.5519126     23     356    6.284153    97.267760
## 5   5.93  6.75  6.34 33  9.0163934     56     343   15.300546    93.715847
## 6   6.75  7.57  7.16 73 19.9453552    129     310   35.245902    84.699454
## 7   7.57  8.39  7.98 95 25.9562842    224     237   61.202186    64.754098
## 8   8.39  9.21  8.80 59 16.1202186    283     142   77.322404    38.797814
## 9   9.21 10.03  9.62 56 15.3005464    339      83   92.622951    22.677596
## 10 10.03 10.85 10.44 27  7.3770492    366      27  100.000000     7.377049

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Temperatura_Mínima), aes(x = Temperatura_Mínima)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  labs(title = "Gráfica No. 1: Distribución de Temperatura Mínima",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Temperatura (°C)",
       y = "Cantidad") +
  scale_x_continuous(breaks = round(cortes, 1)) +
  theme_minimal()

Histograma Frecuencia Absoluta —————–

ggplot(data = data.frame(Temperatura_Mínima), aes(x = Temperatura_Mínima)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 2: Distribución de Temperatura Mínima (Escala N)",
       x = "Temperatura Mínima (°C)",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Temperatura_Mínima), aes(x = Temperatura_Mínima)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3: Distribución Porcentual de Temperatura Mínima",
       x = "Temperatura Mínima (°C)",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

ggplot(data = data.frame(Temperatura_Mínima), aes(x = Temperatura_Mínima)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4: Distribución Porcentual de Temperatura Mínima (Escala 100%)",
       x = "Temperatura Mínima (°C)",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Temperatura Mínima",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Temperatura (°C)",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

boxplot(Temperatura_Mínima, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de Cajas No. 6: Temperatura Mínima",
        xlab = "Temperatura (°C)",
        border = "darkblue",
        pch = 16)

Cálculo de Indicadores Estadístocos

get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

minimo  <- min(Temperatura_Mínima, na.rm = TRUE)
maximo  <- max(Temperatura_Mínima, na.rm = TRUE)
media   <- mean(Temperatura_Mínima, na.rm = TRUE)
mediana <- median(Temperatura_Mínima, na.rm = TRUE)
moda    <- get_moda(Temperatura_Mínima)
sd_val  <- sd(Temperatura_Mínima, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Temperatura_Mínima, na.rm = TRUE)
k_val   <- kurtosis(Temperatura_Mínima, na.rm = TRUE)

indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de la Temperatura Mínima - Antisana")

Indicadores Estadísticos de la Temperatura Mínima - Antisana

Indicador	Valor
Mínimo	2.65
Máximo	10.85
Media	8.05
Mediana	8.00
Moda	7.89
SD (Desv. Est.)	1.37
CV (Coef. Var. %)	17.08
As (Asimetría)	-0.43
K (Curtosis)	0.57

Conclusión

El análisis de la Temperatura Mínima en el Antisana muestra una distribución con una Media y Mediana que reflejan el rigor del clima de alta montaña nocturno. La Curtosis (K) y el CV indican una variabilidad moderada, mientras que la Asimetría (As) negativa sugeriría una tendencia hacia eventos de heladas más intensas. El Diagrama de Cajas permite visualizar que, aunque hay estabilidad, existen valores atípicos inferiores que representan los descensos térmicos más extremos registrados en la estación.

1.3. Precipitation

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)
library(dplyr)

weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Definir la variable de interés: Precipitación
Precipitacion <- weatherdataANTISANA$Precipitation

Tabla de Distribución de Frecuencias

# Aplicar Regla de Sturges
n <- length(Precipitacion)
k <- ceiling(1 + 3.322 * log10(n)) 

# Calcular cortes y límites
cortes <- seq(min(Precipitacion), max(Precipitacion), length.out = k + 1)

# Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Precipitacion, breaks = cortes, include.lowest = TRUE)))

# Calcular columnas
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],
    Ls = cortes[2:(k+1)],
    MC = (Li + Ls) / 2,
    ni = Freq,
    hi_porc = (ni / n) * 100,
    Ni_asc = cumsum(ni),
    Ni_desc = rev(cumsum(rev(ni))),
    Hi_asc_porc = (Ni_asc / n) * 100,
    Hi_desc_porc = (Ni_desc / n) * 100
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

print(tabla_frecuencia)

##        Li     Ls      MC  ni    hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1   0.010  9.481  4.7455 146 39.8907104    146     366    39.89071  100.0000000
## 2   9.481 18.952 14.2165  91 24.8633880    237     220    64.75410   60.1092896
## 3  18.952 28.423 23.6875  55 15.0273224    292     129    79.78142   35.2459016
## 4  28.423 37.894 33.1585  29  7.9234973    321      74    87.70492   20.2185792
## 5  37.894 47.365 42.6295  24  6.5573770    345      45    94.26230   12.2950820
## 6  47.365 56.836 52.1005  10  2.7322404    355      21    96.99454    5.7377049
## 7  56.836 66.307 61.5715   9  2.4590164    364      11    99.45355    3.0054645
## 8  66.307 75.778 71.0425   0  0.0000000    364       2    99.45355    0.5464481
## 9  75.778 85.249 80.5135   0  0.0000000    364       2    99.45355    0.5464481
## 10 85.249 94.720 89.9845   2  0.5464481    366       2   100.00000    0.5464481

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Precipitacion), aes(x = Precipitacion)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  labs(title = "Gráfica No. 1: Distribución de Precipitación",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Precipitación (mm)",
       y = "Cantidad de días") +
  scale_x_continuous(breaks = round(cortes, 1)) +
  theme_minimal()

Histograma Frecuencia Absoluta —————–

ggplot(data = data.frame(Precipitacion), aes(x = Precipitacion)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 2: Distribución de Precipitación (Escala N)",
       x = "Precipitación (mm)",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Precipitacion), aes(x = Precipitacion)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3: Distribución Porcentual de Precipitación",
       x = "Precipitación (mm)",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

ggplot(data = data.frame(Precipitacion), aes(x = Precipitacion)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4: Distribución Porcentual de Precipitación (Escala 100%)",
       x = "Precipitación (mm)",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Precipitación",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Precipitación (mm)",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

boxplot(Precipitacion, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de Cajas No. 6: Precipitación",
        xlab = "Precipitación (mm)",
        border = "darkblue",
        pch = 16)

Cálculo de Indicadores Estadístocos

get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

minimo  <- min(Precipitacion, na.rm = TRUE)
maximo  <- max(Precipitacion, na.rm = TRUE)
media   <- mean(Precipitacion, na.rm = TRUE)
mediana <- median(Precipitacion, na.rm = TRUE)
moda    <- get_moda(Precipitacion)
sd_val  <- sd(Precipitacion, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Precipitacion, na.rm = TRUE)
k_val   <- kurtosis(Precipitacion, na.rm = TRUE)

indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de Precipitación - Antisana")

Indicadores Estadísticos de Precipitación - Antisana

Indicador	Valor
Mínimo	0.01
Máximo	94.72
Media	17.10
Mediana	12.94
Moda	0.01
SD (Desv. Est.)	16.12
CV (Coef. Var. %)	94.21
As (Asimetría)	1.30
K (Curtosis)	1.95

Conclusión

La Precipitación en el Antisana presenta una distribución típica de variables hidrológicas, caracterizada por una fuerte Asimetría Positiva, donde la mayoría de los días presentan valores bajos o nulos (moda en cero), pero con eventos extremos de lluvia intensa reflejados en los valores Máximo y los outliers del Diagrama de Cajas. La Curtosis elevada y un CV significativamente alto indican que el régimen de lluvias es muy variable y está concentrado en eventos específicos, lo cual es fundamental para el aporte hídrico de los glaciares.

1.4. Velocidad del Viento

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)
library(dplyr)

weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Definir la variable de interés: Viento
Viento <- weatherdataANTISANA$Wind

Tabla de Distribución de Frecuencias

# Aplicar Regla de Sturges
n <- length(Viento)
k <- ceiling(1 + 3.322 * log10(n)) 

# Calcular cortes y límites
cortes <- seq(min(Viento), max(Viento), length.out = k + 1)

# Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Viento, breaks = cortes, include.lowest = TRUE)))

# Calcular columnas
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],
    Ls = cortes[2:(k+1)],
    MC = (Li + Ls) / 2,
    ni = Freq,
    hi_porc = (ni / n) * 100,
    Ni_asc = cumsum(ni),
    Ni_desc = rev(cumsum(rev(ni))),
    Hi_asc_porc = (Ni_asc / n) * 100,
    Hi_desc_porc = (Ni_desc / n) * 100
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

print(tabla_frecuencia)

##      Li   Ls   MC ni    hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1  0.59 0.83 0.71  2  0.5464481      2     366   0.5464481  100.0000000
## 2  0.83 1.07 0.95 19  5.1912568     21     364   5.7377049   99.4535519
## 3  1.07 1.31 1.19 43 11.7486339     64     345  17.4863388   94.2622951
## 4  1.31 1.55 1.43 66 18.0327869    130     302  35.5191257   82.5136612
## 5  1.55 1.79 1.67 67 18.3060109    197     236  53.8251366   64.4808743
## 6  1.79 2.03 1.91 57 15.5737705    254     169  69.3989071   46.1748634
## 7  2.03 2.27 2.15 51 13.9344262    305     112  83.3333333   30.6010929
## 8  2.27 2.51 2.39 41 11.2021858    346      61  94.5355191   16.6666667
## 9  2.51 2.75 2.63 18  4.9180328    364      20  99.4535519    5.4644809
## 10 2.75 2.99 2.87  2  0.5464481    366       2 100.0000000    0.5464481

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Viento), aes(x = Viento)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  labs(title = "Gráfica No. 1: Distribución de la Velocidad del Viento",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Velocidad del Viento (km/h)",
       y = "Cantidad de días") +
  scale_x_continuous(breaks = round(cortes, 1)) +
  theme_minimal()

Histograma Frecuencia Absoluta —————–

ggplot(data = data.frame(Viento), aes(x = Viento)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 2: Distribución de Viento (Escala N)",
       x = "Velocidad del Viento (km/h)",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Viento), aes(x = Viento)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3: Distribución Porcentual del Viento",
       x = "Velocidad del Viento (km/h)",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

ggplot(data = data.frame(Viento), aes(x = Viento)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4: Distribución Porcentual del Viento (Escala 100%)",
       x = "Velocidad del Viento (km/h)",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Velocidad del Viento",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Velocidad (km/h)",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

boxplot(Viento, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de Cajas No. 6: Velocidad del Viento",
        xlab = "Velocidad (km/h)",
        border = "darkblue",
        pch = 16)

Cálculo de Indicadores Estadístocos

get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

minimo  <- min(Viento, na.rm = TRUE)
maximo  <- max(Viento, na.rm = TRUE)
media   <- mean(Viento, na.rm = TRUE)
mediana <- median(Viento, na.rm = TRUE)
moda    <- get_moda(Viento)
sd_val  <- sd(Viento, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Viento, na.rm = TRUE)
k_val   <- kurtosis(Viento, na.rm = TRUE)

indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de la Velocidad del Viento - Antisana")

Indicadores Estadísticos de la Velocidad del Viento - Antisana

Indicador	Valor
Mínimo	0.59
Máximo	2.99
Media	1.77
Mediana	1.75
Moda	1.49
SD (Desv. Est.)	0.46
CV (Coef. Var. %)	25.78
As (Asimetría)	0.11
K (Curtosis)	-0.72

Conclusión

La velocidad del Viento en el Antisana muestra una distribución característica de zonas de alta montaña, donde la Media refleja la persistencia de corrientes atmosféricas en el páramo. La Asimetría suele ser positiva, indicando que aunque predominan velocidades moderadas, existen ráfagas intensas representadas por el valor Máximo. El CV elevado sugiere que el viento es una variable con alta fluctuación diaria, mientras que los valores atípicos en el Diagrama de Cajas señalan eventos de vientos fuertes que impactan directamente en la sensación térmica y la evaporación del glaciar.

1.5. Humedad Relativa

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)
library(dplyr)

weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Definir la variable de interés: Humedad Relativa
Humedad_Relativa <- weatherdataANTISANA$`Relative Humidity`

Tabla de Distribución de Frecuencias

# Aplicar Regla de Sturges
n <- length(Humedad_Relativa)
k <- ceiling(1 + 3.322 * log10(n)) 

# Calcular cortes y límites
cortes <- seq(min(Humedad_Relativa), max(Humedad_Relativa), length.out = k + 1)

# Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Humedad_Relativa, breaks = cortes, include.lowest = TRUE)))

# Calcular columnas
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],
    Ls = cortes[2:(k+1)],
    MC = (Li + Ls) / 2,
    ni = Freq,
    hi_porc = (ni / n) * 100,
    Ni_asc = cumsum(ni),
    Ni_desc = rev(cumsum(rev(ni))),
    Hi_asc_porc = (Ni_asc / n) * 100,
    Hi_desc_porc = (Ni_desc / n) * 100
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

print(tabla_frecuencia)

##       Li    Ls     MC  ni    hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1  0.560 0.603 0.5815   2  0.5464481      2     366   0.5464481    100.00000
## 2  0.603 0.646 0.6245  13  3.5519126     15     364   4.0983607     99.45355
## 3  0.646 0.689 0.6675  14  3.8251366     29     351   7.9234973     95.90164
## 4  0.689 0.732 0.7105  20  5.4644809     49     337  13.3879781     92.07650
## 5  0.732 0.775 0.7535  14  3.8251366     63     317  17.2131148     86.61202
## 6  0.775 0.818 0.7965  17  4.6448087     80     303  21.8579235     82.78689
## 7  0.818 0.861 0.8395  25  6.8306011    105     286  28.6885246     78.14208
## 8  0.861 0.904 0.8825  35  9.5628415    140     261  38.2513661     71.31148
## 9  0.904 0.947 0.9255  47 12.8415301    187     226  51.0928962     61.74863
## 10 0.947 0.990 0.9685 179 48.9071038    366     179 100.0000000     48.90710

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Humedad_Relativa), aes(x = Humedad_Relativa)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  labs(title = "Gráfica No. 1: Distribución de Humedad Relativa",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Humedad Relativa (%)",
       y = "Cantidad de días") +
  scale_x_continuous(breaks = round(cortes, 1)) +
  theme_minimal()

Histograma Frecuencia Absoluta —————–

ggplot(data = data.frame(Humedad_Relativa), aes(x = Humedad_Relativa)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 2: Distribución de Humedad Relativa (Escala N)",
       x = "Humedad Relativa (%)",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Humedad_Relativa), aes(x = Humedad_Relativa)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3: Distribución Porcentual de Humedad Relativa",
       x = "Humedad Relativa (%)",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

ggplot(data = data.frame(Humedad_Relativa), aes(x = Humedad_Relativa)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4: Distribución Porcentual de Humedad Relativa (Escala 100%)",
       x = "Humedad Relativa (%)",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Humedad Relativa",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Humedad (%)",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

boxplot(Humedad_Relativa, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de Cajas No. 6: Humedad Relativa",
        xlab = "Humedad (%)",
        border = "darkblue",
        pch = 16)

Cálculo de Indicadores Estadístocos

get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

minimo  <- min(Humedad_Relativa, na.rm = TRUE)
maximo  <- max(Humedad_Relativa, na.rm = TRUE)
media   <- mean(Humedad_Relativa, na.rm = TRUE)
mediana <- median(Humedad_Relativa, na.rm = TRUE)
moda    <- get_moda(Humedad_Relativa)
sd_val  <- sd(Humedad_Relativa, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Humedad_Relativa, na.rm = TRUE)
k_val   <- kurtosis(Humedad_Relativa, na.rm = TRUE)

indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de Humedad Relativa - Antisana")

Indicadores Estadísticos de Humedad Relativa - Antisana

Indicador	Valor
Mínimo	0.56
Máximo	0.99
Media	0.90
Mediana	0.94
Moda	0.98
SD (Desv. Est.)	0.11
CV (Coef. Var. %)	12.14
As (Asimetría)	-1.18
K (Curtosis)	0.20

Conclusión

El análisis de la Humedad Relativa en el Antisana revela un ambiente con alta saturación hídrica, típica de los ecosistemas de páramo y glaciares altoandinos. La Media y Mediana elevadas confirman la persistencia de condiciones húmedas, mientras que la Asimetría Negativa (común en esta variable al estar topada al 100%) indica que la mayoría de los registros se concentran en valores altos. El CV moderado sugiere una estabilidad relativa en la humedad, aunque los valores atípicos inferiores en el Diagrama de Cajas señalan eventos de sequedad atmosférica extrema que pueden incrementar la sublimación del hielo glaciar.

1.6. Solar

Datos

library(readr)
library(tidyr)
library(ggplot2)
library(e1071)
library(dplyr)

weatherdataANTISANA <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Definir la variable de interés: Radiación Solar
Radiacion_Solar <- weatherdataANTISANA$Solar

Tabla de Distribución de Frecuencias

# Aplicar Regla de Sturges
n <- length(Radiacion_Solar)
k <- ceiling(1 + 3.322 * log10(n)) 

# Calcular cortes y límites
cortes <- seq(min(Radiacion_Solar), max(Radiacion_Solar), length.out = k + 1)

# Crear la tabla de frecuencias
tabla_frec <- as.data.frame(table(cut(Radiacion_Solar, breaks = cortes, include.lowest = TRUE)))

# Calcular columnas
tabla_frecuencia <- tabla_frec %>%
  mutate(
    Li = cortes[1:k],
    Ls = cortes[2:(k+1)],
    MC = (Li + Ls) / 2,
    ni = Freq,
    hi_porc = (ni / n) * 100,
    Ni_asc = cumsum(ni),
    Ni_desc = rev(cumsum(rev(ni))),
    Hi_asc_porc = (Ni_asc / n) * 100,
    Hi_desc_porc = (Ni_desc / n) * 100
  ) %>%
  select(Li, Ls, MC, ni, hi_porc, Ni_asc, Ni_desc, Hi_asc_porc, Hi_desc_porc)

print(tabla_frecuencia)

##        Li     Ls      MC ni   hi_porc Ni_asc Ni_desc Hi_asc_porc Hi_desc_porc
## 1   1.260  4.161  2.7105 34  9.289617     34     366    9.289617   100.000000
## 2   4.161  7.062  5.6115 51 13.934426     85     332   23.224044    90.710383
## 3   7.062  9.963  8.5125 50 13.661202    135     281   36.885246    76.775956
## 4   9.963 12.864 11.4135 49 13.387978    184     231   50.273224    63.114754
## 5  12.864 15.765 14.3145 37 10.109290    221     182   60.382514    49.726776
## 6  15.765 18.666 17.2155 21  5.737705    242     145   66.120219    39.617486
## 7  18.666 21.567 20.1165 25  6.830601    267     124   72.950820    33.879781
## 8  21.567 24.468 23.0175 29  7.923497    296      99   80.874317    27.049180
## 9  24.468 27.369 25.9185 45 12.295082    341      70   93.169399    19.125683
## 10 27.369 30.270 28.8195 25  6.830601    366      25  100.000000     6.830601

Histograma Frecuencia Absoluta

ggplot(data = data.frame(Radiacion_Solar), aes(x = Radiacion_Solar)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  labs(title = "Gráfica No. 1: Distribución de Radiación Solar",
       subtitle = "Volcán Antisana (Regla de Sturges)",
       x = "Radiación Solar (W/m² o Ly/day)",
       y = "Cantidad de días") +
  scale_x_continuous(breaks = round(cortes, 1)) +
  theme_minimal()

Histograma Frecuencia Absoluta —————–

ggplot(data = data.frame(Radiacion_Solar), aes(x = Radiacion_Solar)) +
  geom_histogram(breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, n), breaks = seq(0, n, by = 50)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 2: Distribución de Radiación Solar (Escala N)",
       x = "Radiación Solar",
       y = "Cantidad") +
  theme_minimal()

Histograma Porcentaje

ggplot(data = data.frame(Radiacion_Solar), aes(x = Radiacion_Solar)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 3: Distribución Porcentual de Radiación Solar",
       x = "Radiación Solar",
       y = "Porcentaje (%)") +
  theme_minimal()

Histograma Porcentaje ——————

ggplot(data = data.frame(Radiacion_Solar), aes(x = Radiacion_Solar)) +
  geom_histogram(aes(y = after_stat(count)/sum(after_stat(count)) * 100),
                 breaks = cortes, fill = "steelblue", color = "black") + 
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  scale_x_continuous(breaks = round(cortes, 1)) +
  labs(title = "Gráfica No. 4: Distribución Porcentual de Radiación Solar (Escala 100%)",
       x = "Radiación Solar",
       y = "Porcentaje (%)") +
  theme_minimal()

Grafica de Ojivas

df_ojiva <- tabla_frecuencia %>%
  select(MC, Hi_asc_porc, Hi_desc_porc) %>%
  pivot_longer(cols = c(Hi_asc_porc, Hi_desc_porc), 
               names_to = "Tipo_Ojiva", 
               values_to = "Porcentaje")

ggplot(df_ojiva, aes(x = MC, y = Porcentaje, color = Tipo_Ojiva, group = Tipo_Ojiva)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  scale_color_manual(values = c("Hi_asc_porc" = "#2c7bb6", "Hi_desc_porc" = "#d7191c"),
                     labels = c("Ascendente (Menor que)", "Descendente (Mayor que)")) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, by = 10)) +
  labs(title = "Gráfica de Ojivas No. 5: Radiación Solar",
       subtitle = "Estación Antisana (4,048 msnm) - Año 2012",
       x = "Radiación Solar",
       y = "Porcentaje Acumulado (%)",
       color = "Interpretación") +
  theme_minimal() +
  theme(legend.position = "bottom")

Diagrama de Cajas

boxplot(Radiacion_Solar, 
        horizontal = TRUE, 
        col = "skyblue", 
        main = "Diagrama de Cajas No. 6: Radiación Solar",
        xlab = "Radiación Solar",
        border = "darkblue",
        pch = 16)

Cálculo de Indicadores Estadístocos

get_moda <- function(x) {
  uniqv <- unique(x)
  uniqv[which.max(tabulate(match(x, uniqv)))]
}

minimo  <- min(Radiacion_Solar, na.rm = TRUE)
maximo  <- max(Radiacion_Solar, na.rm = TRUE)
media   <- mean(Radiacion_Solar, na.rm = TRUE)
mediana <- median(Radiacion_Solar, na.rm = TRUE)
moda    <- get_moda(Radiacion_Solar)
sd_val  <- sd(Radiacion_Solar, na.rm = TRUE)
cv_val  <- (sd_val / media) * 100
as_val  <- skewness(Radiacion_Solar, na.rm = TRUE)
k_val   <- kurtosis(Radiacion_Solar, na.rm = TRUE)

indicadores <- data.frame(
  Indicador = c("Mínimo", "Máximo", "Media", "Mediana", "Moda", 
                "SD (Desv. Est.)", "CV (Coef. Var. %)", "As (Asimetría)", "K (Curtosis)"),
  Valor = c(minimo, maximo, media, mediana, moda, 
            sd_val, cv_val, as_val, k_val)
)

knitr::kable(indicadores, digits = 2, caption = "Indicadores Estadísticos de Radiación Solar - Antisana")

Indicadores Estadísticos de Radiación Solar - Antisana

Indicador	Valor
Mínimo	1.26
Máximo	30.27
Media	14.44
Mediana	12.66
Moda	7.19
SD (Desv. Est.)	8.33
CV (Coef. Var. %)	57.69
As (Asimetría)	0.30
K (Curtosis)	-1.24

Conclusión

El análisis de la Radiación Solar en el volcán Antisana muestra una distribución con una variabilidad importante, típica de una zona donde la cobertura nubosa cambia drásticamente. La Media y la Mediana indican la energía promedio disponible, mientras que una Asimetría marcada y los valores extremos en el Diagrama de Cajas resaltan los días de “cielo despejado” con alta incidencia de radiación de onda corta. Estos picos de radiación, combinados con la baja presión atmosférica de la altitud, son el factor energético primario que impulsa la ablación glaciar en el sitio de estudio.

1. ESTADÍSTICA INFERENCIAL

1. PRECIPITACIÓN

Datos

1 Identificación y Justificación

La precipitación se define como una variable cuantitativa continua, ya que se expresa mediante valores numéricos reales, generalmente medidos en milímetros, que pueden adoptar cualquier valor decimal dentro de un intervalo determinado. Esta característica de continuidad permite describir con precisión la magnitud de los eventos de lluvia, reflejando de manera adecuada la variabilidad y la intensidad de la precipitación registrada en cada evento.

2 Datos

Importamos el archivo “VOLCAN A.csv” desde una ruta local y lo almacena en el objeto datos, usando espacios o tabulaciones como separador.

datos <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(datos)

3 Extraer Variable

Extraemos la variab precipitación, omitimos las celdas en blanco o valores iguales a cero y verificamos el tamaño muestral

precipitacion <- na.omit(datos$Precipitation)
precipitacion <- precipitacion [precipitacion > 0]
n_total <- length(precipitacion)

4 Tabla de Frecuencia

En la tabla de distribución de frecuencias de la variable Precipitación, se definió un número de 5 clases de forma preestablecida para optimizar la visualización. El ancho de cada intervalo se calculó dividiendo el rango total de los datos entre l número fijo de clases.

# 1. Preparación de la variable
datos_analizar <- datos$Precipitation

# 2. Cálculos base
xmin <- min(datos_analizar, na.rm = TRUE)
xmax <- max(datos_analizar, na.rm = TRUE)
R <- xmax - xmin

# Definición manual de clases
K <- 5
A <- R / K

# 3. Definición de Límites y Marcas de Clase (MC)
cortes <- seq(xmin, xmax, length.out = K + 1)
Li <- head(cortes, -1)
Ls <- tail(cortes, -1)
MC <- (Li + Ls) / 2

# 4. Cálculo de Frecuencias
ni <- as.vector(table(cut(datos_analizar, breaks = cortes, include.lowest = TRUE, right = FALSE)))

# 5. Frecuencias relativas y acumuladas
hi <- ni / sum(ni) * 100
Ni_asc <- cumsum(ni)
Ni_desc <- rev(cumsum(rev(ni)))
Hi_asc <- cumsum(hi)
Hi_desc <- rev(cumsum(rev(hi)))

# 6. Creación del Data Frame final (TDF)
TDF <- data.frame(
  Li = round(Li, 2), 
  Ls = round(Ls, 2), 
  MC = round(MC, 2), 
  ni = ni, 
  hi_porc = round(hi, 2), 
  Ni_asc = Ni_asc, 
  Ni_desc = Ni_desc, 
  Hi_asc_porc = round(Hi_asc, 2), 
  Hi_desc_porc = round(Hi_desc, 2)
)

# 7. Mostrar Tabla
kable(TDF, 
      caption = "Tabla No. 1: Distribución de Frecuencias de Precipitación",
      col.names = c("Lím. Inf.", "Lím. Sup.", "Marca Clase", "ni", "hi (%)", 
                    "Ni Asc.", "Ni Desc.", "Hi Asc. (%)", "Hi Desc. (%)"),
      digits = 2)

Tabla No. 1: Distribución de Frecuencias de Precipitación

Lím. Inf.	Lím. Sup.	Marca Clase	ni	hi (%)	Ni Asc.	Ni Desc.	Hi Asc. (%)	Hi Desc. (%)
0.01	18.95	9.48	237	64.75	237	366	64.75	100.00
18.95	37.89	28.42	84	22.95	321	129	87.70	35.25
37.89	56.84	47.36	34	9.29	355	45	96.99	12.30
56.84	75.78	66.31	9	2.46	364	11	99.45	3.01
75.78	94.72	85.25	2	0.55	366	2	100.00	0.55

5 Histograma

Una vez generada la Tabla de Distribución de Frecuencias, procedemos a visualizar los datos. Esta gráfica es fundamental para identificar la asimetría de la variable continua de precipitación y justificar el uso de modelos probabilísticos posteriores.

# Generación de la Gráfica No. 1
ggplot(TDF, aes(x = MC, y = hi_porc)) +
  
  # Geometría de barras: El ancho 'width' debe ser igual a la amplitud 'A'
  geom_bar(stat = "identity", 
           fill = "steelblue", 
           color = "black", 
           alpha = 0.8,
           width = A) + 
  
  # Configuración del eje X: Mostramos los límites inferiores y el último superior
  scale_x_continuous(breaks = c(Li, Ls[K]), 
                     labels = function(x) round(x, 2)) +
  
  # Configuración del eje Y: Formato de porcentaje y margen superior del 10%
  scale_y_continuous(labels = function(x) paste0(x, "%"),
                     limits = c(0, max(TDF$hi_porc) * 1.1),
                     expand = c(0, 0)) +
  
  # Etiquetas informativas
  labs(
    title = "Gráfica No. 1: Histograma de Distribución Porcentual de Precipitación",
    x = "Precipitación (mm)",
    y = "Porcentaje (%)",
    caption = "Fuente: Datos VOLCAN A"
  ) +
  
  # Estilo visual limpio
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

6 Conjetura de modelo Exponencial

Se selecciono el modelo exponencial para esta vatiable porque visualmente las barras muestran un decaimiento constante este comportamiento “en escalera descendente” es la firma visual característica de una función exponencial .

# 1. Estimación del parámetro Lambda
media_precip <- mean(datos_analizar, na.rm = TRUE)
lambda_exp <- 1 / media_precip

# 2. Generación de los puntos para la curva teórica
x_curva <- seq(xmin, xmax, length.out = 300)

# Escalado de la densidad: dexp() * Amplitud * 100
y_curva <- dexp(x_curva, rate = lambda_exp) * A * 100
df_conjetura <- data.frame(x = x_curva, y = y_curva)

# 3. Construcción de la Gráfica No. 2
ggplot(TDF, aes(x = MC, y = hi_porc)) +
  
  # Histograma (usando la amplitud A calculada en la TDF)
  geom_bar(stat = "identity", 
           fill = "steelblue", 
           color = "black", 
           alpha = 0.8,
           width = A) + 
  
  # Superposición de la CURVA EXPONENCIAL - SE USÓ linewidth PARA EVITAR WARNING
  geom_line(data = df_conjetura, aes(x = x, y = y), 
            color = "red", linewidth = 1.2) +
  
  # Configuración del eje X: Límites de los intervalos
  scale_x_continuous(breaks = c(Li, Ls[K]), 
                     labels = function(x) round(x, 2)) +
  
  # Configuración del eje Y: Formato de porcentaje
  scale_y_continuous(labels = function(x) paste0(x, "%"),
                     limits = c(0, max(c(TDF$hi_porc, df_conjetura$y)) * 1.1),
                     expand = c(0, 0)) +
  
  # Etiquetas y títulos
  labs(
    title = "Gráfica No. 2: Modelo de probabilidad exponencial para precipitación",
    subtitle = paste("Modelo: f(x) = λe^(-λx) | Lambda (λ) =", round(lambda_exp, 5)),
    x = "Precipitación (mm)",
    y = "Porcentaje (%)",
    caption = "Fuente: Datos obtenidos de VOLCAN A.csv"
  ) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

7 Test de PERSON

Para verificar la validez de la conjetura, comparamos las frecuencias observadas obtenidas de los datos con las frecuencias esperadas que predice el modelo exponencial. La correlación de Pearson nos permite cuantificar qué tan bien se ajusta la curva teórica a la distribución real de la precipitación.

# 1. Definir Frecuencia Observada
Fo <- ni 

# 2. Calcular Probabilidades Teóricas (P)
P <- numeric(K)
for (i in 1:K) {
  P[i] <- pexp(cortes[i+1], rate = lambda_exp) - 
          pexp(cortes[i], rate = lambda_exp)
}

# 3. Calcular Frecuencia Esperada (Fe)
# Se multiplica la probabilidad teórica por el tamaño total de la muestra
n_total <- length(datos_analizar)
Fe <- P * n_total

# 4. Cálculo de Correlación de Pearson
Correlacion_Exp <- cor(Fo, Fe) * 100

# 5. Resultados de la Validación
cat("--- CONCLUSIÓN DE LA CONJETURA EXPONENCIAL (PRECIPITACIÓN) ---\n")

## --- CONCLUSIÓN DE LA CONJETURA EXPONENCIAL (PRECIPITACIÓN) ---

cat("Tamaño de la muestra total:", n_total, "registros\n")

## Tamaño de la muestra total: 366 registros

cat("Correlación del modelo:", round(Correlacion_Exp, 2), "%\n")

## Correlación del modelo: 99.92 %

8 Test de CHI-CUADRADO

Para dar rigor estadístico a la conjetura, aplicamos la prueba de Chi-cuadrado. Esta prueba evalúa si las discrepancias entre los porcentajes observados en la muestra y los porcentajes teóricos del modelo exponencial son lo suficientemente pequeñas como para ser atribuidas al azar, o si por el contrario, el modelo debe ser rechazado.

# 1. Grados de libertad 
grados_libertad <- (K - 3) 

# 2. Nivel de significancia 
nivel_significancia <- 0.5

# 3. Preparación de Frecuencias Porcentuales
# Fo_porc: Porcentajes observados (hi de la TDF)
# Fe_porc: Porcentajes teóricos (basados en la integral de la curva)
Fo_porc <- hi
Fe_porc <- P * 100

# 4. Cálculo del estadístico Chi-Cuadrado (X2)
# Sumatoria de (Obs - Esp)^2 / Esp
x2 <- sum((Fe_porc - Fo_porc)^2 / Fe_porc)

# 5. Determinación del Umbral (Valor Crítico)
# Punto de corte en la distribución Chi-cuadrado
umbral_aceptacion <- qchisq(1 - nivel_significancia, grados_libertad)

# 6. Decisión Estadística
modelo_aceptado <- x2 < umbral_aceptacion

# --- Impresión de Resultados ---
cat("--- PRUEBA DE CHI-CUADRADO PARA PRECIPITACIÓN ---\n")

## --- PRUEBA DE CHI-CUADRADO PARA PRECIPITACIÓN ---

cat("Grados de Libertad:", grados_libertad, "\n")

## Grados de Libertad: 2

cat("Estadístico Chi-Cuadrado (X2) calculado:", round(x2, 4), "\n")

## Estadístico Chi-Cuadrado (X2) calculado: 0.7206

cat("Umbral de aceptación (Valor Crítico):", round(umbral_aceptacion, 4), "\n")

## Umbral de aceptación (Valor Crítico): 1.3863

cat("¿El modelo es aceptado estadísticamente?:", ifelse(modelo_aceptado, "SÍ", "NO"), "\n")

## ¿El modelo es aceptado estadísticamente?: SÍ

9 Tabla de resumen

Se presenta la siguiente tabla comparativa que consolida los resultados de la Correlación de Pearson y la prueba de Chi-cuadrado. Este resumen permite confirmar si el modelo propuesto es estadísticamente válido para describir el comportamiento de la precipitación en el Volcán A.

# 1. Definición de la etiqueta del modelo analizado
Modelos <- c("Modelo Exponencial Puro")

# 2. Creación del Dataframe con los resultados consolidados
tabla_resumen_precip <- data.frame(
  Modelo = Modelos,
  Pearson = c(round(Correlacion_Exp, 2)),
  Chi_Sq = c(round(x2, 2)),
  Umbral = c(round(umbral_aceptacion, 2)),
  Resultado = c(ifelse(modelo_aceptado, "ACEPTADO", "RECHAZADO"))
)

# 3. Renombrar columnas para una presentación profesional
colnames(tabla_resumen_precip) <- c("Modelo de Probabilidad", 
                                   "Correlación Pearson (%)", 
                                   "Chi Cuadrado (X2)", 
                                   "Umbral de Aceptación",
                                   "Decisión Final")

# 4. Generación de la tabla con formato kable
library(knitr)
kable(tabla_resumen_precip, 
      format = "markdown", 
      align = "lcccc",
      caption = "Tabla Nro. 2: Resumen de test de bondad de para el modelo exponencial para  Precipitación")

Tabla Nro. 2: Resumen de test de bondad de para el modelo exponencial para Precipitación

Modelo de Probabilidad	Correlación Pearson (%)	Chi Cuadrado (X2)	Umbral de Aceptación	Decisión Final
Modelo Exponencial Puro	99.92	0.72	1.39	ACEPTADO

10 Cálculo de Probabilidades

Una vez validado el modelo exponencial, procedemos a resolver interrogantes clave sobre el comportamiento futuro de la precipitación en el Volcán A, utilizando las funciones de distribución de probabilidad.

¿Cuál es la cantidad de precipitación (en mm) estimada, tal que el 90% de los eventos futuros de lluvia presentarían valores iguales o menores a dicha cantidad?

# 1. Cálculo del cuantil 0.90 usando el parámetro lambda_exp
precip_90 <- qexp(0.90, rate = lambda_exp)

# 2. Resultado
cat("Resultado: ", round(precip_90, 2), "mm")

## Resultado:  39.39 mm

¿Cuál es la probabilidad de que, en un evento futuro de precipitación, el acumulado se encuentre entre 20 mm y 60 mm, de acuerdo con el modelo exponencial estimado?

# 1. Cálculo de la probabilidad (P(20 < X < 60))
probabilidad_rango <- pexp(60, rate = lambda_exp) - pexp(20, rate = lambda_exp)

# 2. Conversión a porcentaje
probabilidad_porc <- probabilidad_rango * 100

# 3. Resultado
cat("Probabilidad calculada:", round(probabilidad_porc, 2), "%")

## Probabilidad calculada: 28.06 %

10.1 Cálculo Gráfico de Probabilidades

# 1. Preparar áreas sombreadas
df_sombra_90 <- data.frame(x = seq(0, precip_90, length.out = 100))
df_sombra_90$y <- dexp(df_sombra_90$x, rate = lambda_exp) * A * 100

df_sombra_rango <- data.frame(x = seq(20, 60, length.out = 100))
df_sombra_rango$y <- dexp(df_sombra_rango$x, rate = lambda_exp) * A * 100

# 2. Curva completa
df_curva_final <- data.frame(x = seq(0, xmax, length.out = 300))
df_curva_final$y <- dexp(df_curva_final$x, rate = lambda_exp) * A * 100

# 3. Gráfica con LEYENDA
ggplot() +
  # Histograma base
  geom_bar(data = TDF, aes(x = MC, y = hi_porc), 
           stat = "identity", fill = "gray90", color = "black", alpha = 0.5, width = A) +
  
  # ÁREAS SOMBREADAS con mapeo para leyenda
  geom_ribbon(data = df_sombra_90, 
              aes(x = x, ymin = 0, ymax = y, fill = "Umbral No Extremo (90%)"), 
              alpha = 0.3) +
  
  geom_ribbon(data = df_sombra_rango, 
              aes(x = x, ymin = 0, ymax = y, fill = "Rango Probabilidad (20-60 mm)"), 
              alpha = 0.6) +
  
  # Línea del modelo
  geom_line(data = df_curva_final, aes(x = x, y = y), color = "darkred", size = 1.1) +
  
  # Configuración manual de colores de la leyenda
  scale_fill_manual(values = c("Umbral No Extremo (90%)" = "red", 
                               "Rango Probabilidad (20-60 mm)" = "#228B22")) +
  
  # Escalas y etiquetas
  scale_y_continuous(labels = function(x) paste0(x, "%"), expand = expansion(mult = c(0, 0.1))) +
  scale_x_continuous(breaks = sort(unique(c(0, 20, 60, round(precip_90, 1), round(xmax, 1))))) +
  
  labs(
    title = "Inferencia de probabilidad sobre el Modelo exponencial de Precipitación",
    x = "Precipitación (mm)", 
    y = "Densidad de probabilidad",
    fill = "Zonas de Inferencia",
    caption = "Fuente: Análisis basado en datos de VOLCAN A.csv"
  ) +
  
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

11 Teorema de límite central

# 1. Definir los datos a analizar (asegúrate de que esta variable exista)
# datos_analizar <- datos$Precipitation 

# Cálculo de estadísticos descriptivos básicos
n_total_p <- length(datos_analizar)           
x_bar_p   <- mean(datos_analizar, na.rm = TRUE)            
sd_p      <- sd(datos_analizar, na.rm = TRUE)                

# 2. Parámetros para el Teorema del Límite Central (Confianza 95%)
z_95 <- 1.96                                   
error_estandar_p <- sd_p / sqrt(n_total_p) 
margen_error_95_p <- z_95 * error_estandar_p

# 3. Estimación de los Intervalos de Confianza (Media Poblacional)
lim_inf_p <- x_bar_p - margen_error_95_p
lim_sup_p <- x_bar_p + margen_error_95_p

# 4. Creación del dataframe para la tabla resumen
tabla_precip_tlc <- data.frame(
  Parametro = "Precipitación Promedio Diaria",
  Lim_Inferior = round(lim_inf_p, 2),
  Media_Muestral = round(x_bar_p, 2),
  Lim_Superior = round(lim_sup_p, 2),
  Margen_Error = paste0("+/- ", round(margen_error_95_p, 2)),
  Confianza = "95% (Z=1.96)"
)

# 5. Generación de la tabla usando kable (evita el error de tab_style)
library(knitr)
kable(tabla_precip_tlc, 
      col.names = c("Parámetro", "Límite Inf. (mm)", "Media (mm)", "Límite Sup. (mm)", "Margen de Error", "Confianza"),
      caption = "ESTIMACIÓN DE LA MEDIA POBLACIONAL DE PRECIPITACIÓN (TLC)")

ESTIMACIÓN DE LA MEDIA POBLACIONAL DE PRECIPITACIÓN (TLC)

Parámetro	Límite Inf. (mm)	Media (mm)	Límite Sup. (mm)	Margen de Error	Confianza
Precipitación Promedio Diaria	15.45	17.1	18.76	+/- 1.65	95% (Z=1.96)

12 Conclusiones

La variable Precipitación, ha sido modelada con éxito mediante una Distribución Exponencial para eventos de alta frecuencia. Con una media aritmética muestral de 17.10 mm, definida por un margen de error de 1.65 mm.Mediante el Teorema del Límite Central (TLC), sabemos que la media aritmética poblacional de la precipitación diaria se encuentra entre [15.45 ; 18.76] mm con un 95% de confianza (\(Z=1.96\)). Esta inferencia estadística permite establecer umbrales de seguridad y previsiones hídricas sólidas, fundamentadas en un modelo matemático validado (\(\mu = 17.10 \pm 1.65\)).

2. VIENTO

1 Identificación y Justificación

La velocidad del viento se define como una variable cuantitativa continua, ya que se expresa mediante valores numéricos reales que pueden adoptar cualquier valor decimal dentro de un intervalo determinado. Esta característica de continuidad es fundamentalpermite registrar con precisión la magnitud del desplazamiento del aire, reflejando fielmente la variabilidad de las ráfagas y la intensidad en cada medición.

2 Datos

Importamos el archivo “VOLCAN A.csv” desde una ruta local y lo almacena en el objeto datos, usando espacios o tabulaciones como separador.

# setwd("C:/Users/ronal/OneDrive/Desktop")
datos <- read_csv("weatherdataANTISANA.csv")

## Rows: 366 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date
## dbl (9): Longitude, Latitude, Elevation, Max Temperature, Min Temperature, P...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(datos)

3 Extraer variable

Extraemos la variable viento, omitimos las celdas en blanco o valores iguales a cero y verificamos el tamaño muestral

viento <- na.omit(datos$Wind)
viento <- viento [viento > 0]
n_total <- length(viento)

4 Tabla de Frecuencia

En la tabla de distribución de frecuencias de la variable viento en el que el número de clases se determinó mediante la regla de Sturges y el ancho de clase se calculó a partir del rango total de los datos.

# 1. Definir el vector de datos
datos_analizar <- viento 

# 2. Cálculos base
n <- length(datos_analizar)
xmin <- min(datos_analizar)
xmax <- max(datos_analizar)
R <- xmax - xmin

# Regla de Sturges para calcular K 
K <- ceiling(1 + log2(n))

# Amplitud (A)
A <- R / K

# 3. Definición de Límites y Marcas de Clase (MC)
cortes <- seq(xmin, xmax, length.out = K + 1)
Li <- head(cortes, -1)
Ls <- tail(cortes, -1)
MC <- (Li + Ls) / 2

# 4. Cálculo de Frecuencias (ni)
ni <- as.vector(table(cut(datos_analizar, breaks = cortes, 
                          include.lowest = TRUE, right = FALSE)))

# 5. Frecuencias relativas y acumuladas
hi <- ni / sum(ni) * 100
Ni_asc <- cumsum(ni)
Ni_desc <- rev(cumsum(rev(ni)))
Hi_asc <- cumsum(hi)
Hi_desc <- rev(cumsum(rev(hi)))

# 6. Creación del Data Frame final
TDF <- data.frame(
  round(Li, 2), 
  round(Ls, 2), 
  round(MC, 2), 
  ni, 
  hi_porc = round(hi, 2), 
  Ni_asc, 
  Ni_desc, 
  Hi_asc_porc = round(Hi_asc, 2), 
  Hi_desc_porc = round(Hi_desc, 2)
)

# 7. Formato de tabla para el reporte
kable(TDF, 
      caption = paste("Tabla No. 1: Distribución de Frecuencias de la variable Viento "),
      col.names = c("Lím. Inf.", "Lím. Sup.", "Marca Clase", "ni", "hi (%)", 
                    "Ni Asc.", "Ni Desc.", "Hi Asc. (%)", "Hi Desc. (%)"),
      digits = 2)

Tabla No. 1: Distribución de Frecuencias de la variable Viento

Lím. Inf.	Lím. Sup.	Marca Clase	ni	hi (%)	Ni Asc.	Ni Desc.	Hi Asc. (%)	Hi Desc. (%)
0.59	0.83	0.71	2	0.55	2	366	0.55	100.00
0.83	1.07	0.95	18	4.92	20	364	5.46	99.45
1.07	1.31	1.19	43	11.75	63	346	17.21	94.54
1.31	1.55	1.43	67	18.31	130	303	35.52	82.79
1.55	1.79	1.67	63	17.21	193	236	52.73	64.48
1.79	2.03	1.91	61	16.67	254	173	69.40	47.27
2.03	2.27	2.15	51	13.93	305	112	83.33	30.60
2.27	2.51	2.39	41	11.20	346	61	94.54	16.67
2.51	2.75	2.63	18	4.92	364	20	99.45	5.46
2.75	2.99	2.87	2	0.55	366	2	100.00	0.55

5 Histograma

Una vez generada la Tabla de Distribución de Frecuencias, procedemos a visualizar los datos. Esta gráfica es fundamental para identificar la asimetría de la variable continua de precipitación y justificar el uso de modelos probabilísticos posteriores.

# 1. Crear el Histograma 
ggplot(TDF, aes(x = round.Li..2., y = hi_porc)) + 
  
  # Geometría de barras: Ancho definido por la amplitud 'A'
  geom_bar(stat = "identity", 
           fill = "skyblue3", 
           color = "black", 
           alpha = 0.8,
           width = A,
           just = 0) + 
  
  # Configuración del eje X: Marcas en los límites de clase
  scale_x_continuous(breaks = c(Li, Ls[K]), 
                     labels = function(x) round(x, 2)) +
  
  # Configuración del eje Y: Formato porcentual y ajuste de escala
  scale_y_continuous(labels = function(x) paste0(x, "%"),
                     limits = c(0, max(TDF$hi_porc) * 1.1),
                     expand = c(0, 0)) +
  
  # Etiquetas del gráfico
  labs(
    title = "Gráfica No. 1: Histograma de Distribución Porcentual del Viento",
    x = "Velocidad del Viento (m/s)",
    y = "Porcentaje (%)",
    caption = "Fuente: Datos analizados del Proyecto Volcán A"
  ) +
  
  # Estilo visual profesional
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 11, color = "darkgray"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

6 Conjetura del modelo LOG-NORMAL

Se seleccionó el modelo log-normal para la variable viento porque el histograma presenta una marcada asimetría positiva, donde las frecuencias no comienzan en el punto máximo, sino que muestran un crecimiento inicial rápido hasta alcanzar un pico y, posteriormente, un decaimiento gradual hacia la derecha.

# Se filtran datos mayores a cero para la transformación logarítmica
datos_log <- log(datos_analizar[datos_analizar > 0])
mu_log <- mean(datos_log, na.rm = TRUE)
sigma_log <- sd(datos_log, na.rm = TRUE)

# 2. Generación de los puntos para la curva teórica
x_curva <- seq(xmin, xmax, length.out = 300)

# Ajuste de escala: densidad * amplitud * 100 para eje porcentual
y_curva <- dlnorm(x_curva, meanlog = mu_log, sdlog = sigma_log) * A * 100
df_conjetura_log <- data.frame(x = x_curva, y = y_curva)

# 3. Gráfica Final: Histograma + Curva Log-Normal
ggplot(TDF, aes(x = MC, y = hi_porc)) +
  
  # Histograma base
  geom_bar(stat = "identity", 
           fill = "steelblue", 
           color = "black", 
           alpha = 0.7,
           width = A) + 
  
  # Superposición de la CURVA LOG-NORMAL
geom_line(data = df_conjetura_log, aes(x = x, y = y), 
          color = "darkgreen", linewidth = 1.2) + 
  
  # Configuración de ejes
  scale_x_continuous(breaks = c(Li, Ls[K]), 
                     labels = function(x) round(x, 2)) +
  
  scale_y_continuous(labels = function(x) paste0(x, "%"),
                     limits = c(0, max(c(TDF$hi_porc, df_conjetura_log$y), na.rm = TRUE) * 1.1),
                     expand = c(0, 0)) +
  
  # Etiquetas y títulos
  labs(
    title = "Gráfica No. 2: Modelo de probabilidad Log-normal para viento",
    subtitle = paste("Parámetros: meanlog (μ) =", round(mu_log, 3), 
                     " | sdlog (σ) =", round(sigma_log, 3)),
    x = "Velocidad del Viento (m/s)",
    y = "Porcentaje (%)",
  ) +
  
  theme_classic() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    plot.subtitle = element_text(hjust = 0.5, size = 10, face = "italic"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Test de PEARSON

Para determinar la validez de la conjetura Log-Normal, se comparan las frecuencias observadas en la muestra de viento contra las frecuencias esperadas. La cercanía entre ambos valores se cuantifica mediante el coeficiente de correlación de Pearson

# 1. Definir Frecuencia Observada (ni extraída de la TDF)
Fo <- ni 

# 2. Calcular Probabilidades Teóricas (P) bajo la curva Log-Normal
P <- numeric(K)
for (i in 1:K) {
  P[i] <- plnorm(cortes[i+1], meanlog = mu_log, sdlog = sigma_log) - 
          plnorm(cortes[i], meanlog = mu_log, sdlog = sigma_log)
}

# 3. Calcular Frecuencia Esperada (Fe)
n_total <- length(datos_analizar)
Fe <- P * n_total

# 4. Cálculo de Correlación de Pearson
Correlacion_Log <- cor(Fo, Fe) * 100

# 5. Salida de resultados de validación
cat("--- CONCLUSIÓN DE LA CONJETURA LOG-NORMAL (VIENTO) ---\n")

## --- CONCLUSIÓN DE LA CONJETURA LOG-NORMAL (VIENTO) ---

cat("Tamaño de la muestra total:", n_total, "registros\n")

## Tamaño de la muestra total: 366 registros

cat("Parámetros estimados: mu =", round(mu_log, 4), ", sigma =", round(sigma_log, 4), "\n")

## Parámetros estimados: mu = 0.5344 , sigma = 0.2725

cat("Correlación del modelo:", round(Correlacion_Log, 2), "%\n")

## Correlación del modelo: 95.73 %

8 Test de CHI-CUADRADO

Para confirmar si la distribución Log-Normal es un modelo estadísticamente válido para la velocidad del viento en el Volcán Antisana, se realiza la prueba de Chi-Cuadrado (\(\chi^2\)). Esta prueba evalúa si las diferencias entre las frecuencias observadas y las esperadas son lo suficientemente pequeñas como para ser atribuidas al azar o si, por el contrario, el modelo debe ser rechazado.

# 1. Determinación de los Grados de Libertad 
grados_libertad <- (K - 3) 

# Verificación de seguridad para grados de libertad
if(grados_libertad <= 0) grados_libertad <- 1

# 2. Nivel de significancia 
nivel_significancia <- 0.5

# 3. Preparación de Frecuencias Porcentuales
Fo_porc <- hi
Fe_porc <- P * 100

# 4. Cálculo del estadístico Chi-Cuadrado (X2)
x2 <- sum((Fe_porc - Fo_porc)^2 / (Fe_porc + 0.5))

# 5. Determinación del Umbral (Valor Crítico)
umbral_aceptacion <- qchisq(1 - nivel_significancia, grados_libertad)

# 6. Decisión Estadística
modelo_aceptado <- x2 < umbral_aceptacion

# --- SALIDA DE RESULTADOS ---
cat("--- PRUEBA DE CHI-CUADRADO PARA VIENTO (LOG-NORMAL) ---\n")

## --- PRUEBA DE CHI-CUADRADO PARA VIENTO (LOG-NORMAL) ---

cat("Grados de Libertad:", grados_libertad, "\n")

## Grados de Libertad: 7

cat("Estadístico Chi-Cuadrado (X2) calculado:", round(x2, 4), "\n")

## Estadístico Chi-Cuadrado (X2) calculado: 5.0692

cat("Umbral de aceptación (Valor Crítico):", round(umbral_aceptacion, 4), "\n")

## Umbral de aceptación (Valor Crítico): 6.3458

cat("¿El modelo es aceptado estadísticamente?:", ifelse(modelo_aceptado, "SÍ", "NO"), "\n")

## ¿El modelo es aceptado estadísticamente?: SÍ

9 Tabla Resumen

En esta sección se sintetizan los resultados de los estadísticos de validación. Se presentan el coeficiente de correlación de Pearson, que mide la relación lineal entre las frecuencias, y el estadístico Chi-Cuadrado, que determina la aceptación formal del modelo probabilístico.

# 1. Definición de las etiquetas de los modelos
Modelos <- c("Modelo Log-Normal")

# 2. Creación del Dataframe con los resultados obtenidos en bloques anteriores
tabla_resumen_viento <- data.frame(
  Modelo = Modelos,
  Pearson = c(round(Correlacion_Log, 2)),
  Chi_Sq = c(round(x2, 2)),
  Umbral = c(round(umbral_aceptacion, 2)),
  Resultado = c(ifelse(modelo_aceptado, "ACEPTADO", "RECHAZADO"))
)

# 3. Renombrar columnas para un formato de reporte profesional
colnames(tabla_resumen_viento) <- c("Modelo de Probabilidad", 
                                   "Correlación Pearson (%)", 
                                   "Chi Cuadrado (X2)", 
                                   "Umbral de Aceptación",
                                   "Decisión Final")

# 4. Generación de la tabla con kable
library(knitr)
kable(tabla_resumen_viento, 
      format = "markdown", 
      align = "lcccc",
      caption = "Tabla Nro. 2: Resumen de test de bondad para el modelo log-normal del Viento")

Tabla Nro. 2: Resumen de test de bondad para el modelo log-normal del Viento

Modelo de Probabilidad	Correlación Pearson (%)	Chi Cuadrado (X2)	Umbral de Aceptación	Decisión Final
Modelo Log-Normal	95.73	5.07	6.35	ACEPTADO

10 Cálculo de probabilidades

Una vez validado el modelo log-normal, procedemos a resolver interrogantes clave sobre el comportamiento futuro del viento en el Volcán Antisana, utilizando las funciones de distribución de probabilidad.

1). Si se toman 200 mediciones futuras, ¿cuántas veces se espera que el viento sople con una velocidad entre 1.79 y 2.27 (m/s)?

# Parámetros extraídos de la conjetura del modelo
mu_viento <- 0.534
sigma_viento <- 0.273

# 1. Cálculo de la probabilidad del intervalo [1.79, 2.27]
P_intervalo_cant <- plnorm(2.27, meanlog = mu_viento, sdlog = sigma_viento) - 
                    plnorm(1.79, meanlog = mu_viento, sdlog = sigma_viento)

# 2. Proyección para n = 200 mediciones
n_futuro <- 200
cantidad_esperada <- P_intervalo_cant * n_futuro

cat("Resultado: Se esperan aproximadamente", round(cantidad_esperada, 0), 
    "mediciones en el rango [1.79 - 2.27].\n")

## Resultado: Se esperan aproximadamente 56 mediciones en el rango [1.79 - 2.27].

2). ¿Cuál será la probabilidad de que la velocidad del viento se encuentre entre 1.3 y 2.3 m/s?

# 1. Definición de los límites del intervalo de interés
lim_inf <- 1.3
lim_sup <- 2.3

# 2. Cálculo de la probabilidad acumulada P(1.3 <= X <= 2.3)
prob_rango <- plnorm(lim_sup, meanlog = mu_viento, sdlog = sigma_viento) - 
              plnorm(lim_inf, meanlog = mu_viento, sdlog = sigma_viento)

# 3. Mostrar resultado porcentual
cat("Resultado: La probabilidad de que el viento esté entre 1.3 y 2.3 unidades es del", 
    round(prob_rango * 100, 2), "%\n")

## Resultado: La probabilidad de que el viento esté entre 1.3 y 2.3 unidades es del 70.34 %

10.1 Cálculo gráfico de probabilidades

Esta sección presenta, una síntesis visual que superpone el modelo Log-Normal sobre el histograma de frecuencias observadas del Volcán Antisana. Se resaltan dos áreas críticas.

### 8. Gráfica Maestra de Inferencia Probabilística

# 1. Parámetros y cálculos de probabilidad para las preguntas
mu_viento <- 0.534
sigma_viento <- 0.273

# Pregunta 1: Rango [1.79 - 2.27]
p_pregunta1 <- plnorm(2.27, meanlog = mu_viento, sdlog = sigma_viento) - 
               plnorm(1.79, meanlog = mu_viento, sdlog = sigma_viento)
cantidad_200 <- p_pregunta1 * 200

# Pregunta 2: Rango [1.3 - 2.3]
p_pregunta2 <- plnorm(2.3, meanlog = mu_viento, sdlog = sigma_viento) - 
               plnorm(1.3, meanlog = mu_viento, sdlog = sigma_viento)

# 2. Preparar datos para el sombreado (geom_ribbon)
# Área para Pregunta 1 (Naranja)
df_sombra_p1 <- data.frame(x = seq(1.79, 2.27, length.out = 100))
df_sombra_p1$y <- dlnorm(df_sombra_p1$x, meanlog = mu_viento, sdlog = sigma_viento) * A * 100

# Área para Pregunta 2 (Verde)
df_sombra_p2 <- data.frame(x = seq(1.3, 2.3, length.out = 100))
df_sombra_p2$y <- dlnorm(df_sombra_p2$x, meanlog = mu_viento, sdlog = sigma_viento) * A * 100

# Curva completa
df_curva_final <- data.frame(x = seq(min(datos_analizar), max(datos_analizar), length.out = 300))
df_curva_final$y <- dlnorm(df_curva_final$x, meanlog = mu_viento, sdlog = sigma_viento) * A * 100

# 3. Generación del Gráfico con Leyenda Manual
library(ggplot2)

ggplot() +
  # Histograma de fondo
  geom_bar(data = TDF, aes(x = MC, y = hi_porc), 
           stat = "identity", fill = "steelblue", color = "black", alpha = 0.1, width = A) +
  
  # Sombreados con 'fill' dentro de aes para generar la leyenda
  geom_ribbon(data = df_sombra_p2, aes(x = x, ymin = 0, ymax = y, fill = "Rango 1.3 - 2.3 (70.34%)"), 
              alpha = 0.4) +
  geom_ribbon(data = df_sombra_p1, aes(x = x, ymin = 0, ymax = y, fill = "Rango 1.79 - 2.27 (28.23%)"), 
              alpha = 0.6) +
  
  # Línea del modelo
  geom_line(data = df_curva_final, aes(x = x, y = y), color = "darkgreen", linewidth = 1.2) +
  
  # Configuración de colores de la leyenda
  scale_fill_manual(name = "Zonas de Inferencia", 
                    values = c("Rango 1.3 - 2.3 (70.34%)" = "forestgreen", 
                               "Rango 1.79 - 2.27 (28.23%)" = "orange")) +
  
  # Escalas y Etiquetas
  scale_y_continuous(labels = function(x) paste0(x, "%"), expand = expansion(mult = c(0, 0.1))) +
  scale_x_continuous(breaks = sort(c(1.3, 1.79, 2.27, 2.3))) +
  
  labs(title = "Inferencia Probabilística: Velocidad del Viento",
       subtitle = "Análisis de rangos operativos y frecuencias esperadas",
       x = "Velocidad del Viento (m/s)", 
       y = "DENSIDAD DE PROBABILIDAD",
       caption = "Modelo Log-Normal aplicado a Datos Volcán A") +
  
  theme_minimal() +
  theme(legend.position = "bottom",
        plot.title = element_text(face = "bold", size = 14, hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))

# 4. Salida de resultados (Responde a las dos preguntas)
cat("--- RESULTADOS DE INFERENCIA ---\n")

## --- RESULTADOS DE INFERENCIA ---

cat("Pregunta 1: En 200 mediciones, se esperan", round(cantidad_200, 0), 
    "veces vientos entre 1.79 y 2.27 m/s.\n")

## Pregunta 1: En 200 mediciones, se esperan 56 veces vientos entre 1.79 y 2.27 m/s.

cat("Pregunta 2: La probabilidad de vientos entre 1.3 y 2.3 m/s es del", 
    round(p_pregunta2 * 100, 2), "%.\n")

## Pregunta 2: La probabilidad de vientos entre 1.3 y 2.3 m/s es del 70.34 %.

11 Teorema de límite central

A continuación, se aplica el Teorema del Límite Central para estimar el valor promedio real de la velocidad del viento con un nivel de confianza del 95%. Esta técnica estadística permite inferir que, si se realizaran múltiples muestreos, el promedio poblacional se encontraría dentro del intervalo calculado en la gran mayoría de los casos.

# Asegúrate de haber definido la variable viento antes:
viento <- datos$Wind

# 1. Cálculo de estadísticos descriptivos básicos
n_total_v <- length(viento)          
x_bar_v   <- mean(viento, na.rm = TRUE)           
sd_v      <- sd(viento, na.rm = TRUE)                

# 2. Parámetros para el Teorema del Límite Central (Confianza 95%)
z_95 <- 1.96                                       
error_estandar_v <- sd_v / sqrt(n_total_v) 
margen_error_95_v <- z_95 * error_estandar_v

# 3. Estimación de los Intervalos de Confianza (Media Poblacional)
lim_inf_v <- x_bar_v - margen_error_95_v
lim_sup_v <- x_bar_v + margen_error_95_v

# 4. Creación del dataframe para la tabla resumen
tabla_viento_tlc <- data.frame(
  Parametro = "Velocidad Promedio del Viento",
  Lim_Inferior = round(lim_inf_v, 2),
  Media_Muestral = round(x_bar_v, 2),
  Lim_Superior = round(lim_sup_v, 2),
  Margen_Error = paste0("+/- ", round(margen_error_95_v, 2)),
  Confianza = "95% (Z=1.96)"
)

# 5. Generación de la tabla usando kable (Librería knitr)
library(knitr)
kable(tabla_viento_tlc, 
      col.names = c("Parámetro", "Límite Inf. (m/s)", "Promedio (m/s)", "Límite Sup. (m/s)", "Margen de Error", "Confianza"),
      caption = "ESTIMACIÓN DE LA MEDIA POBLACIONAL DE VIENTO (TLC)",
      align = "lcccc")

ESTIMACIÓN DE LA MEDIA POBLACIONAL DE VIENTO (TLC)

Parámetro	Límite Inf. (m/s)	Promedio (m/s)	Límite Sup. (m/s)	Margen de Error	Confianza
Velocidad Promedio del Viento	1.72	1.77	1.81	+/- 0.05	95% (Z=1.96)

Conclusiones

La variable Viento ha sido modelada con éxito mediante una Distribución Log-Normal. Con una media aritmética muestral de 1.77 m/s, definida por un margen de error de 0.05 m/s.Mediante la aplicación del Teorema del Límite Central, se ha determinado la media aritmética poblacional de la velocidad del viento se encuentra en el intervalo de [1.72 ; 1.81] m/s con un 95% de confianza (Z=1.96). Esta inferencia estadística es fundamental, permitiendo establecer un modelo matemático validado (\(\mu = 1.77 \pm 0.05\)).

1. REGRECIONES

1.REGRECIONES LINEAL

1 CARGAR DATOS

library(readr)
library(ggplot2)
library(knitr)
library(dplyr)

# Cargamos el archivo directamente como 'datos'
datos <- read_csv("weatherdataANTISANA.csv", show_col_types = FALSE)
# View(datos) # Descomentar si necesitas ver la ventana

2 EXTRAER VARIABLES

## 2 EXTRAER VARIABLES (Sin limpiar todavía para mantener simetría)
# Nota: Quitamos el na.omit individual de aquí para evitar el error de "diferente número de filas"
longitud <- datos$Longitude
latitud <- datos$Latitude
elevacion <- datos$Elevation
temp_maxima <- datos$`Max Temperature`  # Usamos comillas invertidas por el espacio
temp_minima <- datos$`Min Temperature`
precipitacion <- datos$Precipitation
viento <- datos$Wind
humedad_relativa <- datos$`Relative Humidity`
radiacion_solar <- datos$Solar

3 SELECIONAR VARIABLES

# Asignamos variables para el modelo
x_raw <- datos$Solar
y_raw <- datos$`Max Temperature`

4 TABLA DE VALORES

# 3. Crear la tabla de pares de valores usando backticks para el espacio
tabla_regresion <- data.frame(
  Radiacion_Solar = datos$Solar,
  Temp_Maxima = datos$`Max Temperature`
)

# 4. Limpiar filas con valores vacíos (NA)
tabla_regresion <- na.omit(tabla_regresion)

# 5. Configurar R para mostrar las filas necesarias
options(max.print = 1000)

# 6. Mostrar la tabla
print("Tabla de pares de valores (Radiación vs Temp Máxima):")

## [1] "Tabla de pares de valores (Radiación vs Temp Máxima):"

head(tabla_regresion) # Mostramos el inicio para no saturar el reporte

# Actualizamos x e y con los datos ya limpios y pareados
x <- tabla_regresion$Radiacion_Solar
y <- tabla_regresion$Temp_Maxima

cat("Total de filas procesadas:", nrow(tabla_regresion))

## Total de filas procesadas: 366

5 GRÁFICA

plot(x = tabla_regresion$Radiacion_Solar, 
     y = tabla_regresion$Temp_Maxima,
     main = "Dispersión: Temp. Máxima vs Radiación Solar",
     xlab = expression("Radiación Solar (MJ/m"^2*"/día)"), 
     ylab = expression("Temperatura Máxima ("*degree*"C)"), 
     pch = 19, col = "steelblue", cex = 0.8)

6 CONJETURA DEL MODELO LINEAL

# Ajustar el modelo usando los datos limpios de la tabla
modelo <- lm(Temp_Maxima ~ Radiacion_Solar, data = tabla_regresion)

7 CÁLCULO DE PARAMETROS

b0 <- coef(modelo)[1] 
b1 <- coef(modelo)[2] 

cat("--- Parámetros del Modelo ---\n")

## --- Parámetros del Modelo ---

cat("Intercepto (b0): ", round(b0, 4), "\n")

## Intercepto (b0):  11.2237

cat("Pendiente (b1): ", round(b1, 4), "\n")

## Pendiente (b1):  0.3129

cat("Ecuación: Y =", round(b0, 4), "+", round(b1, 4), "* X\n")

## Ecuación: Y = 11.2237 + 0.3129 * X

summary(modelo)

## 
## Call:
## lm(formula = Temp_Maxima ~ Radiacion_Solar, data = tabla_regresion)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5457 -0.8530 -0.0284  0.7079  3.4931 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     11.223741   0.125514   89.42   <2e-16 ***
## Radiacion_Solar  0.312879   0.007534   41.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.199 on 364 degrees of freedom
## Multiple R-squared:  0.8257, Adjusted R-squared:  0.8253 
## F-statistic:  1725 on 1 and 364 DF,  p-value: < 2.2e-16

8 GRÁFICA DE MODELO DE REGRESIÓN LÍNEAL

plot(tabla_regresion$Radiacion_Solar, tabla_regresion$Temp_Maxima,
     main = "Modelo de Regresión Lineal - Antisana",
     xlab = expression("Radiación Solar (MJ/m"^2*"/día)"), 
     ylab = expression("Temperatura Máxima ("*degree*"C)"), 
     pch = 19, col = "steelblue")

abline(modelo, col = "red", lwd = 2)

9 TEST DE PEARSON

resultado_test <- cor.test(x, y, method = "pearson")
print(resultado_test)

## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 41.531, df = 364, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8890033 0.9250430
## sample estimates:
##       cor 
## 0.9087017

r <- resultado_test$estimate
cat("\nEl coeficiente de correlación de Pearson (r) es:", round(r, 4), "\n")

## 
## El coeficiente de correlación de Pearson (r) es: 0.9087

10 COEFICIENTE DE DETERMINACIÓN

resumen <- summary(modelo)
r_cuadrado <- resumen$r.squared

cat("--- Coeficiente de Determinación ---\n")

## --- Coeficiente de Determinación ---

cat("R-squared (R2):", round(r_cuadrado, 4), "\n")

## R-squared (R2): 0.8257

cat("En porcentaje:", round(r_cuadrado * 100, 2), "%\n")

## En porcentaje: 82.57 %

cat("Interpretación: El", round(r_cuadrado * 100, 2), 
    "% de la variación en la Temp. Máxima se explica por la Radiación Solar.")

## Interpretación: El 82.57 % de la variación en la Temp. Máxima se explica por la Radiación Solar.

11 ESTIMACIONES

¿Cuál sería la temperatura máxima probable ante diferentes niveles de radiación solar 5, 15, 25, 30 (MJ/m”^2*“/día) ?

# Valores a predecir
pregunta_sol <- data.frame(Radiacion_Solar = c(5, 15, 25, 30)) 

# Predicción usando el modelo y la tabla de regresión
respuestas <- predict(modelo, newdata = pregunta_sol)

# Resultados
for(i in 1:nrow(pregunta_sol)) {
  cat("Si la radiación es", pregunta_sol$Radiacion_Solar[i], "MJ/m²/día,",
      "la temperatura estimada es", round(respuestas[i], 2), "°C\n")
}

## Si la radiación es 5 MJ/m²/día, la temperatura estimada es 12.79 °C
## Si la radiación es 15 MJ/m²/día, la temperatura estimada es 15.92 °C
## Si la radiación es 25 MJ/m²/día, la temperatura estimada es 19.05 °C
## Si la radiación es 30 MJ/m²/día, la temperatura estimada es 20.61 °C

12 CONCLUSIONES

La conjetura del modelo lineal confirma que la temperatura atmosférica en el Volcán ANTISANA no fluctúa de forma aleatoria, sino que presenta una dependencia crítica y directa de la insolación diaria. Mientras que en niveles de radiación bajos la temperatura se mantiene en un rango base (cercano a los 11.22 °C). Esto demuestra que la radiación solar es el motor termodinámico dominante en la zona, donde cada incremento en la intensidad solar impacta de forma inmediata en el balance calórico del sitio, validando la alta fiabilidad predictiva de nuestro modelo (R² ≈ 82%).

2. REGRCION POTENCIAL

Datos

datasetf <- read_csv("weatherdataANTISANA.csv", show_col_types = FALSE)

2. SELECCIÓN DE DATOS

df_pot <- data.frame(
  Viento = datasetf$Wind,
  Radiacion = datasetf$Solar
)

# 3. LIMPIEZA CRÍTICA PARA MODELO POTENCIAL

df_pot <- na.omit(df_pot)
df_pot <- df_pot[df_pot$Viento > 0 & df_pot$Radiacion > 0, ]

4. MODELO POTENCIAL

# Aplicamos logaritmo natural (log) a AMBAS variables: log(Y) ~ log(X)
modelo_pot <- lm(log(Radiacion) ~ log(Viento), data = df_pot)

5. CÁLCULO DE PARÁMETROS (Y = a * X^b)

a <- exp(coef(modelo_pot)[1]) # Constante de proporcionalidad (intercepto transformado)
b <- coef(modelo_pot)[2]      # Exponente de potencia (elasticidad)

# 6. GRÁFICA DEL MODELO

plot(df_pot$Viento, df_pot$Radiacion, 
     main="Regresión Potencial: Viento vs Radiación Solar",
     xlab="Velocidad del Viento (m/s)", ylab="Radiación Solar (MJ/m2)",
     pch=16, col="darkorange")

# Añadimos la curva potencial predicha ( Y = a * X^b )
x_seq <- seq(min(df_pot$Viento), max(df_pot$Viento), length.out=100)
y_pred <- a * (x_seq ^ b)
lines(x_seq, y_pred, col="blue", lwd=3)

7. ESTIMACIÓN Y RESULTADOS EN CONSOLA

viento_test <- 3.5
prediccion_rad <- a * (viento_test ^ b)

cat("\n--- RESULTADOS DEL MODELO POTENCIAL ---\n")

## 
## --- RESULTADOS DEL MODELO POTENCIAL ---

cat("Exponente o Tasa de Escala (b):", round(b, 4), "\n")

## Exponente o Tasa de Escala (b): 1.7983

cat("Constante inicial (a):", round(a, 4), "\n")

## Constante inicial (a): 4.4622

cat("Estimación para un viento de 3.5 m/s: ", round(prediccion_rad, 2), "MJ/m2\n\n")

## Estimación para un viento de 3.5 m/s:  42.46 MJ/m2

cat("--- RESUMEN ESTADÍSTICO ---\n")

## --- RESUMEN ESTADÍSTICO ---

print(summary(modelo_pot))

## 
## Call:
## lm(formula = log(Radiacion) ~ log(Viento), data = df_pot)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2556 -0.1941  0.1224  0.3296  1.0350 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.49565    0.05986   24.99   <2e-16 ***
## log(Viento)  1.79827    0.09982   18.02   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5197 on 364 degrees of freedom
## Multiple R-squared:  0.4714, Adjusted R-squared:  0.4699 
## F-statistic: 324.6 on 1 and 364 DF,  p-value: < 2.2e-16

Concluciones

modelo de potencias (Y = a * X^b) confirma que a medida que la velocidad del viento incrementa, la radiación solar incidente crece de forma no lineal (potencial). El exponente ‘b’ positivo demuestra que las ráfagas fuertes barren la nubosidad en el Antisana, permitiendo que la radiación solar penetre con mayor intensidad. Al simular vientos fuertes de 3.5 m/s, se estima que la radiación superará los 42 MJ/m2.

3. REGRECION EXPONENCIAL

datasetf <- read_csv("weatherdataANTISANA.csv", show_col_types = FALSE)

2. SELECCIÓN Y LIMPIEZA DE DATOS

# X = Radiación Solar, Y = Temperatura Máxima
df_exp <- data.frame(
  Solar = datasetf$Solar,
  Max_Temp = datasetf$`Max Temperature`
)
# Filtramos valores mayores a 0 para poder aplicar logaritmos
df_exp <- na.omit(df_exp)
df_exp <- df_exp[df_exp$Max_Temp > 0 & df_exp$Solar > 0, ]

3. MODELO EXPONENCIAL

# Aplicamos logaritmo natural (log) a la variable dependiente
modelo_exp <- lm(log(Max_Temp) ~ Solar, data = df_exp)

# 4. CÁLCULO DE PARÁMETROS (Y = a * e^(bX))
a <- exp(coef(modelo_exp)[1]) # Intercepto transformado
b <- coef(modelo_exp)[2]      # Tasa de crecimiento

4. GRÁFICA DEL MODELO

plot(df_exp$Solar, df_exp$Max_Temp, 
     main="Regresión Exponencial: Radiación vs Temperatura",
     xlab="Radiación Solar (MJ/m2)", ylab="Temperatura Máxima (°C)",
     pch=16, col="blue") # Cambié a color azul sólido para evitar errores gráficos

# Añadimos la curva exponencial predicha inmediatamente
x_seq <- seq(min(df_exp$Solar), max(df_exp$Solar), length.out=100)
y_pred <- a * exp(b * x_seq)
lines(x_seq, y_pred, col="red", lwd=3)

5. ESTIMACIÓN Y RESULTADOS

radiacion_test <- 25
prediccion_temp <- a * exp(b * radiacion_test)

cat("\n--- RESULTADOS DEL MODELO ---\n")

## 
## --- RESULTADOS DEL MODELO ---

cat("Coeficiente de crecimiento (b):", round(b, 4), "\n")

## Coeficiente de crecimiento (b): 0.0198

cat("Valor inicial (a):", round(a, 4), "\n")

## Valor inicial (a): 11.6414

cat("Estimación para una radiación de 25: ", round(prediccion_temp, 2), "°C\n")

## Estimación para una radiación de 25:  19.08 °C

conclusion

La estimación realizada confirma que el incremento de temperatura debido a la radiación solar no escala de forma lineal, sino acelerada. Mientras que en niveles bajos de radiación las temperaturas son estables, una vez que la radiación supera un umbral crítico (en este caso, las 25 unidades estimadas), la temperatura explota debido a la alta acumulación térmica en el ambiente.

4.REGRECION LOGARIDMICA

Datos

datasetf <- read_csv("weatherdataANTISANA.csv", show_col_types = FALSE)

2. SELECCIÓN DE DATOS

df_log <- data.frame(
  Precipitacion = datasetf$Precipitation,
  Humedad = datasetf$`Relative Humidity`
)

3. LIMPIEZA CRÍTICA PARA MODELO LOGARÍTMICO

df_log <- na.omit(df_log)
df_log <- df_log[df_log$Precipitacion > 0, ]

4. MODELO LOGARÍTMICO

# Aplicamos logaritmo natural SOLO a la variable independiente (X): Y ~ log(X)
modelo_log <- lm(Humedad ~ log(Precipitacion), data = df_log)

5. CÁLCULO DE PARÁMETROS (Y = a + b * ln(X))

# OJO: En este modelo NO se usa exp(a), se extrae directo.
a <- coef(modelo_log)[1] # Intercepto (punto de base)
b <- coef(modelo_log)[2] # Pendiente logarítmica (tasa de desaceleración)

6. GRÁFICA DEL MODELO

plot(df_log$Precipitacion, df_log$Humedad, 
     main="Regresión Logarítmica: Precipitación vs Humedad",
     xlab="Precipitación o Lluvia (mm)", ylab="Humedad Relativa (0.0 a 1.0)",
     pch=16, col="darkgreen")

# Añadimos la curva logarítmica predicha ( Y = a + b * ln(X) )
x_seq <- seq(min(df_log$Precipitacion), max(df_log$Precipitacion), length.out=100)
y_pred <- a + (b * log(x_seq))
lines(x_seq, y_pred, col="red", lwd=3)

7. ESTIMACIÓN

precip_test <- 50
prediccion_humedad <- a + (b * log(precip_test))

cat("\n--- RESULTADOS DEL MODELO LOGARÍTMICO ---\n")

## 
## --- RESULTADOS DEL MODELO LOGARÍTMICO ---

cat("Constante inicial (a):", round(a, 4), "\n")

## Constante inicial (a): 0.8195

cat("Tasa de cambio (b):", round(b, 4), "\n")

## Tasa de cambio (b): 0.0413

cat("Estimación de humedad relativa para lluvia fuerte de 50 mm: ", round(prediccion_humedad * 100, 2), "%\n\n")

## Estimación de humedad relativa para lluvia fuerte de 50 mm:  98.1 %

cat("--- RESUMEN ESTADÍSTICO ---\n")

## --- RESUMEN ESTADÍSTICO ---

print(summary(modelo_log))

## 
## Call:
## lm(formula = Humedad ~ log(Precipitacion), data = df_log)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.21689 -0.02171  0.01158  0.03625  0.18052 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.819540   0.003910  209.58   <2e-16 ***
## log(Precipitacion) 0.041271   0.001357   30.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05781 on 364 degrees of freedom
## Multiple R-squared:  0.7176, Adjusted R-squared:  0.7169 
## F-statistic: 925.1 on 1 and 364 DF,  p-value: < 2.2e-16

Concluciones

l modelo logarítmico (Y = a + b * ln(X)) describe a la perfección el comportamiento de la humedad relativa. Los primeros milímetros de lluvia disparan la humedad ambiente velozmente, pero a medida que llueve más fuerte (como en nuestra estimación de 50 mm), el incremento se desacelera hasta estancarse contra su límite natural de saturación atmosférica (1.0 o 100%). La curva roja refleja matemáticamente ese ‘techo’ que impide a la humedad crecer infinitamente.

5. REGRECION MULTIPLE

Datos

datasetf <- read_csv("weatherdataANTISANA.csv", show_col_types = FALSE)

2. SELECCIÓN DE MÚLTIPLES VARIABLES

# Z (Eje Vertical) = Temperatura Máxima (Variable a predecir)
# X (Eje Horizontal) = Radiación Solar (Aporta calor)
# Y (Eje Profundidad) = Humedad Relativa (Aporta enfriamiento por nubosidad)
df_mult <- data.frame(
  Max_Temp = datasetf$`Max Temperature`,
  Solar    = datasetf$Solar,
  Humedad  = datasetf$`Relative Humidity`
)

3. LIMPIEZA DE DATOS

# Eliminamos cualquier fila que tenga datos faltantes (NA)
df_mult <- na.omit(df_mult)

4. MODELO DE REGRESIÓN MÚLTIPLE

# Ecuación: Temp = b0 + b1*Solar + b2*Humedad
modelo_multiple <- lm(Max_Temp ~ Solar + Humedad, data = df_mult)

5. CÁLCULO DE PARÁMETROS

b0 <- coef(modelo_multiple)[1] # Intercepto base (beta 0)
b1 <- coef(modelo_multiple)[2] # Peso de la Radiación Solar (beta 1)
b2 <- coef(modelo_multiple)[3] # Peso de la Humedad (beta 2)

6. GRÁFICA DEL MODELO EN 3D

library(scatterplot3d)
# Dibujamos primero la caja 3D y los puntos de los datos reales flotando en el espacio
grafico_3d <- scatterplot3d(x = df_mult$Solar, 
                            y = df_mult$Humedad, 
                            z = df_mult$Max_Temp,
                            pch = 16, color = "steelblue",
                            main = "Plano de Regresión Múltiple 3D",
                            xlab = "Radiación Solar (MJ/m2)",
                            ylab = "Humedad Relativa",
                            zlab = "Temperatura Máxima (°C)",
                            angle = 45, # Ángulo de rotación del cubo
                            grid = TRUE, box = TRUE)

# Le añadimos el "Plano de Regresión" que atraviesa los puntos
# Este plano representa la fórmula b0 + b1*X1 + b2*X2
grafico_3d$plane3d(modelo_multiple, col = "red", lty.box = "solid", lwd = 2)

7. ESTIMACIÓN

solar_test <- 20
humedad_test <- 0.90
estimacion_temp <- b0 + (b1 * solar_test) + (b2 * humedad_test)

cat("\n--- RESULTADOS DEL MODELO MÚLTIPLE ---\n")

## 
## --- RESULTADOS DEL MODELO MÚLTIPLE ---

cat("Intercepto (b0):", round(b0, 4), "\n")

## Intercepto (b0): 18.83

cat("Efecto del Sol (b1):", round(b1, 4), "°C extra por cada MJ/m2 de sol\n")

## Efecto del Sol (b1): 0.2297 °C extra por cada MJ/m2 de sol

cat("Efecto de la Humedad (b2):", round(b2, 4), "°C menos por cada punto de humedad\n")

## Efecto de la Humedad (b2): -7.1563 °C menos por cada punto de humedad

cat("Estimación para un día con 20 de Sol y 90% de Humedad: ", round(estimacion_temp, 2), "°C\n\n")

## Estimación para un día con 20 de Sol y 90% de Humedad:  16.98 °C

cat("--- RESUMEN ESTADÍSTICO ---\n")

## --- RESUMEN ESTADÍSTICO ---

print(summary(modelo_multiple))

## 
## Call:
## lm(formula = Max_Temp ~ Solar + Humedad, data = df_mult)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1425 -0.8037 -0.0192  0.7258  3.7129 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.83003    1.30153  14.468  < 2e-16 ***
## Solar        0.22967    0.01591  14.440  < 2e-16 ***
## Humedad     -7.15626    1.21930  -5.869 9.88e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.147 on 363 degrees of freedom
## Multiple R-squared:  0.8408, Adjusted R-squared:   0.84 
## F-statistic: 958.9 on 2 and 363 DF,  p-value: < 2.2e-16

CONCLUSIÓN

El gráfico en 3D muestra la verdadera naturaleza de la regresión múltiple: ya no es una línea que intenta ajustarse, sino un ‘plano’ (superficie roja) suspendido en 3 dimensiones. Los puntos azules (datos reales) gravitan muy cerca de este plano, lo cual demuestra que el modelo es altamente preciso. Confirmamos visual y matemáticamente que la Temperatura Máxima (eje vertical Z) sube a medida que avanzamos en el eje de Radiación Solar, pero decae abruptamente si avanzamos por el eje de Humedad.