Generalmente se considera que la regresión lineal no es lo suficientemente flexible para abordar datos no lineales. Desde el punto de vista teórico no es capaz de abordarlos. Sin embargo, podemos hacer que funcione para nosotros con cualquier conjunto de datos utilizando mezclas normales finitas en un modelo de regresión. De esta manera, se convierte en una herramienta de aprendizaje automático muy poderosa que se puede aplicar a prácticamente cualquier conjunto de datos, incluso a los muy anormales con dependencias no lineales entre las variables.
Lo que hace que este enfoque sea particularmente interesante es la interpretabilidad. A pesar de un nivel de flexibilidad extremadamente alto, todas las relaciones detectadas pueden interpretarse directamente. El modelo es tan general como una red neuronal, pero aún así no se convierte en una caja negra. Puede leer las relaciones y comprender el impacto de las variables individuales.
En esta publicación, demostramos cómo simular un modelo de mezcla finita para regresión utilizando el muestreo Markov Chain Monte Carlo (MCMC). Generaremos datos con múltiples componentes (grupos) y ajustaremos un modelo mixto para recuperar estos componentes utilizando la inferencia bayesiana. Este proceso involucra modelos de regresión y modelos mixtos, combinándolos con técnicas MCMC para la estimación de parámetros.
Comenzamos cargando las bibliotecas necesarias para trabajar con modelos de regresión, MCMC y distribuciones multivariadas.
# Loading the required libraries for various functions
library("pscl") # For pscl specific functions, like regression models
library("MCMCpack") # For MCMC sampling functions, including posterior distributions
library(mvtnorm) # For multivariate normal distribution functio
- pscl: Se utiliza para diversas funciones estadísticas como modelos de regresión.
- Paquete MCMC: Contiene funciones para la inferencia bayesiana, en particular el muestreo MCMC.
- mvtnorm: Proporciona herramientas para trabajar con distribuciones normales multivariadas.
Simulamos un conjunto de datos donde cada observación pertenece a uno de varios grupos (componentes del modelo mixto) y la variable de respuesta se genera utilizando un modelo de regresión con coeficientes aleatorios.
Consideramos una configuración general para un modelo de regresión que utiliza componentes de mezcla G Normal.
## Generate the observations
# Set the length of the time series (number of observations per group)
N <- 1000
# Set the number of simulations (iterations of the MCMC process)
nSim <- 200
# Set the number of components in the mixture model (G is the number of groups)
G <- 3
- norte: El número de observaciones por grupo.
- nSim: El número de iteraciones de MCMC.
- GRAMO: El número de componentes (grupos) en nuestro modelo de mezcla.
Simulación de datos
Cada grupo se modela utilizando un modelo de regresión univariado, donde las variables explicativas (X) y la variable de respuesta (y) se simulan a partir de distribuciones normales. El betas
representar los coeficientes de regresión para cada grupo, y sigmas
representan la varianza de cada grupo.
# Set the values for the regression coefficients (betas) for each group
betas <- 1:sum(dimG) * 2.5 # Generating sequential betas with a multiplier of 2.5
# Define the variance (sigma) for each component (group) in the mixture
sigmas <- rep(1, G) / 1 # Set variance to 1 for each component, with a fixed divisor of 1
- betas: Estos son los coeficientes de regresión. El coeficiente de cada grupo se asigna secuencialmente.
- sigmas: Representa la varianza de cada grupo en el modelo de mezcla.
En este modelo permitimos que cada componente de la mezcla posea su propio parámetro de varianza y un conjunto de parámetros de regresión.
Asignación de grupos y mezcla
Luego simulamos la asignación grupal de cada observación mediante una asignación aleatoria y mezclamos los datos de todos los componentes.
Aumentamos el modelo con un conjunto de vectores de etiquetas de componentes para
dónde
y así z_gi=1 implica que el i-El individuo se extrae de la gramo-º componente de la mezcla.
Esta asignación aleatoria forma la z_original
vector, que representa el verdadero grupo al que pertenece cada observación.
# Initialize the original group assignments (z_original)
z_original <- matrix(NA, N * G, 1)
# Repeat each group label N times (assign labels to each observation per group)
z_original <- rep(1:G, rep(N, G))
# Resample the data rows by random order
sampled_order <- sample(nrow(data))
# Apply the resampled order to the data
data <- data[sampled_order,]
Establecimos distribuciones previas para los coeficientes de regresión y las varianzas. Estos antecedentes guiarán nuestra estimación bayesiana.
## Define Priors for Bayesian estimation# Define the prior mean (muBeta) for the regression coefficients
muBeta <- matrix(0, G, 1)# Define the prior variance (VBeta) for the regression coefficients
VBeta <- 100 * diag(G) # Large variance (100) as a prior for the beta coefficients# Prior for the sigma parameters (variance of each component)
ag <- 3 # Shape parameter
bg <- 1/2 # Rate parameter for the prior on sigma
shSigma <- ag
raSigma <- bg^(-1)
- mubeta: La media anterior de los coeficientes de regresión. Lo configuramos en 0 para todos los componentes.
- V Beta: La varianza anterior, que es grande (100) para permitir flexibilidad en los coeficientes.
- shSigma y raSigma: Parámetros de forma y tasa para el anterior en la varianza (sigma) de cada grupo.
Para los indicadores y probabilidades de los componentes consideramos la siguiente asignación previa.
El multinomial previo M son las generalizaciones multivariadas del binomio, y el previo de Dirichlet D es una generalización multivariada de la distribución beta.
En esta sección, inicializamos el proceso MCMC configurando matrices para almacenar las muestras de los coeficientes de regresión, las varianzas y las proporciones de mezcla.
## Initialize MCMC sampling# Initialize matrix to store the samples for beta
mBeta <- matrix(NA, nSim, G)# Assign the first value of beta using a random normal distribution
for (g in 1:G) {
mBeta[1, g] <- rnorm(1, muBeta[g, 1], VBeta[g, g])
}# Initialize the sigma^2 values (variance for each component)
mSigma2 <- matrix(NA, nSim, G)
mSigma2[1, ] <- rigamma(1, shSigma, raSigma)# Initialize the mixing proportions (pi), using a Dirichlet distribution
mPi <- matrix(NA, nSim, G)
alphaPrior <- rep(N/G, G) # Prior for the mixing proportions, uniform across groups
mPi[1, ] <- rdirichlet(1, alphaPrior)
- mbeta: Matriz para almacenar muestras de los coeficientes de regresión.
- mSigma2: Matriz para almacenar las varianzas (sigma cuadrado) de cada componente.
- MPi: Matriz para almacenar las proporciones de mezcla, inicializada mediante una distribución de Dirichlet.
Si condicionamos los valores de las variables indicadoras componentes z, la probabilidad condicional se puede expresar como
En el bucle de muestreo de MCMC, actualizamos las asignaciones de grupo (z
), coeficientes de regresión (beta
), y variaciones (sigma
) basado en las distribuciones posteriores. Se calcula la probabilidad de cada asignación de grupo y se selecciona el grupo con la probabilidad posterior más alta.
Se pueden obtener los siguientes condicionales posteriores completos:
dónde
denota todos los parámetros en nuestro posterior excepto incógnita.
y donde negro denota el número de observaciones en el gramo-ésimo componente de la mezcla.
y
El siguiente algoritmo se basa en la serie de distribuciones posteriores anteriores en orden secuencial.
## Start the MCMC iterations for posterior sampling# Loop over the number of simulations
for (i in 2:nSim) {
print(i) # Print the current iteration number# For each observation, update the group assignment (z)
for (t in 1:(N*G)) {
fig <- NULL
for (g in 1:G) {
# Calculate the likelihood of each group and the corresponding posterior probability
fig[g] <- dnorm(y[t, 1], X[t, ] %*% mBeta[i-1, g], sqrt(mSigma2[i-1, g])) * mPi[i-1, g]
}
# Avoid zero likelihood and adjust it
if (all(fig) == 0) {
fig <- fig + 1/G
}
# Sample a new group assignment based on the posterior probabilities
z[i, t] <- which(rmultinom(1, 1, fig/sum(fig)) == 1)
}
# Update the regression coefficients for each group
for (g in 1:G) {
# Compute the posterior mean and variance for beta (using the data for group g)
DBeta <- solve(t(X[z[i, ] == g, ]) %*% X[z[i, ] == g, ] / mSigma2[i-1, g] + solve(VBeta[g, g]))
dBeta <- t(X[z[i, ] == g, ]) %*% y[z[i, ] == g, 1] / mSigma2[i-1, g] + solve(VBeta[g, g]) %*% muBeta[g, 1]
# Sample a new value for beta from the multivariate normal distribution
mBeta[i, g] <- rmvnorm(1, DBeta %*% dBeta, DBeta)
# Update the number of observations in group g
ng[i, g] <- sum(z[i, ] == g)
# Update the variance (sigma^2) for each group
mSigma2[i, g] <- rigamma(1, ng[i, g]/2 + shSigma, raSigma + 1/2 * sum((y[z[i, ] == g, 1] - (X[z[i, ] == g, ] * mBeta[i, g]))^2))
}
# Reorder the group labels to maintain consistency
reorderWay <- order(mBeta[i, ])
mBeta[i, ] <- mBeta[i, reorderWay]
ng[i, ] <- ng[i, reorderWay]
mSigma2[i, ] <- mSigma2[i, reorderWay]
# Update the mixing proportions (pi) based on the number of observations in each group
mPi[i, ] <- rdirichlet(1, alphaPrior + ng[i, ])
}
Este bloque de código realiza los pasos clave en MCMC:
- Actualización de tareas de grupo: Para cada observación, calculamos la probabilidad de que los datos pertenezcan a cada grupo y actualizamos la asignación del grupo en consecuencia.
- Actualización del coeficiente de regresión: Los coeficientes de regresión para cada grupo se actualizan utilizando la media y la varianza posteriores, que se calculan en función de los datos observados.
- Actualización de variación: La varianza de la variable de respuesta para cada grupo se actualiza utilizando la distribución gamma inversa.
Finalmente, visualizamos los resultados del muestreo de MCMC. Trazamos las distribuciones posteriores para cada coeficiente de regresión, las comparamos con los valores verdaderos y trazamos las asignaciones de grupo más probables.
# Plot the posterior distributions for each beta coefficient
par(mfrow=c(G,1))
for (g in 1:G) {
plot(density(mBeta[5:nSim, g]), main = 'True parameter (vertical) and the distribution of the samples') # Plot the density for the beta estimates
abline(v = betas[g]) # Add a vertical line at the true value of beta for comparison
}
Este gráfico muestra cómo las muestras de MCMC (distribución posterior) para los coeficientes de regresión convergen a los valores verdaderos (betas
).
A través de este proceso, demostramos cómo se pueden usar mezclas normales finitas en un contexto de regresión, combinadas con MCMC para la estimación de parámetros. Al simular datos con agrupaciones conocidas y recuperar los parámetros mediante inferencia bayesiana, podemos evaluar qué tan bien nuestro modelo captura la estructura subyacente de los datos.
A menos que se indique lo contrario, todas las imágenes son del autor.