Distribución de la cantidad de anomalías geoquímicas

ANÁLISIS INFERENCIAL

CARGA DE DATOS Y LIBRERÍAS

#Limpiar entorno
rm(list = ls())

#Cargar librerías
if (!require("readr")) install.packages("readr")
if (!require("dplyr")) install.packages("dplyr")
if (!require("knitr")) install.packages("knitr")
if (!require("moments")) install.packages("moments")

library(readr)
library(dplyr)
library(knitr)
library(moments)
library(gt)

#Cargar datos
ruta <- "D:/geoquimica datos.csv"
datos <- read_csv(ruta)

#====================================================================
# LIMPIEZA DE LA VARIABLE
#====================================================================

anomalias <- as.numeric(datos$geoquimic_anomaly_count)

anomalias <- na.omit(anomalias)

GRAFICA Y TABLA DE DISTRIBUCIÓN DE FRECUENCIAS

#====================================================================
# GRAFICA Y TABLA DE DISTRIBUCIÓN DE FRECUENCIAS
#====================================================================

histograma_anom <- hist(
  anomalias,
  breaks = seq(
    min(anomalias)-0.5,
    max(anomalias)+0.5,
    by = 1
  ),
  main = "Grafica Nº1: Distribución de la cantidad de anomalías geoquímicas",
  xlab = "Cantidad de anomalías geoquímicas",
  ylab = "Frecuencia",
  col = "gray"
)

#====================================================================
# TABLA DE DISTRIBUCIÓN DE FRECUENCIAS
#====================================================================

tabla_freq <- table(anomalias)

# LÍMITE INFERIOR

lis <- as.numeric(names(tabla_freq))

# LÍMITE SUPERIOR

lss <- lis

# MARCA DE CLASE

MC_prof <- lis

# FRECUENCIA ABSOLUTA

ni_prof <- as.numeric(tabla_freq)

# FRECUENCIA RELATIVA

hi_prof <- (ni_prof/sum(ni_prof))*100

# TABLA BASE

TDFprof <- round(
  data.frame(
    lis,
    lss,
    MC_prof,
    ni_prof,
    hi_prof
  ),
  2
)

# FILA TOTAL

fila_total_prof <- data.frame(
  lis = "TOTAL",
  lss = "",
  MC_prof = "",
  ni_prof = sum(ni_prof),
  hi_prof = round(sum(hi_prof),2)
)

TDFprof_total <- rbind(
  TDFprof,
  fila_total_prof
)

# TABLA FINAL

tabla_anomalias <- TDFprof_total %>%
  gt() %>%
  tab_header(
    title = md("*Tabla Nº1*"),
    subtitle = md(
      "Tabla de distribución de frecuencias de la cantidad de anomalías geoquímicas"
    )
  ) %>%
  cols_label(
    lis = "Valor",
    lss = "Valor",
    MC_prof = "Conteo",
    ni_prof = "Frecuencia absoluta",
    hi_prof = "Frecuencia relativa (%)"
  ) %>%
  tab_source_note(
    source_note = md("Autor: Grupo 2")
  ) %>%
  tab_options(
    table.border.top.color = "black",
    table.border.bottom.color = "black",
    table.border.top.style = "solid",
    table.border.bottom.style = "solid",
    column_labels.border.top.color = "black",
    column_labels.border.bottom.color = "black",
    column_labels.border.bottom.width = px(2),
    row.striping.include_table_body = TRUE,
    heading.border.bottom.color = "black",
    heading.border.bottom.width = px(2),
    table_body.hlines.color = "gray",
    table_body.border.bottom.color = "black"
  )

tabla_anomalias

