Agrupamiento Multivariado
Juan Sosa PhD
Webpage https://sites.google.com/view/juansosa/
YouTube https://www.youtube.com/c/JuanSosa1702
GitHub https://github.com/jstats1702
Rpubs https://rpubs.com/jstats1702
1 Modelos
Considere los siguientes modelos para realizar tareas de agrupamiento de observaciones multivariadas \(\boldsymbol{y}_1,\ldots,\boldsymbol{y}_n\) con \(\boldsymbol{y}_i = (y_{i,1},\ldots,y_{i,p})\) para \(i=1,\ldots,n\) y \(p\geq 2\).
1.1 Modelo 1: \(\textsf{M}_1\)
Distribución muestral: \[ \boldsymbol{y}_i\mid\xi_i,\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma}\stackrel{\text {ind}}{\sim}\small{\textsf{N}}(\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma})\,,\qquad i=1,\ldots,n\,, \] donde \(\xi_i\) es una variable categórica que asume valores en el conjunto de enteros positivos \(\{1,\ldots,H\}\) con pro-ba-bi-li-da-des \(\textsf{Pr}(\xi_i = h\mid\omega_h) = \omega_h\) para \(i=1,\ldots,n\) y \(h=1,\ldots,H\), tales que \(0<\omega_h<1\) y \(\sum_{h=1}^H \omega_h = 1\), \(\boldsymbol{\theta}_h\) es el vector de medias específico del grupo \(h\) para \(h=1,\ldots,H\), y \(\mathbf{\Sigma}\) es la matriz de covarianza común a todos los grupos.
Distribución previa: \[ p(\boldsymbol{\xi},\boldsymbol{\omega},\boldsymbol{\theta},\mathbf{\Sigma}) = p(\boldsymbol{\xi}\mid\boldsymbol{\omega})\,p(\boldsymbol{\omega})\,p(\boldsymbol{\theta})\,p(\mathbf{\Sigma}) \] donde \(\boldsymbol{\xi}=(\xi_1,\ldots,\xi_n)\), \(\boldsymbol{\omega}=(\omega_1,\ldots,\omega_H)\), y \(\boldsymbol{\theta}= (\boldsymbol{\theta}_1,\ldots,\boldsymbol{\theta}_H)\) con \[ \begin{align*} \boldsymbol{\xi}\mid \boldsymbol{\omega}&\sim \small{\textsf{Categorica}}(\boldsymbol{\omega}) \\ \boldsymbol{\omega}&\sim \small{\textsf{Dirichlet}}(\boldsymbol{\alpha}_0)\\ \boldsymbol{\theta}_h &\stackrel{\text {iid}}{\sim}\small{\textsf{N}}(\boldsymbol{\mu}_0,\mathbf{\Lambda}_0) \\ \mathbf{\Sigma}&\sim \small{\textsf{WI}}\left(\nu_0,\mathbf{S}_0^{-1}\right) \end{align*} \] y \(\boldsymbol{\alpha}_0\), \(\boldsymbol{\mu}_0\), \(\mathbf{\Lambda}_0\), \(\nu_0\), \(\mathbf{S}_0\) son hiperparámetros.
1.2 Modelo 2: \(\textsf{M}_2\)
Distribución muestral: \[ \boldsymbol{y}_i\mid\xi_i,\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma}_{\xi_i} \stackrel{\text {ind}}{\sim}\small{\textsf{N}}(\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma}_{\xi_i})\,,\qquad i=1,\ldots,n\,, \] donde \(\xi_i\) es una variable categórica que asume valores en el conjunto de enteros positivos \(\{1,\ldots,H\}\) con pro-ba-bi-li-da-des \(\textsf{Pr}(\xi_i = h\mid\omega_h) = \omega_h\) para \(i=1,\ldots,n\) y \(h=1,\ldots,H\), tales que \(0<\omega_h<1\) y \(\sum_{h=1}^H \omega_h = 1\), \(\boldsymbol{\theta}_h\) y \(\mathbf{\Sigma}_{h}\) son el vector de medias y la matriz de covarianza específicos del grupo \(h\) para \(h=1,\ldots,H\), respectivamente.
Distribución previa: \[ p(\boldsymbol{\xi},\boldsymbol{\omega},\boldsymbol{\theta},\mathbf{\Sigma}) = p(\boldsymbol{\xi}\mid\boldsymbol{\omega})\,p(\boldsymbol{\omega})\,p(\boldsymbol{\theta})\,p(\mathbf{\Sigma}) \] donde \(\boldsymbol{\xi}=(\xi_1,\ldots,\xi_n)\), \(\boldsymbol{\omega}=(\omega_1,\ldots,\omega_H)\), \(\boldsymbol{\theta}= (\boldsymbol{\theta}_1,\ldots,\boldsymbol{\theta}_H)\), y \(\mathbf{\Sigma}= (\mathbf{\Sigma}_1,\ldots,\mathbf{\Sigma}_H)\) con \[ \begin{align*} \boldsymbol{\xi}\mid \boldsymbol{\omega}&\sim \small{\textsf{Categorica}}(\boldsymbol{\omega}) \\ \boldsymbol{\omega}&\sim \small{\textsf{Dirichlet}}(\boldsymbol{\alpha}_0)\\ \boldsymbol{\theta}_h &\stackrel{\text {iid}}{\sim}\small{\textsf{N}}(\boldsymbol{\mu}_0,\mathbf{\Lambda}_0) \\ \mathbf{\Sigma}_h &\stackrel{\text {iid}}{\sim}\small{\textsf{WI}}\left(\nu_0,\mathbf{S}_0^{-1}\right) \end{align*} \] y \(\boldsymbol{\alpha}_0\), \(\boldsymbol{\mu}_0\), \(\mathbf{\Lambda}_0\), \(\nu_0\), \(\mathbf{S}_0\) son hiperparámetros.
2 Datos
Las bases de datos \(\textsf{D}_1,\textsf{D}_2,\textsf{D}_3,\textsf{D}_4,\textsf{D}_5,\textsf{D}_6\)
corresponden a seis (6) conjuntos de datos sintéticos con diferentes
características para probar las bondades de agrupamiento de \(\textsf{M}_1\) y \(\textsf{M}_2\). Estos datos están
inspirados en el artículo Clustering benchmark datasets exploiting
the fundamental clustering problems (https://www.sciencedirect.com/science/article/pii/S2352340920303954),
donde se presenta una colección de conjuntos de datos de fácil acceso en
R para comparar varios métodos de agrupamiento.
Los datos se encuentran disponibles en Moodle por medio de archivos que contienen dos objetos:
Y: matriz de datos de tamaño \(n\times p\) que contiene las observaciones \(p\)-variadas \(\boldsymbol{y}_1,\ldots,\boldsymbol{y}_n\) por filas.xi: concatenación de tamaño \(n\) que contiene la partición de los datos verdadera \(\xi_1,\ldots,\xi_n\) con la que se generaron los datos.
Se observa que la partición dada en xi sirve para
reconocer la número de grupos verdadero \(H\), así como punto de referencia para
evaluar la calidad de la segmentación inferida por medio de los modelos
bajo consideración.
3 Modelamiento
3.1 Visualización de los datos
Hacer una visualización de cada \(\textsf{D}_k\), para \(k=1,\ldots,6\), usando una convención de colores y/o símbolos adecuada que permita identificar fácilmente los grupos verdaderos. Presentar todas las visualizaciones en un solo panel de \(2\times3\) o \(3\times2\) para reconocer las características propias de cada conjunto de datos.
rango <- c(-1,1)*6
col <- RColorBrewer::brewer.pal(n = 9, name = "Set1")
pch <- c(16,15,17,18)
par(mfrow = c(2,3), mar = c(2.75,2.75,1.5,0.5), mgp = c(1.7,0.7,0))
for (j in 1:6) {
load(paste0("D",j,".RData"))
plot(Y, col = adjustcolor(col[xi], 0.5), pch = pch[xi], xlim = rango, ylim = rango, xlab = expression(y[1]), ylab = expression(y[2]), main = paste0("D",j))
abline(h = 0, v = 0, col = "lightgray")
}
3.2 DAGs de los modelos
Representar cada modelo por medio de un DAG. Presentar los dos esquemas en un solo panel de \(2\times1\) para reconocer las características propias de cada modelo.
3.3 Distribuciones condicionales completas para \(\textsf{M}_1\)
Bajo \(\textsf{M}_1\), demostrar las distribuciones condicionales completas de \(\boldsymbol{\theta}_h\), \(\mathbf{\Sigma}\), \(\xi_i\), y \(\boldsymbol{\omega}\).
3.3.1 Distribución condicional completa de \(\boldsymbol{\theta}_h\)
La distribución condicional completa de \(\boldsymbol{\theta}_h\), para \(h = 1,\ldots,H\), es \(\boldsymbol{\theta}_h\mid\text{resto}\sim\small{\textsf{N}}(\boldsymbol{\mu}_n,\mathbf{\Lambda}_n)\) con: \[ \boldsymbol{\mu}_n = \left( \mathbf{\Lambda}_0^{-1} + n_h\mathbf{\Sigma}^{-1} \right)^{-1} \left( \mathbf{\Lambda}_0^{-1}\boldsymbol{\mu}_0 + n_h\mathbf{\Sigma}^{-1}\bar{\boldsymbol{y}}_h \right) \qquad\text{y}\qquad \mathbf{\Lambda}_n = \left( \mathbf{\Lambda}_0^{-1} + n_h\mathbf{\Sigma}^{-1} \right)^{-1} \] donde \(n_h = |\{i:\xi_i=h\}|\) es el número de individuos del grupo \(h\) y \[ \bar{\boldsymbol{y}}_h = \frac{1}{n_h}\sum_{i:\xi_i=h}\boldsymbol{y}_i \] es la media de las observaciones del grupo \(h\). Finalmente, si \(n_h = 0\), entonces \(\boldsymbol{\mu}_n = \boldsymbol{\mu}_0\) y \(\mathbf{\Lambda}_n = \mathbf{\Lambda}_0\).
sample_Theta_M1 <- function (H, nh, Ybh, L0, m0, iL0, iLm0, Theta, Sigma)
{
# Theta: H x p
for (h in 1:H) {
if (nh[h] > 0) {
iS <- solve(Sigma)
Vh <- solve(iL0 + nh[h]*iS)
Theta[h,] <- c(mvtnorm::rmvnorm(1, Vh%*%(iLm0 + nh[h]*(iS%*%Ybh[h,])), Vh))
} else {
Theta[h,] <- c(mvtnorm::rmvnorm(1, m0, L0))
}
}
return(Theta)
}3.3.2 Distribución condicional completa de \(\mathbf{\Sigma}\)
La distribución condicional completa de \(\mathbf{\Sigma}\) es \(\mathbf{\Sigma}\mid\text{resto}\sim\textsf{WI}(\nu_n,\mathbf{S}_n^{-1})\) con: \[ \nu_n = \nu_0 + n \qquad\text{y}\qquad \mathbf{S}_n = \mathbf{S}_0 + \sum_{i=1}^n(\boldsymbol{y}_i - \boldsymbol{\theta}_{\xi_i})(\boldsymbol{y}_i - \boldsymbol{\theta}_{\xi_i})^{\textsf{T}}\,. \]
sample_Sigma_M1 <- function(H, nh, Ybh, Sh, nun, S0, Theta)
{
# Sigma: p x p
tmp <- 0
for(h in 1:H) {
if (nh[h] > 0) {
if (nh[h] == 1) {
tmp <- tmp + S0 + ((Ybh[h,] - Theta[h,])%*%t(Ybh[h,] - Theta[h,]))
} else {
tmp <- tmp + (nh[h]-1)*Sh[,,h] + nh[h]*((Ybh[h,] - Theta[h,])%*%t(Ybh[h,] - Theta[h,]))
}
}
}
return(solve(rWishart::rWishart(1, nun, solve(S0 + tmp))[,,1]))
}3.3.3 Distribución condicional completa de \(\xi_i\)
La distribución condicional completa de \(\xi_i\), para \(i=1,\ldots,n\), es una distribución de probabilidad discreta tal que: \[ p(\xi_i = h\mid\text{resto})\propto \omega_h\, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_h,\mathbf{\Sigma}) \] y por lo tanto \[ \textsf{Pr}(\xi_i = h\mid\text{resto}) = \frac{ \omega_h\, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_h,\mathbf{\Sigma}) }{ \sum_{h=1}^H \omega_h, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_{h},\mathbf{\Sigma}) } \] para \(h = 1,\ldots,H\).
sample_xi_M1 <- function (n, H, xi, omega, Theta, Sigma, Y)
{
# xi: n
for (i in 1:n) {
lp <- NULL
for (h in 1:H) {
lp[h] <- log(omega[h]) + mvtnorm::dmvnorm(Y[i,], Theta[h,], Sigma, log = T)
}
xi[i] <- sample(x = 1:H, size = 1, replace = F, prob = exp(lp - max(lp)))
}
return(xi)
}3.3.4 Distribución condicional completa de \(\boldsymbol{\omega}\)
La distribución condicional completa de \(\boldsymbol{\omega}\) es: \[ \boldsymbol{\omega}\mid\text{resto} \sim \small{\textsf{Dirichlet}}(\alpha^0_1+n_1,\ldots,\alpha^0_{H}+n_H)\,. \]
sample_omega_M1 <- function (nh, al0)
{
# omega: H
return(c(gtools::rdirichlet(1, al0 + nh)))
}3.3.5 Log-verosimilitud
La log-verosimilitud de \(\textsf{M}_1\) es:
\[ \ell = \sum_{i=1}^n\log (\boldsymbol{y}_i\mid\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma})\,. \]
ll_M1 <- function(Y, Theta, Sigma, xi)
{
# escalar
out <- 0
for (i in 1:nrow(Y))
out <- out + mvtnorm::dmvnorm(Y[i,], Theta[xi[i],], Sigma, T)
return(out)
}3.4 Distribuciones condicionales completas para \(\textsf{M}_2\)
Bajo \(\textsf{M}_2\), demostrar las distribuciones condicionales completas de \(\boldsymbol{\theta}_h\), \(\mathbf{\Sigma}_h\), \(\xi_i\), y \(\boldsymbol{\omega}\).
3.4.1 Distribución condicional completa de \(\boldsymbol{\theta}_h\)
La distribución condicional completa de \(\boldsymbol{\theta}_h\), para \(h = 1,\ldots,H\), es \(\boldsymbol{\theta}_h\mid\text{resto}\sim\small{\textsf{N}}(\boldsymbol{\mu}_n,\mathbf{\Lambda}_n)\) con: \[ \boldsymbol{\mu}_n = \left( \mathbf{\Lambda}_0^{-1} + n_h\mathbf{\Sigma}_h^{-1} \right)^{-1} \left( \mathbf{\Lambda}_0^{-1}\boldsymbol{\mu}_0 + n_h\mathbf{\Sigma}_h^{-1}\bar{\boldsymbol{y}}_h \right) \qquad\text{y}\qquad \mathbf{\Lambda}_n = \left( \mathbf{\Lambda}_0^{-1} + n_h\mathbf{\Sigma}_h^{-1} \right)^{-1} \] donde \(n_h = |\{i:\xi_i=h\}|\) es el número de individuos del grupo \(h\) y \[ \bar{\boldsymbol{y}}_h = \frac{1}{n_h}\sum_{i:\xi_i=h}\boldsymbol{y}_i \] es la media de las observaciones del grupo \(h\). Finalmente, si \(n_h = 0\), entonces \(\boldsymbol{\mu}_n = \boldsymbol{\mu}_0\) y \(\mathbf{\Lambda}_n = \mathbf{\Lambda}_0\).
sample_Theta_M2 <- function (H, nh, Ybh, L0, m0, iL0, iLm0, Theta, Sigma)
{
# Theta: H x p
for (h in 1:H) {
if (nh[h] > 0) {
iS <- solve(Sigma[,,h])
Vh <- solve(iL0 + nh[h]*iS)
Theta[h,] <- c(mvtnorm::rmvnorm(1, Vh%*%(iLm0 + nh[h]*(iS%*%Ybh[h,])), Vh))
} else {
Theta[h,] <- c(mvtnorm::rmvnorm(1, m0, L0))
}
}
return(Theta)
}3.4.2 Distribución condicional completa de \(\mathbf{\Sigma}_h\)
La distribución condicional completa de \(\mathbf{\Sigma}_h\) es \(\mathbf{\Sigma}_h\mid\text{resto}\sim\textsf{WI}(\nu_n,\mathbf{S}_n^{-1})\) con: \[ \nu_n = \nu_0 + n_h \qquad\text{y}\qquad \mathbf{S}_n = \mathbf{S}_0 + \sum_{i:\xi_i=h}(\boldsymbol{y}_i - \boldsymbol{\theta}_{h})(\boldsymbol{y}_i - \boldsymbol{\theta}_{h})^{\textsf{T}}\,. \]
sample_Sigma_M2 <- function(H, nh, Ybh, Sh, nu0, S0, iS0, Theta, Sigma)
{
# Sigma: p x p x H
for(h in 1:H) {
if (nh[h] > 0) {
if (nh[h] == 1) {
Sigma[,,h] <- solve(rWishart::rWishart(1, nu0 + 1, solve(S0 + ((Ybh[h,] - Theta[h,])%*%t(Ybh[h,] - Theta[h,]))))[,,1])
} else {
Sigma[,,h] <- solve(rWishart::rWishart(1, nu0 + nh[h], solve(S0 + (nh[h]-1)*Sh[,,h] + nh[h]*((Ybh[h,] - Theta[h,])%*%t(Ybh[h,] - Theta[h,]))))[,,1])
}
} else {
Sigma[,,h] <- solve(rWishart::rWishart(1, nu0, iS0)[,,1])
}
}
return(Sigma)
}3.4.3 Distribución condicional completa de \(\xi_i\)
La distribución condicional completa de \(\xi_i\), para \(i=1,\ldots,n\), es una distribución de probabilidad discreta tal que: \[ p(\xi_i = h\mid\text{resto})\propto \omega_h\, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_h,\mathbf{\Sigma}_h) \] y por lo tanto \[ \textsf{Pr}(\xi_i = h\mid\text{resto}) = \frac{ \omega_h\, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_h,\mathbf{\Sigma}_h) }{ \sum_{h=1}^H \omega_h, \small{\textsf{N}}(\boldsymbol{y}_i\mid\boldsymbol{\theta}_{h},\mathbf{\Sigma}_h) } \] para \(h = 1,\ldots,H\).
sample_xi_M2 <- function (n, H, xi, omega, Theta, Sigma, Y)
{
# xi: n
for (i in 1:n) {
lp <- NULL
for (h in 1:H) {
lp[h] <- log(omega[h]) + mvtnorm::dmvnorm(Y[i,], Theta[h,], Sigma[,,h], log = T)
}
xi[i] <- sample(x = 1:H, size = 1, replace = F, prob = exp(lp - max(lp)))
}
return(xi)
}3.4.4 Distribución condicional completa de \(\boldsymbol{\omega}\)
La distribución condicional completa de \(\boldsymbol{\omega}\) es: \[ \boldsymbol{\omega}\mid\text{resto} \sim \small{\textsf{Dirichlet}}(\alpha^0_1+n_1,\ldots,\alpha^0_{H}+n_H)\,. \]
sample_omega_M2 <- function (nh, al0)
{
# omega: H
c(gtools::rdirichlet(1, al0 + nh))
}3.4.5 Log-verosimilitud
La log-verosimilitud de \(\textsf{M}_2\) es:
\[ \ell = \sum_{i=1}^n\log (\boldsymbol{y}_i\mid\boldsymbol{\theta}_{\xi_i},\mathbf{\Sigma}_{\xi_i})\,. \]
ll_M2 <- function(Y, Theta, Sigma, xi)
{
# escalar
out <- 0
for (i in 1:nrow(Y))
out <- out + mvtnorm::dmvnorm(Y[i,], Theta[xi[i],], Sigma[,,xi[i]], T)
return(out)
}3.5 Ajuste del modelo e inferencia
Ajustar cada uno de los modelos por medio de un muestreador de Gibbs con 1000 iteraciones de calentamiento y 1000 iteraciones efectivas espaciadas cada 10 iteraciones (para un total de 10000 iteraciones) usando el valor de \(H\) dado en la partición de referencia y los hiperparámetros \(\boldsymbol{\alpha}_0=(1/H,\ldots,1/H)\), \(\boldsymbol{\mu}_0=\frac{1}{n}\sum_{i=1}^n \boldsymbol{y}_i\), \(\mathbf{\Lambda}_0=\frac{1}{n-1}\sum_{i=1}^n (\boldsymbol{y}_i-\bar{\boldsymbol{y}})(\boldsymbol{y}_i-\bar{\boldsymbol{y}})^{\textsf{T}}\), \(\nu_0 = p+2\), y \(\mathbf{S}_0 = \mathbf{\Lambda}_0\).
Con el fin de inicializar el algoritmo de manera conveniente para que la cadena sea estacionaria rápidamente, implementar el método de agrupamiento de \(K\)-medias con \(H\) grupos, y fijar \(\boldsymbol{\xi}^{(0)}\) a la partición resultante, \(\boldsymbol{\omega}^{(0)}\) a las frecuencias relativas de los grupos, \(\boldsymbol{\theta}^{(0)}\) a los centros de los grupos, y \(\mathbf{\Sigma}^{(0)}\) al promedio de las matrices de covarianza de los grupos en el caso de \(\textsf{M}_1\) y a las matrices de covarianza de los grupos en el caso de \(\textsf{M}_2\).
Para cada combinación de modelo y conjunto de datos, calcular la media posterior el Índice de Rand Ajustado (ARI, por sus siglas en inglés) junto con el Criterio de Información de Watanabe-Akaike (WAIC, por sus siglas en inglés). Presentar todos los resultados organizadamente en una sola tabla para facilitar la interpretación de los resultados.
3.5.1 Algorito para ajustar \(\textsf{M}_1\)
# muestreador de Gibbs para ajustar M1
MCMC_M1 <- function (niter, nburn, nskip, Y, H, xi0, verbose = TRUE)
{
# numero iteraciones
B <- nburn + nskip*niter
ncat <- floor(0.1*B)
# hiperparametros
n <- nrow(Y)
p <- ncol(Y)
m0 <- colMeans(Y)
L0 <- cov(Y)
nu0 <- p+2
S0 <- L0
al0 <- rep(1/H, H)
# cantidades y valores iniciales
iL0 <- solve(L0)
iLm0 <- solve(L0)%*%m0
nun <- nu0 + n
set.seed(1)
km <- stats::kmeans(x = Y, centers = H)
xi <- c(km$cluster)
omega <- as.numeric(table(factor(xi, levels = 1:H))/n)
Theta <- km$centers
# almacenamiento
THETA <- SIGMA <- OMEGA <- XI <- LP <- ARI <- NULL
# cadena
set.seed(1)
for (b in 1:B) {
# actualizar estadisticos
nh <- as.numeric(table(factor(xi, levels = 1:H)))
Ybh <- array(NA, c(H, p))
Sh <- array(NA, c(p, p, H))
for (h in 1:H) {
if (nh[h] > 0) {
index_h <- xi == h
if (nh[h] == 1) {
Ybh[h,] <- Y[index_h,]
} else {
Ybh[h,] <- colMeans(Y[index_h,])
Sh[,,h] <- cov(Y[index_h,])
}
}
}
# actualizar parametros
Sigma <- sample_Sigma_M1(H, nh, Ybh, Sh, nun, S0, Theta)
Theta <- sample_Theta_M1(H, nh, Ybh, L0, m0, iL0, iLm0, Theta, Sigma)
xi <- sample_xi_M1 (n, H, xi, omega, Theta, Sigma, Y)
omega <- sample_omega_M1(nh, al0)
# almacenar y log-verosimilitud
if (b > nburn) {
if (b%%nskip == 0) {
XI <- rbind(XI, c(xi))
THETA <- rbind(THETA, c(Theta))
SIGMA <- rbind(SIGMA, c(Sigma))
OMEGA <- rbind(OMEGA, c(omega))
LP <- rbind(LP, c(ll_M1(Y, Theta, Sigma, xi)))
ARI <- rbind(ARI, c(aricode::RI (xi, xi0)))
}
}
# progreso
if (verbose) if (b%%ncat == 0) cat(100*round(b/B, 1), "% completado \n", sep = "")
}
# return
return(list(XI = XI, THETA = THETA, SIGMA = SIGMA, OMEGA = OMEGA, LP = LP, ARI = ARI))
}EL WAIC se calcula como: \[ \text{WAIC} = -2\text{lppd} + 2p_{\text{WAIC}} \] donde \[ \begin{align*} \text{lppd} &= \text{log}\prod_{i=1}^n p(\boldsymbol{y}_i\mid\mathbf{Y})\\ &\approx \sum_{i=1}^n\text{log}\left(\frac{1}{B}\sum_{b=1}^B p\left(\boldsymbol{y}_i\mid\boldsymbol{\theta}^{(b)}\right)\right) \end{align*} \] es la log densidad predictiva puntual (que como la devianza resume la precisión predictiva del modelo ajustado a los datos), y \[ \begin{align*} p_{\text{WAIC}} &= 2\sum_{i=1}^n\left( \text{log}\,\textsf{E}(p(\boldsymbol{y}_i\mid\boldsymbol{\theta})\mid\mathbf{Y}) - \textsf{E}(\text{log}\,p(\boldsymbol{y}_i\mid\boldsymbol{\theta})\mid\mathbf{Y}) \right)\\ &\approx 2\sum_{i=1}^n\left(\text{log}\left( \frac{1}{B}\sum_{b=1}^B p\left(\boldsymbol{y}_i\mid\boldsymbol{\theta}^{(b)}\right) \right) - \frac{1}{B}\sum_{b=1}^B \text{log}\,p\left(\boldsymbol{y}_i\mid \boldsymbol{\theta}^{(b)} \right)\right) \end{align*} \] es el número efectivo de parámetros.
WAIC_M1 <- function (Y, muestras)
{
B <- nrow(muestras$XI)
n <- nrow(Y)
p <- ncol(Y)
lppd <- pWAIC <- 0
for (i in 1:n) {
tmp1 <- tmp2 <- 0
for (b in 1:B) {
xi <- muestras$XI[b,]
Theta <- matrix(data = muestras$THETA[b,], nrow = H, ncol = p)
Sigma <- matrix(data = muestras$SIGMA[b,], nrow = p, ncol = p)
aux <- mvtnorm::dmvnorm(x = Y[i,], mean = Theta[xi[i],], sigma = Sigma, log = T)
tmp1 <- tmp1 + exp(aux)/B
tmp2 <- tmp2 + aux/B
}
lppd <- lppd + log(tmp1)
pWAIC <- pWAIC + 2*(log(tmp1) - tmp2)
}
# return
-2*lppd + 2*pWAIC
}3.5.2 Algorito para ajustar \(\textsf{M}_2\)
# muestreador de Gibbs para ajustar M2
MCMC_M2 <- function (niter, nburn, nskip, Y, H, xi0, verbose = TRUE)
{
# numero iteraciones
B <- nburn + nskip*niter
ncat <- floor(0.1*B)
# hiperparametros
n <- nrow(Y)
p <- ncol(Y)
m0 <- colMeans(Y)
L0 <- cov(Y)
nu0 <- p+2
S0 <- L0
al0 <- rep(1/H, H)
iL0 <- solve(L0)
iLm0 <- solve(L0)%*%m0
iS0 <- solve(S0)
# valores iniciales
set.seed(1)
km <- stats::kmeans(x = Y, centers = H)
xi <- c(km$cluster)
omega <- as.numeric(table(factor(xi, levels = 1:H))/n)
Theta <- km$centers
Sigma <- array(data = NA, dim = c(p,p,H))
for (h in 1:H)
Sigma[,,h] <- cov(Y[xi == h,])
# almacenamiento
THETA <- SIGMA <- OMEGA <- XI <- LP <- ARI <- NULL
# cadena
set.seed(1)
for (b in 1:B) {
# actualizar estadisticos
nh <- as.numeric(table(factor(xi, levels = 1:H)))
Ybh <- array(NA, c(H, p))
Sh <- array(NA, c(p, p, H))
for (h in 1:H) {
if (nh[h] > 0) {
index_h <- xi == h
if (nh[h] == 1) {
Ybh[h,] <- Y[index_h,]
} else {
Ybh[h,] <- colMeans(Y[index_h,])
Sh[,,h] <- cov(Y[index_h,])
}
}
}
# actualizar parametros
Sigma <- sample_Sigma_M2(H, nh, Ybh, Sh, nu0, S0, iS0, Theta, Sigma)
Theta <- sample_Theta_M2(H, nh, Ybh, L0, m0, iL0, iLm0, Theta, Sigma)
xi <- sample_xi_M2 (n, H, xi, omega, Theta, Sigma, Y)
omega <- sample_omega_M2(nh, al0)
# almacenar y log-verosimilitud
if (b > nburn) {
if (b%%nskip == 0) {
XI <- rbind(XI, c(xi))
THETA <- rbind(THETA, c(Theta))
SIGMA <- rbind(SIGMA, c(Sigma))
OMEGA <- rbind(OMEGA, c(omega))
LP <- rbind(LP, ll_M2(Y, Theta, Sigma, xi))
ARI <- rbind(ARI, c(aricode::RI (xi, xi0)))
}
}
# progreso
if (verbose) if (b%%ncat == 0) cat(100*round(b/B, 1), "% completado \n", sep = "")
}
# return
return(list(XI = XI, THETA = THETA, SIGMA = SIGMA, OMEGA = OMEGA, LP = LP, ARI = ARI))
}WAIC_M2 <- function (Y, muestras)
{
B <- nrow(muestras$XI)
n <- nrow(Y)
p <- ncol(Y)
lppd <- pWAIC <- 0
for (i in 1:n) {
tmp1 <- tmp2 <- 0
for (b in 1:B) {
xi <- muestras$XI[b,]
Theta <- array(data = muestras$THETA[b,], dim = c(H,p))
Sigma <- array(data = muestras$SIGMA[b,], dim = c(p,p,H))
aux <- mvtnorm::dmvnorm(x = Y[i,], mean = Theta[xi[i],], sigma = Sigma[,,xi[i]], log = T)
tmp1 <- tmp1 + exp(aux)/B
tmp2 <- tmp2 + aux/B
}
lppd <- lppd + log(tmp1)
pWAIC <- pWAIC + 2*(log(tmp1) - tmp2)
}
# return
-2*lppd + 2*pWAIC
}3.5.3 Ajuste de \(\textsf{M}_1\) y \(\textsf{M}_2\)
# iteraciones
niter <- 1000
nburn <- 1000
nskip <- 10# ajuste
for (j in 1:6) {
# datos
load(file = paste0("D",j,".RData"))
H <- max(xi)
for (k in 1:2) {
nombre <- paste0("muestras","_",j,"_","M",k)
mcmc <- get(x = paste0("MCMC_M",k))
waic <- get(x = paste0("WAIC_M",k))
assign(x = nombre, value = mcmc(niter, nburn, nskip, Y, H, xi, verbose = FALSE))
muestras <- get(x = nombre)
criterio <- waic(Y, muestras)
save(muestras, criterio, file = paste0(nombre,".RData"))
}
}3.5.4 Consolidación de resultados ARI y WAIC para \(\textsf{M}_1\) y \(\textsf{M}_2\)
tab <- NULL
for (j in 1:6) {
tmp <- NULL
for (k in 1:2) {
load(file = paste0("muestras","_",j,"_","M",k,".RData"))
tmp <- c(tmp, mean(muestras$ARI), criterio)
}
tab <- rbind(tab, tmp)
}
colnames(tab) <- c("M1 : ARI", "M1 : WAIC", "M2 : ARI", "M2 : WAIC")
rownames(tab) <- c(paste0("D",1:6))
knitr::kable(x = tab, digits = c(3,1), align = "c")| M1 : ARI | M1 : WAIC | M2 : ARI | M2 : WAIC | |
|---|---|---|---|---|
| D1 | 0.914 | 929.1 | 0.904 | 936.9 |
| D2 | 0.922 | 918.1 | 0.942 | 847.4 |
| D3 | 0.898 | 635.3 | 0.940 | 588.0 |
| D4 | 0.912 | 526.7 | 0.966 | 481.3 |
| D5 | 0.828 | 1303.6 | 0.844 | 1271.4 |
| D6 | 0.981 | 478.2 | 0.978 | 484.1 |
3.6 Matrices de incidencia bajo \(\textsf{M}_1\)
Bajo \(\textsf{M}_1\), hacer un diagrama de calor de cada matriz de incidencia a posteriori ordenada de acuerdo con la partición de referencia. Presentar todas las visualizaciones en un solo panel de \(2\times3\) o \(3\times2\) para reconocer las características propias de cada matriz de incidencia.
get_incidence_matrix <- function(XI)
{
# matriz de incidencia
n <- ncol(XI)
B <- nrow(XI)
A <- matrix(data = 0, nrow = n, ncol = n)
for (b in 1:B) {
for (i in 1:(n-1)) {
for (j in (i+1):n) {
if (XI[b,i] == XI[b,j]) {
A[i,j] <- A[i,j] + 1/B
}
}
}
}
A <- A + t(A)
diag(A) <- 1
return(A)
}heat.plot0 <- function (mat, show.grid = FALSE, cex.axis, tick, labs, col.axis,...)
{
# funcion para graficar las probabilidades por medio de un diagrama de calor
JJ <- dim(mat)[1]
colorscale <- c("white", rev(heat.colors(100)))
if(missing(labs)) labs <- 1:JJ
if(missing(col.axis)) col.axis <- rep("black", JJ)
if(missing(cex.axis)) cex.axis <- 0.5
if(missing(tick)) tick <- TRUE
## adjacency matrix
image(seq(1, JJ), seq(1, JJ), mat, axes = FALSE, xlab = "", ylab = "",
col = colorscale[seq(floor(100*min(mat)), floor(100*max(mat)))], ...)
for(j in 1:JJ){
axis(1, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j])
axis(2, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j])
}
box()
if(show.grid) grid(nx = JJ, ny = JJ)
}par(mfrow = c(2,3), mar = c(2.75,2.75,1.5,0.5), mgp = c(1.7,0.7,0))
for (j in 1:6) {
load(file = paste0("muestras","_",j,"_","M1.RData"))
A <- get_incidence_matrix(muestras$XI)
heat.plot0(mat = A, labs = NA, tick = F, main = paste("M1 : D",j))
}
3.7 Matrices de incidencia bajo \(\textsf{M}_2\)
Bajo \(\textsf{M}_2\), hacer un diagrama de calor de cada matriz de incidencia a posteriori ordenada de acuerdo con la partición de referencia. Presentar todas las visualizaciones en un solo panel de \(2\times3\) o \(3\times2\) para reconocer las características propias de cada matriz de incidencia.
par(mfrow = c(2,3), mar = c(2.75,2.75,1.5,0.5), mgp = c(1.7,0.7,0))
for (j in 1:6) {
load(file = paste0("muestras","_",j,"_","M2.RData"))
A <- get_incidence_matrix(muestras$XI)
heat.plot0(mat = A, labs = NA, tick = F, main = paste("M2 : D",j))
}
4 Conclusiones
Interpretar los hallazgos obtenidos en los numerales 5., 6., y 7.. Utilizar máximo 1000 palabras.