INTRODUCCIÓN

Este trabajo tiene como objetivo analizar diferentes distribuciones de probabilidad discretas y continuas utilizando R Markdown, a partir de datos reales de la ciudad de Guadalajara de Buga.

Se emplean bases de datos de áreas como salud, educación y economía, las cuales son procesadas en R para construir tablas con kable(), gráficos y análisis estadístico.

FUNCIONES BÁSICAS PROGRAMADAS MANUALMENTE

library(knitr)
library(kableExtra)

# Suma de un vector
sumamanu <- function(vector) {
  suma <- 0
  for (i in vector) {
    suma <- suma + i
  }
  return(suma)
}

# Media de un vector
mediamanu <- function(vector) {
  n <- length(vector)
  return(sumamanu(vector) / n)
}

# Raiz cuadrada por Newton-Raphson
raizmanu <- function(numero) {
  raiz <- numero / 2
  for (i in 1:100) {
    raiz <- (raiz + numero / raiz) / 2
  }
  return(raiz)
}

# Exponencial por serie de Taylor
expomanual <- function(x) {
  resultado <- 1
  termino <- 1
  for (i in 1:150) {
    termino <- termino * x / i
    resultado <- resultado + termino
  }
  return(resultado)
}

# Logaritmo natural por serie
logmanual <- function(x) {
  y <- (x - 1) / (x + 1)
  resultado <- 0
  potencia <- y
  i <- 1
  while (i <= 199) {
    resultado <- resultado + potencia / i
    potencia <- potencia * y * y
    i <- i + 2
  }
  return(2 * resultado)
}

# Factorial
factorialmanu <- function(n) {
  if (n == 0) return(1)
  resultado <- 1
  for (i in 1:n) {
    resultado <- resultado * i
  }
  return(resultado)
}

# Desviación estándar
desviamanu <- function(vector) {
  media <- mediamanu(vector)
  n <- length(vector)
  suma <- 0
  for (i in vector) {
    dife <- i - media
    suma <- suma + dife * dife
  }
  return(raizmanu(suma / (n - 1)))
}

# Mínimo de un vector
minmanu <- function(vector) {
  m <- vector[1]
  for (i in 1:length(vector)) {
    if (vector[i] < m) {
      m <- vector[i]
    }
  }
  return(m)
}

# Máximo de un vector
maxmanu <- function(vector) {
  m <- vector[1]
  for (i in 1:length(vector)) {
    if (vector[i] > m) {
      m <- vector[i]
    }
  }
  return(m)
}

# Secuencia de n puntos entre dos valores
seqmanu <- function(desde, hasta, npuntos) {
  paso <- (hasta - desde) / (npuntos - 1)
  resultado <- c()
  for (i in 0:(npuntos - 1)) {
    resultado <- c(resultado, desde + i * paso)
  }
  return(resultado)
}

# Constante pi
numeropi <- 3.141592653589793

DISTRIBUCIÓN LOGNORMAL

Surge en procesos donde el cambio es multiplicativo y no aditivo. Se define a partir de una variable \(X\) tal que su logaritmo natural, \(\ln(X)\), sigue una distribución normal. Este modelo es esencial para variables que solo pueden tomar valores positivos y que presentan una asimetría pronunciada hacia la derecha.

La densidad se expresa como:

\[f(x) = \frac{1}{x\,\sigma\sqrt{2\pi}} \exp\!\left( -\frac{(\ln x - \mu)^2}{2\sigma^2} \right)\]

  • \(\ln x\): Es el logaritmo natural del valor observado.
  • \(\mu\) y \(\sigma\): Son el promedio y la desviación de los logaritmos de los datos.
  • Relación: El promedio real se calcula con \(e^{\mu + \frac{\sigma^2}{2}}\), por lo que sube si hay más desigualdad (\(\sigma^2\)).

A diferencia de la normal, su media siempre es mayor que su mediana. Esto la hace adecuada para modelar fenómenos con frecuentes valores bajos pero con presencia de valores extremos muy altos, como la distribución de ingresos, el precio de activos financieros o la duración de procesos biológicos.

EJEMPLO

Para el análisis de la distribución Lognormal se utilizaron datos del impuesto predial de Guadalajara de Buga. La variable seleccionada fue el avalúo catastral de los predios, debido a que representa valores positivos y presenta asimetría hacia la derecha, característica típica de una distribución lognormal.

predio <- c("PREDIO URBANO","CORREGIMIENTO 06","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO RURAL","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "CORREGIMIENTO 06","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO RURAL","PREDIO URBANO","PREDIO URBANO",
            "PREDIO URBANO","PREDIO URBANO","PREDIO URBANO")

estrato <- c("1","SIN INFORMACION","4","1","3","SIN INFORMACION","1","1",
             "SIN INFORMACION","1","3","3","2","2","2","3","3","1","1","1",
             "SIN INFORMACION","1","3","2","1","1","1","3","1","3","1","1",
             "1","SIN INFORMACION","2","1","2","2","2")