Tabla Nº1
Tabla de distribución de frecuencias de la cantidad de anomalías geoquímicas
Valor	Valor	Conteo	Frecuencia absoluta	Frecuencia relativa (%)
1	1	1	1	0.04
3	3	3	3	0.12
4	4	4	12	0.48
5	5	5	28	1.12
6	6	6	57	2.28
7	7	7	112	4.48
8	8	8	171	6.84
9	9	9	213	8.52
10	10	10	285	11.40
11	11	11	285	11.40
12	12	12	292	11.68
13	13	13	262	10.48
14	14	14	232	9.28
15	15	15	159	6.36
16	16	16	143	5.72
17	17	17	89	3.56
18	18	18	68	2.72
19	19	19	34	1.36
20	20	20	27	1.08
21	21	21	11	0.44
22	22	22	6	0.24
23	23	23	5	0.20
24	24	24	2	0.08
25	25	25	2	0.08
26	26	26	1	0.04
TOTAL			2500	100.00
Autor: Grupo 2

CONJETURA MODELO

#====================================================================
# HISTOGRAMA Y MODELO POISSON
#====================================================================

prof_1 <- anomalias

Histograma_1 <- hist(
  prof_1,
  freq = FALSE,
  breaks = seq(
    min(prof_1)-0.5,
    max(prof_1)+0.5,
    by = 1
  ),
  main = "Grafica Nº2: Comparación de la realidad con el modelo de probabilidad
Poisson de la cantidad de anomalías geoquímicas",
  ylab = "Probabilidad",
  xlab = "Cantidad de anomalías geoquímicas",
  col = "lightgray",
  border = "black"
)

#====================================================================
# CALCULO DEL PARAMETRO POISSON
#====================================================================

lambda_1 <- mean(prof_1)

lambda_1

## [1] 11.98

#====================================================================
# VALORES POSIBLES
#====================================================================

x <- min(prof_1):max(prof_1)

#====================================================================
# CURVA POISSON
#====================================================================

points(
  x,
  dpois(x, lambda_1),
  type = "b",
  pch = 19,
  col = "blue"
)

TEST DE APROBACIÓN

#====================================================================
# TAMAÑO MUESTRAL
#====================================================================

n1 <- length(prof_1)

n1

## [1] 2500

#====================================================================
# FRECUENCIAS OBSERVADAS
#====================================================================

Fo_tabla <- table(prof_1)

Fo_1 <- as.numeric(Fo_tabla)

Fo_1

##  [1]   1   3  12  28  57 112 171 213 285 285 292 262 232 159 143  89  68  34  27
## [20]  11   6   5   2   2   1

#====================================================================
# VALORES OBSERVADOS
#====================================================================

Valores_1 <- as.numeric(names(Fo_tabla))

Valores_1

##  [1]  1  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

#====================================================================
# PROBABILIDADES POISSON
#====================================================================

P1 <- dpois(
  Valores_1,
  lambda_1
)

P1

##  [1] 7.509464e-05 1.796269e-03 5.379825e-03 1.289006e-02 2.573715e-02
##  [6] 4.404730e-02 6.596083e-02 8.780120e-02 1.051858e-01 1.145569e-01
## [11] 1.143660e-01 1.053927e-01 9.018602e-02 7.202857e-02 5.393139e-02
## [16] 3.800577e-02 2.529495e-02 1.594913e-02 9.553530e-03 5.450061e-03
## [21] 2.967806e-03 1.545840e-03 7.716317e-04 3.697659e-04 1.703768e-04

#====================================================================
# NUMERO DE CATEGORIAS
#====================================================================

h1 <- length(Valores_1)

h1

## [1] 25

#====================================================================
# FRECUENCIA ESPERADA
#====================================================================

Fe_1 <- P1 * n1

Fe_1

##  [1]   0.1877366   4.4906718  13.4495621  32.2251509  64.3428846 110.1182511
##  [7] 164.9020811 219.5029923 262.9645848 286.3923387 285.9150182 263.4816860
## [13] 225.4650427 180.0714141 134.8284713  95.0144169  63.2373730  39.8728278
## [19]  23.8838239  13.6251528   7.4195151   3.8645996   1.9290793   0.9244148
## [25]   0.4259419

#====================================================================
# TEST DE PEARSON
#====================================================================

# GUARDAR FRECUENCIAS ABSOLUTAS

Fo_abs <- Fo_1

Fe_abs <- Fe_1

# EXPRESAR EN PORCENTAJES PARA GRAFICA

Fo_1 <- (Fo_abs / n1) * 100

Fo_1

