Capitulo
María Guzmán Martínez
2025-12-04
# Stan
library(cmdstanr) # cmdstan_model
library(posterior) # mcmc_hist, as_draws_df,
library(bayesplot) # mcmc_trace
#-- Gráficos
library(ggpmisc) # stat_poly_line
library(ggplot2) # ggtitle
library(patchwork)
library(GGally)
library(nortest) # muestra las funciones del paquete
library(coda)
library(tinytex)
library(flextable) # flextable
library(tidyverse) # %>%
library(dplyr)
library(parameters) # model_parametersBase de datos
Modelo de regresión lineal múltiple
Registro de 50 clientes que hicieron sus compras en línea
PersonID: ID del cliente
Sex: Un valor binario que representa el género del cliente (0: mujer, 1: hombre)
Income: Ingreso anual,
Y: La proporción de compras del año
d <-read.csv(file="C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Facultad_Matematicas/1. Programa_Licenciatura en Matematicas/Cursos de la LM/Modelacion_bayesiana_LM/Libros del Curso/Codigo_libro_Bayesian/chap05/input/data-shopping-1.csv") # Datos
d<-d %>%
mutate(Income.100=Income/100,
Sex=factor(Sex)) %>%
mutate(Sex=fct_collapse(Sex, "Mujer"="0","Hombre"="1"))%>%
glimpse() # Muestra las variables## Rows: 50
## Columns: 5
## $ PersonID <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ Sex <fct> Mujer, Hombre, Mujer, Hombre, Hombre, Mujer, Hombre, Hombre…
## $ Income <dbl> 85.1, 64.6, 97.6, 72.8, 51.1, 86.9, 86.3, 84.6, 80.7, 76.6,…
## $ Y <dbl> 0.56, 0.59, 0.56, 0.65, 0.55, 0.55, 0.67, 0.63, 0.69, 0.46,…
## $ Income.100 <dbl> 0.851, 0.646, 0.976, 0.728, 0.511, 0.869, 0.863, 0.846, 0.8…N <- nrow(d) Nota: La proporción de compras se calcula como la relación entre el número de veces que un cliente ha comprado en la tienda y el número total de visitas de dicho cliente en ese año. Aunque el número total de visitas varía según cada cliente, suponemos que se sitúa entre 100 y 200.
Usar el paquete forcast
Análisis descriptivo
tabla_estad <- d[,-c(1:2)] %>%
summarise(across(everything(),
list( Min = ~ min(.x, na.rm = TRUE),
Max = ~ max(.x, na.rm = TRUE),
Media = ~ mean(.x, na.rm = TRUE),
SD = ~ sd(.x, na.rm = TRUE),
CV = ~ (sd(.x, na.rm = TRUE) / mean(.x, na.rm = TRUE)) * 100))) %>%
tidyr::pivot_longer( cols = everything(), # Formato largo
names_to = c("Variable", "Estadistico"),
names_sep = "_") %>%
tidyr::pivot_wider(names_from = Estadistico, # Formato ancho
values_from = value)
tabla_estadEstadísticas de los resultados
flextable(tabla_estad) %>%
set_header_labels( # Cambio de nombre de las variables
Variable = "Variable",
Min = "Mínimo",
Max = "Máximo",
Media = "Media",
SD = "Desviación Est.",
CV = "C.V. (%)") %>%
colformat_double(digits = 3) %>%
theme_booktabs() %>%
autofit() %>%
set_caption("Estadísticas descriptivas de las variables del estudio.")Variable | Mínimo | Máximo | Media | Desviación Est. | C.V. (%) |
|---|---|---|---|---|---|
Income | 38.700 | 105.500 | 66.786 | 17.608 | 26.365 |
Y | 0.300 | 0.770 | 0.509 | 0.111 | 21.800 |
Income.100 | 0.387 | 1.055 | 0.668 | 0.176 | 26.365 |
ggpairs(d[,-1]) 
ggpairs(d[,-1],
aes(color = Sex, alpha = 0.5),
lower = list(combo = "count"))
Modelación clásica
mod_fre<-lm(Y ~Sex+Income.100,d)
summary(mod_fre)##
## Call:
## lm(formula = Y ~ Sex + Income.100, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.081838 -0.023525 -0.000055 0.021086 0.089602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.14821 0.02128 6.964 9.27e-09 ***
## SexHombre 0.13800 0.01047 13.177 < 2e-16 ***
## Income.100 0.43750 0.03004 14.563 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03702 on 47 degrees of freedom
## Multiple R-squared: 0.8934, Adjusted R-squared: 0.8889
## F-statistic: 197 on 2 and 47 DF, p-value: < 2.2e-16Intervalos de confianza de los coeficientes de regresión:
confint(mod_fre)## 2.5 % 97.5 %
## (Intercept) 0.1053985 0.1910282
## SexHombre 0.1169288 0.1590662
## Income.100 0.3770643 0.4979332Estimación de sigma.
sigma(mod_fre) # SD## [1] 0.03701899A continuación se presentan las estimaciones de los coeficientes de regresión, así como sus intervalos de confianza; los cuales podemos observar que todos resultan ser estadísticamente significativos a un nivel de confianza del \(\alpha=0.05\).
model_parameters(mod_fre)str(d)## 'data.frame': 50 obs. of 5 variables:
## $ PersonID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sex : Factor w/ 2 levels "Mujer","Hombre": 1 2 1 2 2 1 2 2 2 1 ...
## $ Income : num 85.1 64.6 97.6 72.8 51.1 86.9 86.3 84.6 80.7 76.6 ...
## $ Y : num 0.56 0.59 0.56 0.65 0.55 0.55 0.67 0.63 0.69 0.46 ...
## $ Income.100: num 0.851 0.646 0.976 0.728 0.511 0.869 0.863 0.846 0.807 0.766 ...Modelación bayesiana
Modelo 5-3
Modelo estadístico \[ y=b_1+b_2 Sex+b_3Income+e \] donde:
\(b_0\), \(b_1\) y \(b_2\) son los coeficientes de regresión.
\(e\) son iid con \(N(0,\sigma^2)\).
Así \[ Y[n] = b_1 + b_2 \text{Sex}[n] + b_3 \text{Income}[n] + e[n] \]
\(n = 1, \dots, N\), \(e[n] \sim \mathcal{N}(0, \sigma^2)\; n = 1, \dots, N\)
El cual se puede escribir como: \[ Y[n] \sim \mathcal{N}(b_1 + b_2 \text{Sex}[n] + b_3 \text{Income}[n], \sigma), \quad n = 1, \dots, N \] Es decir, \[ \mu[n] = b_1 + b_2\text{Sex}[n] + b_3\text{Income}[n],\; n = 1, \dots, N \] \[ Y[n] \sim \mathcal{N}(\mu[n], \sigma),\; n = 1, \dots, N \] Declaración del modelo en Stan:
Modelo estimación
d <-read.csv(file="C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Facultad_Matematicas/1. Programa_Licenciatura en Matematicas/Cursos de la LM/Modelacion_bayesiana_LM/Libros del Curso/Codigo_libro_Bayesian/chap05/input/data-shopping-1.csv") # Datos
d$Income <-d$Income/100
N <-nrow(d)
data<-list(N=N,
Sex=d$Sex,
Income=d$Income,
Y=d$Y)
str(d)## 'data.frame': 50 obs. of 4 variables:
## $ PersonID: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Sex : int 0 1 0 1 1 0 1 1 1 0 ...
## $ Income : num 0.851 0.646 0.976 0.728 0.511 0.869 0.863 0.846 0.807 0.766 ...
## $ Y : num 0.56 0.59 0.56 0.65 0.55 0.55 0.67 0.63 0.69 0.46 ...model <-cmdstan_model(stan_file="C:/Documetos MGM_OS_Dell/Documentos Maria_Dell/Facultad_Matematicas/1. Programa_Licenciatura en Matematicas/Cursos de la LM/Modelacion_bayesiana_LM/Libros del Curso/Codigo_libro_Bayesian/chap05/model/model5-3.stan") fit<-model$sample(data=data,
seed=123,
parallel_chains=4,
refresh = 0 # APAGA EL PROGRESO EN PANTALLA
)## Running MCMC with 4 parallel chains...## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.8 seconds.estimaciones_d.f<-fit$summary()
estimaciones_d.f[2:5,] |>
flextable()|>
colformat_double(digits = 3) |>
set_caption("Estimaciones de los parámetros del modelo ")variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
b[1] | 0.148 | 0.148 | 0.022 | 0.022 | 0.112 | 0.185 | 1.001 | 1,776.047 | 1,937.880 |
b[2] | 0.138 | 0.138 | 0.011 | 0.010 | 0.120 | 0.155 | 1.004 | 2,166.281 | 1,768.315 |
b[3] | 0.438 | 0.438 | 0.032 | 0.032 | 0.385 | 0.490 | 1.000 | 1,820.358 | 1,784.544 |
sigma | 0.038 | 0.038 | 0.004 | 0.004 | 0.032 | 0.045 | 1.000 | 2,406.273 | 2,268.216 |
Histogramas de las cadenas de los parámetros estimados.
library(bayesplot)
muestras_all_df<-fit$draws(format = "df") # Muestras MCMC
p1 <- mcmc_hist(muestras_all_df, pars = "b[1]")
p2 <- mcmc_hist(muestras_all_df, pars = "b[2]")
p3 <- mcmc_hist(muestras_all_df, pars = "b[3]")
p4 <- mcmc_hist(muestras_all_df, pars = "sigma")
p1 / p2 / p3 / p4
Gráficos de las cadenas
t1 <- mcmc_trace(muestras_all_df, pars = "b[1]") +
ggtitle("Cadena de 'a'")
t2 <- mcmc_trace(muestras_all_df, pars = "b[1]") +
ggtitle("Cadena de 'b'")
t3 <- mcmc_trace(muestras_all_df, pars = "b[1]") +
ggtitle("Cadena de 'b'")
t4 <- mcmc_trace(muestras_all_df, pars = "sigma") +
ggtitle("Cadena de 'sigma'")
t1 / t2 / t3 / t4
muestras_df <- fit$draws( variables = c("b", "sigma", "mu[1]", "mu[50]"),
format = "df")
parametros_interes <- muestras_df[, c("b[1]", "b[2]", "b[3]", "sigma", "mu[1]", "mu[50]")]
ggpairs(
parametros_interes,
upper = list(continuous = wrap("cor", method = "spearman", color = "black")),
diag = list(continuous = wrap("barDiag", fill = "white", color = "black")),
lower = list(continuous = wrap("points", alpha = 0.3, size = 0.6))
) +
theme_bw() 
Modelo predicho
set.seed(123)
Y_sim_mean <-0.2
Y_sim <-rnorm(N,
mean=Y_sim_mean,
sd=0.1)
# Base de datos
data_sim <-list(N=N,
Sex=d$Sex,
Income=d$Income,
Y=Y_sim)
fit_sim <-model$sample(data=data_sim,
seed=123,
parallel_chains=4) ## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
## Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
## Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
## Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
## Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
## Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
## Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
## Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
## Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
## Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
## Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.