avaluo <- c(9046000,5764000,15224000,80079000,57205000,10996000,30171000,
            23126000,296831000,59354000,49835000,48282000,25099000,31636000,
            74738000,108501000,107093000,17326000,494000,31550000,590000,
            53381000,79756000,25459000,494000,494000,66862000,64780000,
            52993000,57205000,47999000,494000,52601000,29256000,62715000,
            71094000,33128000,45378000,63625000)

datoslognormal <- data.frame(
  Tipo_Predio = predio,
  Estrato    = estrato,
  Avaluo     = avaluo
)

kable(datoslognormal,
      caption = "Tabla 1. Avalúos catastrales - Impuesto predial Buga") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  scroll_box(height = "300px")
Tabla 1. Avalúos catastrales - Impuesto predial Buga
Tipo_Predio Estrato Avaluo
PREDIO URBANO 1 9046000
CORREGIMIENTO 06 SIN INFORMACION 5764000
PREDIO URBANO 4 15224000
PREDIO URBANO 1 80079000
PREDIO URBANO 3 57205000
PREDIO RURAL SIN INFORMACION 10996000
PREDIO URBANO 1 30171000
PREDIO URBANO 1 23126000
PREDIO URBANO SIN INFORMACION 296831000
PREDIO URBANO 1 59354000
PREDIO URBANO 3 49835000
PREDIO URBANO 3 48282000
PREDIO URBANO 2 25099000
PREDIO URBANO 2 31636000
PREDIO URBANO 2 74738000
PREDIO URBANO 3 108501000
PREDIO URBANO 3 107093000
PREDIO URBANO 1 17326000
PREDIO URBANO 1 494000
PREDIO URBANO 1 31550000
CORREGIMIENTO 06 SIN INFORMACION 590000
PREDIO URBANO 1 53381000
PREDIO URBANO 3 79756000
PREDIO URBANO 2 25459000
PREDIO URBANO 1 494000
PREDIO URBANO 1 494000
PREDIO URBANO 1 66862000
PREDIO URBANO 3 64780000
PREDIO URBANO 1 52993000
PREDIO URBANO 3 57205000
PREDIO URBANO 1 47999000
PREDIO URBANO 1 494000
PREDIO URBANO 1 52601000
PREDIO RURAL SIN INFORMACION 29256000
PREDIO URBANO 2 62715000
PREDIO URBANO 1 71094000
PREDIO URBANO 2 33128000
PREDIO URBANO 2 45378000
PREDIO URBANO 2 63625000

Aunque el data frame contiene variables de contexto como el tipo de predio y el estrato social, la distribución Lognormal se aplicó únicamente sobre la variable avalúo.

# Calculamos el logaritmo natural de cada avaluo
logavaluo <- c()
for (i in 1:length(avaluo)) {
  logavaluo <- c(logavaluo, logmanual(avaluo[i]))
}

# Parametros de la distribucion lognormal
mu    <- mediamanu(logavaluo)
sigma <- desviamanu(logavaluo)

tablaparam <- data.frame(
  Parametro = c("mu (media de los logaritmos)", "sigma (desv. de los logaritmos)"),
  Valor     = c(round(mu, 4), round(sigma, 4))
)

kable(tablaparam,
      caption = "Parámetros estimados - Distribución lognormal") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Parámetros estimados - Distribución lognormal
Parametro Valor
mu (media de los logaritmos) 6.5686
sigma (desv. de los logaritmos) 0.0003

Se calcularon los parámetros de la distribución Lognormal usando el logaritmo natural de los avalúos.

# Funcion de densidad lognormal programada manualmente
lognormalmanu <- function(x, mu, sigma) {
  parte1    <- 1 / (x * sigma * raizmanu(2 * numeropi))
  exponente <- -((logmanual(x) - mu)^2) / (2 * sigma^2)
  parte2    <- expomanual(exponente)
  return(parte1 * parte2)
}

Comparamos nuestra función con la de R:

valorprueba <- 9046000

comparacionln <- data.frame(
  Valor  = valorprueba,
  Manual = lognormalmanu(valorprueba, mu, sigma),
  R      = dlnorm(valorprueba, meanlog = mu, sdlog = sigma)
)

kable(comparacionln,
      caption = "Tabla 2. Comparación función manual vs R - Lognormal",
      digits = 15) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 2. Comparación función manual vs R - Lognormal
Valor Manual R
9046000 0.0001629094 0

Tenemos como resultado que las dos funciones arrojan el mismo valor.

Para la distribución Lognormal la gráfica más apropiada es el histograma de densidad con la curva teórica superpuesta, ya que permite ver visualmente la asimetría hacia la derecha que caracteriza esta distribución.