##  [1]  0.04  0.12  0.48  1.12  2.28  4.48  6.84  8.52 11.40 11.40 11.68 10.48
## [13]  9.28  6.36  5.72  3.56  2.72  1.36  1.08  0.44  0.24  0.20  0.08  0.08
## [25]  0.04

Fe_1 <- (Fe_abs / n1) * 100

Fe_1

##  [1]  0.007509464  0.179626873  0.537982486  1.289006036  2.573715385
##  [6]  4.404730045  6.596083242  8.780119694 10.518583393 11.455693550
## [11] 11.436600727 10.539267439  9.018601709  7.202856565  5.393138853
## [16]  3.800576674  2.529494920  1.594913113  0.955352954  0.545006114
## [21]  0.296780602  0.154583983  0.077163172  0.036976592  0.017037676

#====================================================================
# CORRELACIÓN DE FRECUENCIAS
#====================================================================

plot(
  Fo_1,
  Fe_1,
  main = "Grafica Nº3: Correlación de frecuencias observadas y esperadas
del modelo Poisson de anomalías geoquímicas",
  xlab = "Frecuencia observada (%)",
  ylab = "Frecuencia esperada (%)",
  col = "blue3",
  pch = 19
)

abline(
  a = 0,
  b = 1,
  col = "red",
  lwd = 2
)

Correlacion_1 <- cor(
  Fo_1,
  Fe_1
) * 100

Correlacion_1

## [1] 99.76214

#====================================================================
# CHI CUADRADO
#====================================================================

# ELIMINAR CATEGORIAS CON FE < 5

indices_validos <- Fe_abs >= 5

Fo_chi <- Fo_abs[indices_validos]

Fe_chi <- Fe_abs[indices_validos]

# GRADOS DE LIBERTAD

# k - 1 - parametros estimados
# en Poisson se estima lambda

grados_libertad_1 <-
length(Fo_chi) - 2

grados_libertad_1

## [1] 17

nivel_significancia <- 0.95

Chi_cuadrado <- sum(
  ((Fo_chi - Fe_chi)^2) / Fe_chi
)

Chi_cuadrado

## [1] 9.928306

umbral_aceptacion_1 <- qchisq(
  nivel_significancia,
  grados_libertad_1
)

umbral_aceptacion_1

## [1] 27.58711

#====================================================================
# TABLA RESUMEN
#====================================================================

Variable <- c(
  "Cantidad de anomalías geoquímicas"
)

tabla_resumen_1 <- data.frame(
  Variable,
  round(Correlacion_1,2),
  round(Chi_cuadrado,2),
  round(umbral_aceptacion_1,2)
)

colnames(tabla_resumen_1) <- c(
  "Variable",
  "Test Pearson (%)",
  "Chi Cuadrado",
  "Umbral de aceptación"
)

kable(
  tabla_resumen_1,
  format = "markdown",
  caption = "Tabla resumen del modelo Poisson"
)

Tabla resumen del modelo Poisson

Variable	Test Pearson (%)	Chi Cuadrado	Umbral de aceptación
Cantidad de anomalías geoquímicas	99.76	9.93	27.59

CALCULO DE PROBABILIDADES

#====================================================================
# PROBABILIDAD
#====================================================================

#¿Cuál es la probabilidad de que una muestra futura
#presente entre 10 y 15 anomalías geoquímicas?

Probabilidad_1 <-
(
  ppois(15, lambda_1) -
  ppois(9, lambda_1)
) * 100

Probabilidad_1

## [1] 60.1716

#====================================================================
# GRAFICA DE PROBABILIDAD
#====================================================================

x <- min(prof_1):max(prof_1)

y <- dpois(
  x,
  lambda_1
)

plot(
  x,
  y,
  type = "h",
  lwd = 3,
  col = "skyblue3",
  main = "Grafica Nº4: Cálculo de probabilidades",
  ylab = "Probabilidad",
  xlab = "Cantidad de anomalías geoquímicas"
)

#=========================================================
# RANGO DE INTERES
#=========================================================

x_section <- 10:15

y_section <- dpois(
  x_section,
  lambda_1
)

points(
  x_section,
  y_section,
  pch = 19,
  col = "red"
)

