Carga de Datos y Librerías

1. CARGA DE LIBRERIAS

La variable año de inicio es una variable discreta pero debido a su gran cantidad de diferentes años se decidio agruparlos en intervalos, La cual la trabajaremos como una variable continua

# cargar librerías
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(gt)
library(knitr)

1.2. CARGA DE DATOS 

#cargar datos
datos <- read.csv("C:\\Users\\joeja\\Desktop\\Proyecto Estadística\\Depositos_sulfuro.csv", 
                  header = TRUE, 
                  sep = ";", 
                  dec = ".")

Tabla de distribucion de probabilidad

2. TABLAS DE DISTRIBUCIÓN DE FRECUENCIA 

# LIMPIEZA DE LA VARIABLE AÑO DE INICIO
año_inicio <- as.numeric(datos$startdate)
## Warning: NAs introducidos por coerción
año_inicio <- na.omit(año_inicio)

# SEPARAR OUTLIERS
caja <- boxplot(año_inicio, plot = FALSE)

limite_sup <- caja$stats[5]
limite_inf <- caja$stats[1]

año_outliers <- año_inicio[año_inicio < limite_inf | año_inicio > limite_sup]
año_sin_outliers <- año_inicio[año_inicio >= limite_inf & año_inicio <= limite_sup]


# RESUMEN
cat("Cantidad con outliers:", length(año_inicio), "\n")
## Cantidad con outliers: 518
cat("Cantidad de outliers:", length(año_outliers), "\n")
## Cantidad de outliers: 42
cat("Cantidad sin outliers:", length(año_sin_outliers), "\n")
## Cantidad sin outliers: 476
#TABLA DE DISTRIBUCION DE FRECUENCIA

#Crear un histograma para extraer información
histograma <- hist(año_sin_outliers,
                   plot = FALSE)
# Frecuencia absoluta (ni)
ni <- histograma$counts
# Frecuencia relativa (hi)
hi <- ni / sum(ni)*100
# Intervalos
intervalos <- paste0(
  "[", round(histograma$breaks[-length(histograma$breaks)], 2),
  ", ",
  round(histograma$breaks[-1], 2),
  ")"
)

# TABLA FINAL
tabla_frecuencias <- data.frame(
  Intervalo = intervalos,
  ni = ni,
  hi = round(hi, 2)
)

# Mostrar la tabla
tabla_frecuencias
##       Intervalo  ni    hi
## 1  [1800, 1820)   4  0.84
## 2  [1820, 1840)   3  0.63
## 3  [1840, 1860)   8  1.68
## 4  [1860, 1880)  31  6.51
## 5  [1880, 1900)  38  7.98
## 6  [1900, 1920)  89 18.70
## 7  [1920, 1940)  60 12.61
## 8  [1940, 1960)  76 15.97
## 9  [1960, 1980) 109 22.90
## 10 [1980, 2000)  53 11.13
## 11 [2000, 2020)   5  1.05

TABLA DE DISTRIBUCIÓN DE PROBABILIDADES