hist(datoslognormal$Avaluo,
     probability = TRUE,
     main   = "Distribución Lognormal - Avalúo catastral",
     xlab   = "Avalúo (COP)",
     col    = "#AED6F1",
     border = "white",
     breaks = 12)

xmin <- minmanu(datoslognormal$Avaluo) * 0.5
xmax <- maxmanu(datoslognormal$Avaluo) * 1.1
xseq <- seqmanu(xmin, xmax, 300)

yseq <- c()
for (i in 1:length(xseq)) {
  yseq <- c(yseq, lognormalmanu(xseq[i], mu, sigma))
}

lines(xseq, yseq, col = "#154360", lwd = 3)

CONCLUSIÓN

Los avalúos prediales presentan características de una distribución Lognormal: datos estrictamente positivos y clara asimetría hacia la derecha. La función programada manualmente presentó resultados equivalentes a la función incorporada en R.

DISTRIBUCIÓN GAUSSIANA (NORMAL)

Este modelo se fundamenta en el Teorema del Límite Central, el cual establece que la suma de variables independientes con varianza finita tiende a esta distribución a medida que el número de observaciones crece. Se caracteriza por ser la distribución de máxima incertidumbre para una media y varianza dadas.

Su función de densidad de probabilidad es:

\[f(x) = \frac{1}{\sigma\sqrt{2\pi}}\; e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}\]

  • \(x\): Es el valor de la variable aleatoria.
  • \(\mu\) (media): Determina el centro de la campana.
  • \(\sigma\) (desviación estándar): Controla la amplitud de la campana.
  • \(\sigma^2\) (varianza): Es el cuadrado de \(\sigma\). Si cambia \(\sigma\), la campana se ensancha o se encoge.
  • \(\pi\) y \(e\): Constantes que garantizan que el área total bajo la curva sea 1.

EJEMPLO

Para el análisis de la distribución Normal se utilizaron datos de nacidos vivos en Guadalajara de Buga durante 2025. La variable principal fue el peso al nacer (kg).

sexo <- c("MASCULINO","FEMENINO","MASCULINO","MASCULINO","FEMENINO",
          "MASCULINO","MASCULINO","MASCULINO","FEMENINO","MASCULINO",
          "FEMENINO","FEMENINO","MASCULINO","FEMENINO","FEMENINO",
          "MASCULINO","MASCULINO","FEMENINO","FEMENINO")

peso <- c(3.695, 2.530, 3.605, 3.275, 2.675,
          3.170, 3.255, 2.940, 2.730, 2.870,
          2.670, 2.445, 3.060, 3.620, 2.905,
          3.700, 3.795, 3.250, 3.480)

datosgaussiana <- data.frame(Sexo = sexo, PesoNac = peso)

kable(datosgaussiana,
      caption = "Tabla 3. Peso al nacer - Nacidos vivos Buga 2025") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  scroll_box(height = "250px")
Tabla 3. Peso al nacer - Nacidos vivos Buga 2025
Sexo PesoNac
MASCULINO 3.695
FEMENINO 2.530
MASCULINO 3.605
MASCULINO 3.275
FEMENINO 2.675
MASCULINO 3.170
MASCULINO 3.255
MASCULINO 2.940
FEMENINO 2.730
MASCULINO 2.870
FEMENINO 2.670
FEMENINO 2.445
MASCULINO 3.060
FEMENINO 3.620
FEMENINO 2.905
MASCULINO 3.700
MASCULINO 3.795
FEMENINO 3.250
FEMENINO 3.480

Aunque el data frame incluye la variable sexo, la distribución Normal se aplicó únicamente sobre el peso al nacer.

mediapeso <- mediamanu(datosgaussiana$PesoNac)
sigmapeso <- desviamanu(datosgaussiana$PesoNac)

tablaparamnorm <- data.frame(
  Parametro = c("mu (media del peso en kg)", "sigma (desv. est. en kg)"),
  Valor     = c(round(mediapeso, 4), round(sigmapeso, 4))
)

kable(tablaparamnorm,
      caption = "Parámetros estimados - Distribución normal") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Parámetros estimados - Distribución normal
Parametro Valor
mu (media del peso en kg) 3.1405
sigma (desv. est. en kg) 0.4277 
# Funcion de densidad normal programada manualmente
distrinormal <- function(x, media, sigma) {
  parte1    <- 1 / (sigma * raizmanu(2 * numeropi))
  dife      <- x - media
  exponente <- -((dife * dife) / (2 * sigma * sigma))
  parte2    <- expomanual(exponente)
  return(parte1 * parte2)
}

Ahora realizamos la comparación de la función de la distribución Gaussiana realizada manualmente con la función que viene incluida en RStudio.

valorpruebag <- 3.695