segments(
  x_section,
  0,
  x_section,
  y_section,
  col = "red",
  lwd = 3
)

legend(
  "topright",
  legend = c(
    "Modelo Poisson",
    "Área de Probabilidad"
  ),
  col = c(
    "skyblue3",
    "red"
  ),
  lwd = 3,
  cex = 0.7
)

texto_prob <- paste0(
  "Probabilidad = ",
  round(
    Probabilidad_1,
    2
  ),
  " %"
)

text(
  x = min(x)+3,
  y = max(y)*0.90,
  labels = texto_prob,
  col = "black",
  cex = 0.8,
  font = 2
)

#====================================================================
# CANTIDAD ESPERADA
#====================================================================

#¿De 300 futuras muestras cuántas presentarán
#entre 10 y 15 anomalías geoquímicas?

cantidad_1 <-
(
  ppois(15, lambda_1) -
  ppois(9, lambda_1)
) * 300

cantidad_1

## [1] 180.5148

INTERVALOS DE CONFIANZA

#====================================================================
# INTERVALOS DE CONFIANZA
#====================================================================

# Parametro lambda

x <- mean(prof_1)

x

## [1] 11.98

# Tamaño muestral

n <- length(prof_1)

n

## [1] 2500

# Error estándar Poisson

e <- sqrt(x/n)

e

## [1] 0.06922427

# Intervalo de confianza 95%

li <- x - 1.96*e

li

## [1] 11.84432

ls <- x + 1.96*e

ls

## [1] 12.11568

tabla_media <- data.frame(
  round(li,2),
  Variable,
  round(ls,2),
  round(e,4)
)

colnames(tabla_media) <- c(
  "Limite inferior",
  "Lambda poblacional",
  "Limite superior",
  "Error estándar"
)

kable(
  tabla_media,
  format = "markdown",
  caption = "Tabla Nro.2: Intervalo de confianza del parámetro λ"
)

Tabla Nro.2: Intervalo de confianza del parámetro λ

Limite inferior	Lambda poblacional	Limite superior	Error estándar
11.84	Cantidad de anomalías geoquímicas	12.12	0.0692

CONCLUSION

"La variable cantidad de anomalías geoquímicas presenta un comportamiento compatible con un modelo de distribución Poisson, debido a que representa datos de conteo dentro de las muestras analizadas. El parámetro λ obtenido fue de 11.98 anomalías por muestra, calculado a partir de una muestra total de 2500 registros.

El ajuste del modelo fue adecuado, ya que el Chi cuadrado calculado fue 9.93, menor que el umbral de aceptación de 27.59, y el Test de Pearson alcanzó un valor de 99.76 %, evidenciando una alta correspondencia entre las frecuencias observadas y esperadas.

Además, se calculó que la probabilidad de que una muestra futura presente entre 10 y 15 anomalías geoquímicas es de 60.17 %, por lo que, de 300 futuras muestras, aproximadamente 181 estarían dentro de ese rango. Finalmente, con un 95 % de confianza, el parámetro λ poblacional se encuentra entre 11.84 y 12.12 anomalías geoquímicas."

## [1] "La variable cantidad de anomalías geoquímicas presenta un comportamiento compatible con un modelo de distribución Poisson, debido a que representa datos de conteo dentro de las muestras analizadas. El parámetro λ obtenido fue de 11.98 anomalías por muestra, calculado a partir de una muestra total de 2500 registros.\n\nEl ajuste del modelo fue adecuado, ya que el Chi cuadrado calculado fue 9.93, menor que el umbral de aceptación de 27.59, y el Test de Pearson alcanzó un valor de 99.76 %, evidenciando una alta correspondencia entre las frecuencias observadas y esperadas.\n\nAdemás, se calculó que la probabilidad de que una muestra futura presente entre 10 y 15 anomalías geoquímicas es de 60.17 %, por lo que, de 300 futuras muestras, aproximadamente 181 estarían dentro de ese rango. Finalmente, con un 95 % de confianza, el parámetro λ poblacional se encuentra entre 11.84 y 12.12 anomalías geoquímicas."