# Mejorar la Tabla
tabla_Añoinicio_gt <- tabla_frecuencias %>%
  gt() %>%
  tab_header(
    title = md("**Tabla N° 1**"),
    subtitle = md("**Distribución de probabilidad de los Depósitos Masivos<br>
                     de Sulfuros Volcánicos por Año de inicio**")
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 2")
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black",
    row.striping.include_table_body = TRUE
  )

tabla_Añoinicio_gt
Tabla N° 1
Distribución de probabilidad de los Depósitos Masivos
de Sulfuros Volcánicos por Año de inicio
Intervalo ni hi
[1800, 1820) 4 0.84
[1820, 1840) 3 0.63
[1840, 1860) 8 1.68
[1860, 1880) 31 6.51
[1880, 1900) 38 7.98
[1900, 1920) 89 18.70
[1920, 1940) 60 12.61
[1940, 1960) 76 15.97
[1960, 1980) 109 22.90
[1980, 2000) 53 11.13
[2000, 2020) 5 1.05
Autor: Grupo 2

Gráficas de distribución de probabilidad

#Gráfica de la variable
histograma<-hist(año_sin_outliers,
                 freq = FALSE,
                 main="Gráfica 1. Distribucion de densidad
                 de probabilidad de los años de 
                 inicio",
                 xlab="Años de inicio",
                 ylab="Densidad de probabilidad",
                 col="blue")

## Conjetura del modelo Debido a la similitud de las barras asociamos con el modelo de probabilidad log-normal

#Gráfica de modelo de probabilidad Log-normal
histograma<-hist(año_sin_outliers,
                 freq = FALSE,
                 main="Gráfica 2.Comparación de la realidad 
                 con el modelo log-normal de los años de
                 inicio de los depositos masivos 
                 de sulfuros volcanicos",
                 xlab="Años de inicio",
                 ylab="Densidad de probabilidad",
                 col="blue")

# Parámetros Log-normales

h<-length(histograma$counts)
medialog <- mean(log(año_sin_outliers)) 
sd_log<-sd(log(año_sin_outliers)) 
sd_log
## [1] 0.02064737
medialog
## [1] 7.568395
x <- seq(min(año_sin_outliers),max(año_sin_outliers),0.01) 
curve(dlnorm(x,meanlog = medialog,sdlog = sd_log),add = TRUE, col=("black"), lwd=3)

#FECUENCIAS OBSERVADAS
Fo <- histograma$counts
Fo
##  [1]   4   3   8  31  38  89  60  76 109  53   5
#FRECUENCIAS ESPERADAS

P<-c(0) 
for (i in 1:h) {P[i] <-(plnorm(histograma$breaks[i+1],medialog,sd_log)- plnorm(histograma$breaks[i],medialog,sd_log))}
Fe<-P*length(año_sin_outliers)

Fe
##  [1]  0.5573708  2.6134605  9.1853289 24.4148307 49.4992402 77.1751779
##  [7] 93.2594247 88.0052637 65.3224908 38.4037423 18.0031616

Test de aprobación

Test de Pearson 

#Tamaño muestral sin outliers
n<-length(año_sin_outliers)
n
## [1] 476
#Representar la frecuencia observada y esperada en porcentaje
Fo<-(Fo/n)*100
Fo
##  [1]  0.8403361  0.6302521  1.6806723  6.5126050  7.9831933 18.6974790
##  [7] 12.6050420 15.9663866 22.8991597 11.1344538  1.0504202
Fe<-(Fe/n)*100
Fe
##  [1]  0.1170947  0.5490463  1.9296909  5.1291661 10.3990000 16.2132727
##  [7] 19.5923161 18.4885008 13.7232124  8.0680131  3.7821768
#Correlacionar Fo y Fe
plot(Fo,Fe,main="Gráfica 3: Correlación de frecuencias 
     en el modelo log-normal del año de 
     inicio",
     xlab="Frecuencia Observada(%)",
     ylab="Frecuencia esperada(%)",
     col="blue3")

abline(lm(Fe ~ Fo), col="red",lwd=2)

Correlación<-cor(Fo,Fe)*100

Correlación
## [1] 85.27661

APRUEBA EL TEST PEARSON

Test de Chi-Cuadrado 

grados_libertad <- (length(histograma$counts)-1)
grados_libertad
## [1] 10
nivel_significancia <- 0.999   #Subir el nivel de significancia para aprobar el test de chi-cuadrado

x2<-sum((Fe-Fo)^2/Fe)
x2
## [1] 17.93621
umbral_aceptacion <- qchisq(nivel_significancia, grados_libertad)
umbral_aceptacion
## [1] 29.5883
x2<umbral_aceptacion
## [1] TRUE

APRUEBA EL TEST DE CHI-CUADRADO

TABLA DE RESUMEN

Variable<-c("Año de incio")
tabla_resumen<-data.frame(Variable,round(Correlación,2),round(x2,2),round(umbral_aceptacion,2))
colnames(tabla_resumen)<-c("Variable","Test Pearson (%)","Chi Cuadrado","Umbral de aceptación")
kable(tabla_resumen, format = "markdown", caption = "Tabla.Resumen de test de bondad al modelo de probabilidad")
Tabla.Resumen de test de bondad al modelo de probabilidad
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Año de incio 85.28 17.94 29.59

Cálculo de probabilidades

¿Cuál es la probabilidad de que un yacimiento haya iniciado entre los años 1850 y 1900?

# PROBABILIDAD ENTRE 1850 y 1900
probabilidad_Año <- plnorm(1900, meanlog = medialog, sdlog = sd_log) -
  plnorm(1850, meanlog = medialog, sdlog = sd_log)

# En porcentaje
probabilidad_Año * 100
## [1] 16.75973
# Rango para la curva
x <- seq(min(año_sin_outliers), max(año_sin_outliers), 0.01)

# Curva log-normal
plot(x, dlnorm(x, meanlog = medialog, sdlog = sd_log),
     col = "skyblue3",
     lwd = 2,
     main = "Gráfica 4. Cálculo de probabilidades del año de inicio
     de los depósitos masivos de sulfuros volcánicos",
     ylab = "Densidad de probabilidad",
     xlab = "Año de inicio")

# Área de probabilidad
x_area <- seq(1850, 1900, 0.01)
y_area <- dlnorm(x_area, meanlog = medialog, sdlog = sd_log)

lines(x_area, y_area, col = "red", lwd = 2)

polygon(c(x_area, rev(x_area)),
        c(y_area, rep(0, length(y_area))),
        col = rgb(1, 0, 0, 0.5), border = "black")

legend("topright",
       legend = c("Modelo Log-normal", "Área de Probabilidad"),
       col = c("skyblue3", "red"),
       lwd = 2,
       pch = c(NA, 15),
       cex = 0.8)

texto_prob <- paste0("Probabilidad = ",
                     round(probabilidad_Año * 100, 2), " %")

text(x = 1875,
     y = max(dlnorm(x, medialog, sd_log)) * 0.8,
     labels = texto_prob,
     font = 2)

Conclusión