comparaciong <- data.frame(
  Valor  = valorpruebag,
  Manual = distrinormal(valorpruebag, mediapeso, sigmapeso),
  R      = dnorm(valorpruebag, mean = mediapeso, sd = sigmapeso)
)

kable(comparaciong,
      caption = "Tabla 4. Comparación función manual vs R - Normal",
      digits = 10) %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 4. Comparación función manual vs R - Normal
Valor Manual R
3.695 0.4025269 0.4025269

Tenemos como resultado que las dos funciones arrojan el mismo valor.

Para la distribución Normal la gráfica más apropiada es el histograma de densidad con curva de campana superpuesta, ya que permite observar visualmente la simetría alrededor de la media y la forma característica de campana de Gauss.

hist(datosgaussiana$PesoNac,
     probability = TRUE,
     main   = "Distribución Normal - Peso al nacer (kg)",
     xlab   = "Peso (kg)",
     col    = "#AED6F1",
     border = "white",
     breaks = 8)

xming <- minmanu(datosgaussiana$PesoNac) - 0.5
xmaxg <- maxmanu(datosgaussiana$PesoNac) + 0.5
xseqg <- seqmanu(xming, xmaxg, 300)

yseqg <- c()
for (i in 1:length(xseqg)) {
  yseqg <- c(yseqg, distrinormal(xseqg[i], mediapeso, sigmapeso))
}

lines(xseqg, yseqg, col = "#154360", lwd = 3)

CONCLUSIÓN

El peso al nacer presentó un comportamiento cercano a una distribución Normal, concentrándose alrededor de la media calculada. La función programada manualmente presentó resultados equivalentes a la función de R.

DISTRIBUCIÓN CHI-CUADRADO

Es una distribución de muestreo que se obtiene al sumar los cuadrados de \(k\) variables aleatorias normales estándar independientes. Su comportamiento está regido por el parámetro \(k\), denominado grados de libertad.

A medida que \(k\) aumenta, la distribución gana simetría y se aproxima a una normal. Su utilidad principal reside en las pruebas de independencia, donde se mide cuánto se desvían los datos observados de los esperados bajo una hipótesis.

\[\chi^2 = \sum \frac{(O - E)^2}{E}\]

  • \(O\): Frecuencia observada en cada celda.
  • \(E\): Frecuencia esperada bajo independencia.
  • Grados de libertad: \((filas - 1) \times (columnas - 1)\).

EJEMPLO

Para la distribución Chi-cuadrado se utilizaron datos de instituciones educativas públicas de Guadalajara de Buga. Se analizó la relación entre zona geográfica (Urbana/Rural) y nivel educativo ofrecido.

zona <- c("Urbana","Urbana","Urbana","Urbana",
          "Rural","Rural","Rural","Rural","Rural","Rural",
          "Rural","Rural","Rural","Rural","Rural","Rural",
          "Rural","Rural","Rural","Rural")

nivel <- c("Media","Media","Secundaria","Secundaria",
           "Media","Primaria","Primaria","Primaria","Primaria","Completo",
           "Primaria","Secundaria","Primaria","Completo","Secundaria",
           "Media","Secundaria","Secundaria","Primaria","Completo")

datoschi <- data.frame(Zona = zona, Nivel = nivel)

kable(datoschi,
      caption = "Tabla 5. Instituciones educativas - Zona y nivel") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  scroll_box(height = "250px")
Tabla 5. Instituciones educativas - Zona y nivel
Zona Nivel
Urbana Media
Urbana Media
Urbana Secundaria
Urbana Secundaria
Rural Media
Rural Primaria
Rural Primaria
Rural Primaria
Rural Primaria
Rural Completo
Rural Primaria
Rural Secundaria
Rural Primaria
Rural Completo
Rural Secundaria
Rural Media
Rural Secundaria
Rural Secundaria
Rural Primaria
Rural Completo

Contamos las frecuencias de cada combinación zona-nivel de forma manual, recorriendo el data frame fila por fila:

ntotal <- length(datoschi$Zona)

# Conteos por combinacion
urbmedia      <- 0
urbsecundaria <- 0
urbprimaria   <- 0
urbcompleto   <- 0

ruralmedia      <- 0
ruralsecundaria <- 0
ruralprimaria   <- 0
ruralcompleto   <- 0

for (i in 1:ntotal) {
  if (datoschi$Zona[i] == "Urbana") {
    if (datoschi$Nivel[i] == "Media")      { urbmedia      <- urbmedia + 1 }
    if (datoschi$Nivel[i] == "Secundaria") { urbsecundaria <- urbsecundaria + 1 }
    if (datoschi$Nivel[i] == "Primaria")   { urbprimaria   <- urbprimaria + 1 }
    if (datoschi$Nivel[i] == "Completo")   { urbcompleto   <- urbcompleto + 1 }
  }
  if (datoschi$Zona[i] == "Rural") {
    if (datoschi$Nivel[i] == "Media")      { ruralmedia      <- ruralmedia + 1 }
    if (datoschi$Nivel[i] == "Secundaria") { ruralsecundaria <- ruralsecundaria + 1 }
    if (datoschi$Nivel[i] == "Primaria")   { ruralprimaria   <- ruralprimaria + 1 }
    if (datoschi$Nivel[i] == "Completo")   { ruralcompleto   <- ruralcompleto + 1 }
  }
}

# Tabla de frecuencias observadas
obs <- data.frame(
  Zona       = c("Urbana", "Rural"),
  Media      = c(urbmedia,      ruralmedia),
  Secundaria = c(urbsecundaria, ruralsecundaria),
  Primaria   = c(urbprimaria,   ruralprimaria),
  Completo   = c(urbcompleto,   ruralcompleto)
)

kable(obs,
      caption = "Tabla 6. Frecuencias observadas") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 6. Frecuencias observadas
Zona Media Secundaria Primaria Completo
Urbana 2 2 0 0
Rural 2 4 7 3 
# Totales por fila
totalurbana <- urbmedia + urbsecundaria + urbprimaria + urbcompleto
totalrural  <- ruralmedia + ruralsecundaria + ruralprimaria + ruralcompleto

# Totales por columna
totalmedia      <- urbmedia      + ruralmedia
totalsecundaria <- urbsecundaria + ruralsecundaria
totalprimaria   <- urbprimaria   + ruralprimaria
totalcompleto   <- urbcompleto   + ruralcompleto

# Frecuencias esperadas
eurmedia      <- (totalurbana * totalmedia)      / ntotal
eursecundaria <- (totalurbana * totalsecundaria) / ntotal
eurprimaria   <- (totalurbana * totalprimaria)   / ntotal
eurcompleto   <- (totalurbana * totalcompleto)   / ntotal

erumedia      <- (totalrural * totalmedia)      / ntotal
erusecundaria <- (totalrural * totalsecundaria) / ntotal
eruprimaria   <- (totalrural * totalprimaria)   / ntotal
erucompleto   <- (totalrural * totalcompleto)   / ntotal

esp <- data.frame(
  Zona       = c("Urbana", "Rural"),
  Media      = c(round(eurmedia,      3), round(erumedia,      3)),
  Secundaria = c(round(eursecundaria, 3), round(erusecundaria, 3)),
  Primaria   = c(round(eurprimaria,   3), round(eruprimaria,   3)),
  Completo   = c(round(eurcompleto,   3), round(erucompleto,   3))
)

kable(esp,
      caption = "Tabla 7. Frecuencias esperadas") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 7. Frecuencias esperadas
Zona Media Secundaria Primaria Completo
Urbana 0.8 1.2 1.4 0.6
Rural 3.2 4.8 5.6 2.4 
# Estadistico chi-cuadrado: suma de (O-E)^2 / E para cada celda
chimanual <- 0
chimanual <- chimanual + (urbmedia      - eurmedia)^2      / eurmedia
chimanual <- chimanual + (urbsecundaria - eursecundaria)^2 / eursecundaria
chimanual <- chimanual + (urbprimaria   - eurprimaria)^2   / eurprimaria
chimanual <- chimanual + (urbcompleto   - eurcompleto)^2   / eurcompleto
chimanual <- chimanual + (ruralmedia      - erumedia)^2      / erumedia
chimanual <- chimanual + (ruralsecundaria - erusecundaria)^2 / erusecundaria
chimanual <- chimanual + (ruralprimaria   - eruprimaria)^2   / eruprimaria
chimanual <- chimanual + (ruralcompleto   - erucompleto)^2   / erucompleto

# Grados de libertad: (2 filas - 1) * (4 columnas - 1) = 3
gl <- 3

resultadochi <- data.frame(
  Estadistico = c("Chi-cuadrado manual", "Grados de libertad"),
  Valor       = c(round(chimanual, 6), gl)
)

kable(resultadochi,
      caption = "Tabla 8. Resultado prueba Chi-cuadrado") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 8. Resultado prueba Chi-cuadrado
Estadistico Valor
Chi-cuadrado manual 5.416667
Grados de libertad 3.000000

Ahora comparamos la función elaborada manualmente con la función propia de RStudio.

chiR <- chisq.test(table(datoschi$Zona, datoschi$Nivel),
                   simulate.p.value = FALSE)$statistic
## Warning in chisq.test(table(datoschi$Zona, datoschi$Nivel), simulate.p.value =
## FALSE): Chi-squared approximation may be incorrect
comparacionchi <- data.frame(
  Manual = round(chimanual, 6),
  R      = round(chiR, 6)
)

kable(comparacionchi,
      caption = "Tabla 9. Comparación Chi-cuadrado manual vs R") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 9. Comparación Chi-cuadrado manual vs R
Manual R
X-squared 5.416667 5.416667

R presenta una advertencia por el tamaño reducido de la muestra; sin embargo, ambos métodos arrojan el mismo estadístico.

Para la distribución Chi-cuadrado la gráfica más apropiada es el gráfico de barras agrupadas, ya que permite comparar visualmente las frecuencias de dos variables categóricas (zona y nivel educativo) una al lado de la otra.

vurbana <- c(urbmedia, urbsecundaria, urbprimaria, urbcompleto)
vrural  <- c(ruralmedia, ruralsecundaria, ruralprimaria, ruralcompleto)

barplot(rbind(vurbana, vrural),
        beside    = TRUE,
        names.arg = c("Media", "Secundaria", "Primaria", "Completo"),
        main      = "Zona geográfica vs Nivel educativo",
        xlab      = "Nivel educativo",
        ylab      = "Frecuencia",
        col       = c("#AED6F1", "#154360"),
        legend.text = c("Urbana", "Rural"),
        args.legend = list(x = "topright"))

CONCLUSIÓN

La prueba Chi-cuadrado permitió analizar la relación entre zona geográfica y niveles educativos ofrecidos. Los resultados obtenidos manualmente fueron equivalentes a los calculados por R.

DISTRIBUCIÓN DE POISSON

Este modelo describe el número de veces que ocurre un evento en un intervalo determinado de tiempo o espacio. Se deriva como el límite de una distribución Binomial cuando el número de ensayos es muy elevado y la probabilidad de éxito es muy pequeña, manteniendo una tasa promedio constante \(\lambda\).

Su función de masa de probabilidad es:

\[P(X = k) = \frac{e^{-\lambda}\, \lambda^k}{k!}\]

  • \(k\): Número de eventos que queremos calcular (0, 1, 2, 3…).
  • \(\lambda\): Tasa promedio de ocurrencia del evento.
  • \(k!\): Factorial — penaliza probabilidades para eventos muy numerosos.
  • **Dato clave*: En Poisson, la media y la varianza son iguales a \(\lambda\).

EJEMPLO

Para la distribución de Poisson se utilizaron datos de accidentalidad vial de Guadalajara de Buga. La variable de análisis es el número de accidentes por día durante el mes de enero.

accidentesdia <- c(1, 0, 2, 3, 3, 1, 2,
                   3, 0, 1, 2, 2, 1, 0,
                   1, 0, 1, 2, 2, 0, 2,
                   0, 2, 0, 0, 2, 0, 0,
                   0, 1, 0)
datospoisson <- data.frame(
  Dia           = 1:31,
  AccidentesDia = accidentesdia
)

kable(datospoisson,
      caption = "Tabla 10. Accidentes viales diarios - Buga, enero") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  scroll_box(height = "250px")
Tabla 10. Accidentes viales diarios - Buga, enero
Dia AccidentesDia
1 1
2 0
3 2
4 3
5 3
6 1
7 2
8 3
9 0
10 1
11 2
12 2
13 1
14 0
15 1
16 0
17 1
18 2
19 2
20 0
21 2
22 0
23 2
24 0
25 0
26 2
27 0
28 0
29 0
30 1
31 0 
lam <- mediamanu(datospoisson$AccidentesDia)

tablaparamP <- data.frame(
  Parametro = "lambda (tasa media diaria de accidentes)",
  Valor     = round(lam, 4)
)

kable(tablaparamP,
      caption = "Parámetro estimado - Distribución Poisson") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Parámetro estimado - Distribución Poisson
Parametro Valor
lambda (tasa media diaria de accidentes) 1.0968 
# Funcion de masa de probabilidad Poisson programada manualmente
poissonmanu <- function(k, lambda) {
  eneglambda <- expomanual(-lambda)
  lambdak    <- 1
  if (k > 0) {
    for (i in 1:k) {
      lambdak <- lambdak * lambda
    }
  }
  kfact <- factorialmanu(k)
  return(eneglambda * lambdak / kfact)
}

Calculamos la probabilidad para cada valor posible de \(k\) (de 0 a 6 accidentes):

kvals    <- 0:6
probmanu <- c()
probr    <- c()

for (k in kvals) {
  probmanu <- c(probmanu, poissonmanu(k, lam))
  probr    <- c(probr,    dpois(k, lambda = lam))
}

tablapoisson <- data.frame(
  k      = kvals,
  Manual = round(probmanu, 6),
  R      = round(probr,    6)
)

kable(tablapoisson,
      caption = "Tabla 11. Comparación P(X=k) manual vs R - Poisson") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 11. Comparación P(X=k) manual vs R - Poisson
k Manual R
0 0.333947 0.333947
1 0.366264 0.366264
2 0.200854 0.200854
3 0.073431 0.073431
4 0.020134 0.020134
5 0.004417 0.004417
6 0.000807 0.000807

Tenemos como resultado que las dos funciones arrojan los mismos valores.

Para la distribución de Poisson la gráfica más apropiada es el gráfico de barras comparando frecuencias observadas vs teóricas, ya que Poisson es una distribución discreta (cuenta eventos enteros) y las barras representan mejor cada valor posible de \(k\).

frecobs <- c()
for (k in kvals) {
  conteo <- 0
  for (i in 1:length(datospoisson$AccidentesDia)) {
    if (datospoisson$AccidentesDia[i] == k) {
      conteo <- conteo + 1
    }
  }
  frecobs <- c(frecobs, conteo)
}

barplot(rbind(frecobs, probmanu * length(datospoisson$AccidentesDia)),
        beside      = TRUE,
        names.arg   = kvals,
        main        = "Distribución Poisson - Accidentes diarios en Buga",
        xlab        = "Número de accidentes (k)",
        ylab        = "Frecuencia",
        col         = c("#AED6F1", "#154360"),
        legend.text = c("Observado", "Teórico Poisson"),
        args.legend = list(x = "topright"))

CONCLUSIÓN

El número de accidentes viales diarios en Buga durante enero sigue una distribución de Poisson. La función programada manualmente presentó resultados equivalentes a la función de R.

DISTRIBUCIÓN EXPONENCIAL

Esta distribución es la contraparte continua de la de Poisson y modela el tiempo de espera entre dos eventos sucesivos. Su característica analítica más importante es la falta de memoria, lo cual implica que la probabilidad de que ocurra un evento en el futuro es totalmente independiente del tiempo que ya haya transcurrido.

Su función de densidad es:

\[f(x) = \lambda\, e^{-\lambda x}, \quad x \geq 0\]

  • \(x\): Tiempo de espera entre eventos.
  • \(\lambda\): Tasa de ocurrencia (eventos por unidad de tiempo).
  • \(e^{-\lambda x}\): A medida que aumenta \(x\), la probabilidad de seguir esperando disminuye.
  • Falta de memoria: \(P(X > s + t \mid X > s) = P(X > t)\).

EJEMPLO

Para la distribución Exponencial se utilizaron datos de defunciones registradas en Guadalajara de Buga durante enero de 2025. La variable de análisis es el tiempo transcurrido entre fallecimientos consecutivos (en horas).

tiemposentre <- c(
  1.17, 0.17, 29.33, 0.58, 1.75, 5.08, 2.83,
  27.15, 5.85, 0.42, 0.25, 17.83, 0.63, 1.9,
  0.48, 6.57, 4.3, 2.45, 2.92, 2.92, 6.42,
  8.28, 1.22, 4.67, 0.92, 2.58, 2.17, 6.67,
  5.5, 4.17, 13.5, 6.33, 4.0, 2.95, 1.28,
  1.6, 1.17, 11.75, 6.25, 10.0, 2.58, 2.17,
  4.15, 0.1, 27.83, 6.17, 5.0, 14.77, 0.4,
  8.83, 0.5, 3.75, 0.5, 4.92, 0.83, 3.05,
  18.78, 1.63, 1.53, 0.33, 5.12, 7.1, 15.82,
  3.13, 0.15, 6.35, 5.25, 5.25, 1.58, 0.03,
  0.5, 4.82, 1.4, 8.68, 8.65, 1.17, 3.42,
  5.33, 0.42, 3.53, 7.18, 0.53, 0.18, 0.57,
  0.25, 2.83, 2.25, 3.5, 0.43, 11.38, 10.35,
  1.67, 6.17, 1.42, 14.17, 1.33, 15.25, 3.9,
  13.93, 8.37, 33.8, 5.5, 5.0, 0.3, 12.53,
  3.63, 2.93, 6.18, 9.77, 1.92, 3.9, 7.58,
  3.4, 8.68, 0.5, 12.92, 4.0, 2.45, 5.3,
  17.33, 14.03, 0.88, 4.25, 0.5, 6.63, 9.37,
  3.5, 2.5, 0.87, 5.3, 13.0, 6.33, 0.5
)

datosintervalos <- data.frame(
  Intervalo = 1:length(tiemposentre),
  TiempoHoras = tiemposentre
)

kable(datosintervalos,
      caption = "Tabla 12. Tiempo entre fallecimientos consecutivos - Buga, enero 2025 (horas)",
      col.names = c("Intervalo N", "Tiempo (horas)")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, font_size = 13) %>%
  scroll_box(height = "300px")
Tabla 12. Tiempo entre fallecimientos consecutivos - Buga, enero 2025 (horas)
Intervalo N Tiempo (horas)
1 1.17
2 0.17
3 29.33
4 0.58
5 1.75
6 5.08
7 2.83
8 27.15
9 5.85
10 0.42
11 0.25
12 17.83
13 0.63
14 1.90
15 0.48
16 6.57
17 4.30
18 2.45
19 2.92
20 2.92
21 6.42
22 8.28
23 1.22
24 4.67
25 0.92
26 2.58
27 2.17
28 6.67
29 5.50
30 4.17
31 13.50
32 6.33
33 4.00
34 2.95
35 1.28
36 1.60
37 1.17
38 11.75
39 6.25
40 10.00
41 2.58
42 2.17
43 4.15
44 0.10
45 27.83
46 6.17
47 5.00
48 14.77
49 0.40
50 8.83
51 0.50
52 3.75
53 0.50
54 4.92
55 0.83
56 3.05
57 18.78
58 1.63
59 1.53
60 0.33
61 5.12
62 7.10
63 15.82
64 3.13
65 0.15
66 6.35
67 5.25
68 5.25
69 1.58
70 0.03
71 0.50
72 4.82
73 1.40
74 8.68
75 8.65
76 1.17
77 3.42
78 5.33
79 0.42
80 3.53
81 7.18
82 0.53
83 0.18
84 0.57
85 0.25
86 2.83
87 2.25
88 3.50
89 0.43
90 11.38
91 10.35
92 1.67
93 6.17
94 1.42
95 14.17
96 1.33
97 15.25
98 3.90
99 13.93
100 8.37
101 33.80
102 5.50
103 5.00
104 0.30
105 12.53
106 3.63
107 2.93
108 6.18
109 9.77
110 1.92
111 3.90
112 7.58
113 3.40
114 8.68
115 0.50
116 12.92
117 4.00
118 2.45
119 5.30
120 17.33
121 14.03
122 0.88
123 4.25
124 0.50
125 6.63
126 9.37
127 3.50
128 2.50
129 0.87
130 5.30
131 13.00
132 6.33
133 0.50 
mediatiempo <- mediamanu(datosintervalos$TiempoHoras)
lambdae <- 1 / mediatiempo

print(mediatiempo)
## [1] 5.507293
print(lambdae)
## [1] 0.1815774
# f(x) = lambda * e^(-lambda * x)
exponencialmanu <- function(x, lambda) {
  return(lambda * expomanual(-lambda * x))
}

Comparamos nuestra función con la de R para varios valores de \(x\):

xprueba <- c( 1.17, 0.17, 29.33, 0.58, 1.75)

manuvals <- c()
rvals <- c()
for (xi in xprueba) {
  manuvals <- c(manuvals, exponencialmanu(xi, lambdae))
  rvals <- c(rvals, dexp(xi, rate = lambdae))
}

comparacionexp <- data.frame(
  x = xprueba,
  fmanual = round(manuvals, 8),
  fR = round(rvals, 8)
)

kable(comparacionexp,
      caption = "Tabla 13. Comparación f(x) manual vs R - Exponencial") %>%
  kable_styling(bootstrap_options = c("striped", "hover"),
                full_width = FALSE, font_size = 13)
Tabla 13. Comparación f(x) manual vs R - Exponencial
x fmanual fR
1.17 0.1468243 0.1468243
0.17 0.1760581 0.1760581
29.33 0.0008834 0.0008834
0.58 0.1634271 0.1634271
1.75 0.1321479 0.1321479

Tenemos como resultado que las dos distribuciones arrojan el mismo resultado.

Para la distribución exponencial la gráfica más apropiada es el histograma de densidad con curva de decaimiento superpuesta, ya que permite observar la caída rápida al inicio (muchos tiempos de espera cortos) y la cola hacia la derecha (pocos tiempos muy largos).

hist(datosintervalos$TiempoHoras,
     probability = TRUE,
     main = "Distribución Exponencial - Tiempo entre defunciones (horas)",
     xlab = "Tiempo (horas)",
     col = "#AED6F1",
     border = "white",
     breaks = 8)

inicioe <- 0
fine <- maxmanu(datosintervalos$TiempoHoras) * 1.1
xseqe <- seqmanu(inicioe, fine, 300)
yseqe <- c()
for (xi in xseqe) {
  yseqe <- c(yseqe, exponencialmanu(xi, lambdae))
}
lines(xseqe, yseqe, col = "#154360", lwd = 3)

CONCLUSIÓN

La función programada manualmente presentó resultados equivalentes a la función de R.

CONCLUSIÓN GENERAL

A lo largo de este trabajo se analizaron cinco distribuciones de probabilidad :lognormal, normal, Chi-cuadrado, Poisson y exponencial, aplicadas a datos reales de Guadalajara de Buga. En todos los casos:

  1. Se implementaron manualmente las funciones matemáticas sin usar las funciones integradas de R.
  2. Los resultados manuales fueron equivalentes a los obtenidos con las funciones nativas de R.
  3. Se utilizó kableExtra para presentar los datos de forma organizada y visualmente clara.

REFERENCIAS