Introducción

El análisis de agrupamiento (clustering) es una técnica de aprendizaje no supervisado cuyo objetivo es agrupar observaciones similares en un mismo conjunto. Uno de los métodos más utilizados es el algoritmo K-means, el cual se basa en minimizar la variabilidad dentro de cada grupo utilizando una métrica de distancia.

En este documento se presentan:

Métrica inducida por producto interno

Sea \(V\) un espacio vectorial real. Un producto interno es una función

\[ \langle \cdot , \cdot \rangle : V \times V \to \mathbb{R} \]

que satisface:

  1. Simetría

\[ \langle x,y \rangle = \langle y,x \rangle \]

  1. Linealidad

\[ \langle ax + by, z \rangle = a\langle x,z \rangle + b\langle y,z \rangle \]

  1. Definida positiva

\[ \langle x,x \rangle \ge 0 \]

y

\[ \langle x,x \rangle =0 \iff x=0 \]

A partir del producto interno se define la norma

\[ ||x|| = \sqrt{\langle x,x \rangle} \]

y la distancia inducida

\[ d(x,y)=||x-y||=\sqrt{\langle x-y,x-y\rangle} \]

En \(\mathbb{R}^p\) el producto interno usual es

\[ \langle x,y\rangle = x^Ty \]

lo que induce la distancia euclidiana.

Distancia de Mahalanobis

La distancia de Mahalanobis generaliza la distancia euclidiana incorporando la estructura de covarianza de los datos.

Sea \(X \in \mathbb{R}^p\) con matriz de covarianza \(\Sigma\) definida positiva.

La distancia de Mahalanobis entre dos vectores \(x\) y \(y\) se define como

\[ d_M(x,y)=\sqrt{(x-y)^T \Sigma^{-1} (x-y)} \]

Propiedades importantes:

  1. Tiene en cuenta la correlación entre variables.
  2. Es invariante ante transformaciones lineales.
  3. Si \(\Sigma = I\), se reduce a la distancia euclidiana.

También puede interpretarse como una métrica inducida por el producto interno

\[ \langle x,y \rangle_{\Sigma^{-1}} = x^T\Sigma^{-1}y \]

Algoritmo K-means

Sea un conjunto de datos

\[ X=\{x_1,\dots,x_n\}, \quad x_i \in \mathbb{R}^p \]

El objetivo es particionar los datos en \(K\) grupos

\[ C_1,\dots,C_K \]

minimizando la suma de cuadrados dentro del grupo:

\[ \min_{C_1,\dots,C_K} \sum_{k=1}^{K}\sum_{x_i\in C_k} ||x_i-\mu_k||^2 \]

donde

\[ \mu_k=\frac{1}{|C_k|}\sum_{x_i\in C_k} x_i \]

es el centroide del cluster.

Algoritmo iterativo

El algoritmo procede de la siguiente forma:

  1. Inicializar \(K\) centroides.
  2. Asignación

\[ x_i \rightarrow C_k \quad \text{si} \quad ||x_i-\mu_k|| \le ||x_i-\mu_j|| \]

para todo \(j\).

  1. Actualización

\[ \mu_k = \frac{1}{|C_k|}\sum_{x_i \in C_k} x_i \]

  1. Repetir hasta convergencia.

El algoritmo converge a un mínimo local de la función objetivo.

K-means con distancia de Mahalanobis

Si se usa una métrica inducida por una matriz definida positiva \(A\), la función objetivo se vuelve

\[ \sum_{k=1}^{K}\sum_{x_i\in C_k} (x_i-\mu_k)^T A (x_i-\mu_k) \]

Cuando

\[ A=\Sigma^{-1} \]

se obtiene una versión basada en distancia de Mahalanobis, lo que permite clusters elípticos en lugar de esféricos.

Ejemplo en R

set.seed(123)

x <- matrix(rnorm(200), ncol=2)

k <- kmeans(x, centers=3)

plot(x, col=k$cluster, pch=19)

x <- matrix(rnorm(100), ncol=2)

S <- cov(x)
mu <- colMeans(x)

mahalanobis(x, center=mu, cov=S)
##  [1] 0.404402444 0.568227294 2.192849900 0.298207527 1.442613139 0.505167388
##  [7] 4.182298315 3.461421983 0.707409389 0.752495347 2.374952797 1.497054146
## [13] 5.808750994 5.860433970 0.299057906 4.815718293 2.088430156 1.861019226
## [19] 0.720030414 1.214795720 2.095295537 3.870466031 0.513579529 0.556090272
## [25] 0.485188679 0.526724269 1.292158666 1.142990008 9.767816117 0.272358318
## [31] 0.148579199 0.798392949 0.310524977 1.953492766 0.009820559 2.467099292
## [37] 0.476031020 0.615802932 1.617174001 0.966884846 4.142029025 2.097192323
## [43] 3.440335968 5.296549827 2.916784181 4.743751468 0.466279235 1.533329236
## [49] 1.086363682 1.335578737
points(k$centers, col=1:3, pch=8, cex=2)

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
data(iris)

X <- iris[,1:4]

set.seed(123)
k <- kmeans(X, centers = 3)

iris$cluster <- as.factor(k$cluster)

mu <- colMeans(X)
S <- cov(X)

iris$mahalanobis <- mahalanobis(X, center = mu, cov = S)

p <- ncol(X)
iris$chi_value <- pchisq(iris$mahalanobis, df = p)

centers <- as.data.frame(k$centers)
centers$cluster <- as.factor(1:3)

g <- ggplot(iris, aes(Petal.Length, Petal.Width, color = cluster)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_point(data = centers,
             aes(Petal.Length, Petal.Width),
             color = "black",
             size = 6,
             shape = 8) +
  labs(
    title = "K-means con Iris y distancia de Mahalanobis",
    x = "Petal Length",
    y = "Petal Width",
    color = "Cluster"
  ) +
  theme_minimal()

print(g)

head(iris[,c("cluster","mahalanobis","chi_value")])
##   cluster mahalanobis chi_value
## 1       1    2.134468 0.2889571
## 2       1    2.849119 0.4166159
## 3       1    2.081339 0.2791985
## 4       1    2.452382 0.3468176
## 5       1    2.462155 0.3485748
## 6       1    3.883418 0.5779866
library(ggplot2)

data(iris)

X <- iris[,1:4]

set.seed(123)
k <- kmeans(X, centers = 3)

iris$cluster <- as.factor(k$cluster)

mu <- colMeans(X)
S <- cov(X)

iris$mahal <- mahalanobis(X, center = mu, cov = S)

p <- ncol(X)
iris$chi_prob <- pchisq(iris$mahal, df = p)

centers <- as.data.frame(k$centers)
centers$cluster <- as.factor(1:3)

g <- ggplot(iris, aes(Petal.Length, Petal.Width, color = cluster)) +
  geom_point(size = 3, alpha = 0.8) +
  stat_ellipse(aes(group = cluster), level = 0.95, linewidth = 1.2) +
  geom_point(data = centers,
             aes(Petal.Length, Petal.Width),
             color = "black",
             shape = 8,
             size = 5) +
  labs(
    title = "K-means en Iris con elipses tipo Mahalanobis",
    x = "Petal Length",
    y = "Petal Width",
    color = "Cluster"
  ) +
  theme_minimal()

print(g)

head(iris[,c("cluster","mahal","chi_prob")])
##   cluster    mahal  chi_prob
## 1       1 2.134468 0.2889571
## 2       1 2.849119 0.4166159
## 3       1 2.081339 0.2791985
## 4       1 2.452382 0.3468176
## 5       1 2.462155 0.3485748
## 6       1 3.883418 0.5779866

Introducción

En estadística multivariada la distribución normal multivariada juega un papel central. Muchas técnicas estadísticas, como análisis discriminante, análisis de componentes principales y métodos de clustering, dependen de sus propiedades.

Una de las cantidades fundamentales asociadas a esta distribución es la distancia de Mahalanobis, la cual mide la distancia entre un vector aleatorio y su media teniendo en cuenta la estructura de covarianza.

En este documento se presenta:

Distribución normal multivariada

Definición

Un vector aleatorio

\[ X = (X_1, X_2, \dots, X_p)^T \]

tiene distribución normal multivariada con media \(\mu\) y matriz de covarianza \(\Sigma\) si su función de densidad es

\[ f(x) = \frac{1}{(2\pi)^{p/2} |\Sigma|^{1/2}} \exp\left( -\frac{1}{2}(x-\mu)^T \Sigma^{-1} (x-\mu) \right) \]

donde

  • \(\mu \in \mathbb{R}^p\) es el vector de medias
  • \(\Sigma\) es una matriz simétrica definida positiva

Se denota

\[ X \sim N_p(\mu,\Sigma) \]

Propiedades de la normal multivariada

Esperanza

\[ E(X) = \mu \]

Covarianza

\[ Cov(X) = \Sigma \]

Transformación lineal

Si

\[ Y = AX + b \]

entonces

\[ Y \sim N_q(A\mu + b, A\Sigma A^T) \]

donde \(A\) es una matriz \(q\times p\).

Estandarización

Sea

\[ X \sim N_p(\mu,\Sigma) \]

y sea

\[ Z = \Sigma^{-1/2}(X-\mu) \]

Entonces

\[ Z \sim N_p(0,I) \]

donde \(I\) es la matriz identidad.

Esta transformación se llama blanqueamiento o whitening transformation.

Distancia de Mahalanobis

La distancia de Mahalanobis entre un punto \(x\) y la media \(\mu\) se define como

\[ D^2 = (x-\mu)^T \Sigma^{-1} (x-\mu) \]

Esta cantidad mide la distancia teniendo en cuenta:

Si \(\Sigma = I\), entonces

\[ D^2 = ||x-\mu||^2 \]

que corresponde a la distancia euclidiana al cuadrado.

Teorema fundamental

Teorema

Si

\[ X \sim N_p(\mu,\Sigma) \]

entonces

\[ D^2 = (X-\mu)^T \Sigma^{-1} (X-\mu) \]

sigue una distribución

\[ D^2 \sim \chi^2_p \]

Demostración

Partimos de la estandarización

\[ Z = \Sigma^{-1/2}(X-\mu) \]

Sabemos que

\[ Z \sim N_p(0,I) \]

Ahora observamos que

\[ D^2 = (X-\mu)^T\Sigma^{-1}(X-\mu) \]

pero

\[ \Sigma^{-1} = (\Sigma^{-1/2})^T \Sigma^{-1/2} \]

entonces

\[ D^2 = (X-\mu)^T (\Sigma^{-1/2})^T \Sigma^{-1/2}(X-\mu) \]

lo cual puede escribirse como

\[ D^2 = Z^T Z \]

Ahora

\[ Z^T Z = \sum_{i=1}^p Z_i^2 \]

donde

\[ Z_i \sim N(0,1) \]

independientes.

Por definición de la distribución chi-cuadrado

\[ \sum_{i=1}^p Z_i^2 \sim \chi^2_p \]

por lo tanto

\[ D^2 \sim \chi^2_p \]

Regiones elípticas

Las superficies de nivel de la densidad normal multivariada están dadas por

\[ (x-\mu)^T \Sigma^{-1}(x-\mu) = c \]

Estas ecuaciones representan elipsoides en \(\mathbb{R}^p\).

Si

\[ c = \chi^2_{p,\alpha} \]

entonces

\[ P(D^2 \le c) = \alpha \]

lo que define regiones de confianza elípticas.

Ejemplo en R

data(iris)

X <- iris[,1:4]

mu <- colMeans(X)
S <- cov(X)

d <- mahalanobis(X, center = mu, cov = S)

head(d)
## [1] 2.134468 2.849119 2.081339 2.452382 2.462155 3.883418

Introducción

En estadística multivariada muchas técnicas suponen que el vector aleatorio sigue una distribución normal multivariada. Entre estos métodos se encuentran:

Por esta razón es importante verificar el supuesto de normalidad multivariada.

Sea

\[ X = (X_1,\dots,X_p)^T \]

un vector aleatorio. Decimos que

\[ X \sim N_p(\mu,\Sigma) \]

si su densidad es

\[ f(x)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}} \exp\left(-\frac12 (x-\mu)^T\Sigma^{-1}(x-\mu)\right) \]

donde

Distancia de Mahalanobis

Una herramienta fundamental para estudiar la normalidad multivariada es la distancia de Mahalanobis.

Para una observación \(x_i\) se define

\[ D_i^2 = (x_i-\mu)^T\Sigma^{-1}(x_i-\mu) \]

Si

\[ X \sim N_p(\mu,\Sigma) \]

entonces

\[ D_i^2 \sim \chi^2_p \]

Esto permite construir diagnósticos gráficos y pruebas estadísticas.

Gráfico Chi-cuadrado

Una forma visual de evaluar la normalidad multivariada consiste en ordenar las distancias de Mahalanobis

\[ D_{(1)}^2 \le D_{(2)}^2 \le \dots \le D_{(n)}^2 \]

y compararlas con los cuantiles teóricos

\[ \chi^2_{p,\frac{i-0.5}{n}} \]

Si los datos provienen de una normal multivariada, los puntos deben alinearse aproximadamente en una recta.

Prueba de Mardia

La prueba de Mardia se basa en asimetría multivariada y curtosis multivariada.

Sea

\[ S \]

la matriz de covarianza muestral.

Asimetría multivariada

La medida de asimetría de Mardia es

\[ b_{1,p} = \frac{1}{n^2}\sum_{i=1}^n\sum_{j=1}^n [(x_i-\bar x)^TS^{-1}(x_j-\bar x)]^3 \]

Bajo normalidad multivariada

\[ n b_{1,p}/6 \sim \chi^2_{p(p+1)(p+2)/6} \]

Curtosis multivariada

La curtosis de Mardia es

\[ b_{2,p} = \frac{1}{n}\sum_{i=1}^n [(x_i-\bar x)^TS^{-1}(x_i-\bar x)]^2 \]

Bajo normalidad

\[ E(b_{2,p}) = p(p+2) \]

y el estadístico

\[ z = \frac{b_{2,p}-p(p+2)}{\sqrt{\frac{8p(p+2)}{n}}} \]

es aproximadamente normal estándar.

Prueba de Henze–Zirkler

La prueba de Henze–Zirkler es una de las más utilizadas porque tiene buena potencia frente a muchas alternativas.

Se basa en comparar la función característica empírica con la de una normal multivariada.

El estadístico general tiene la forma

\[ HZ = n \left[ \frac{1}{n^2}\sum_{i,j} e^{-\beta^2 d_{ij}/2} -2(1+\beta^2)^{-p/2} \frac{1}{n}\sum_i e^{-\frac{\beta^2}{2(1+\beta^2)}d_i} +(1+2\beta^2)^{-p/2} \right] \]

donde

Prueba de Royston

La prueba de Royston es una extensión multivariada del test de Shapiro-Wilk.

Consiste en aplicar transformaciones a las variables y combinar los estadísticos univariados para obtener un test global de normalidad.

Es particularmente útil cuando el tamaño de muestra es pequeño.

Ejemplo en R

data(iris)

X <- iris[,1:4]

mu <- colMeans(X)
S <- cov(X)

d <- mahalanobis(X, center = mu, cov = S)

head(d)
## [1] 2.134468 2.849119 2.081339 2.452382 2.462155 3.883418
p <- ncol(X)
n <- nrow(X)

d_sorted <- sort(d)
chi <- qchisq((1:n-0.5)/n, df=p)

plot(chi, d_sorted,
     xlab="Cuantiles Chi-cuadrado",
     ylab="Distancia Mahalanobis",
     main="QQ plot para normalidad multivariada")

abline(0,1,col="red")

library(MVN)

mvn(data = X, mvnTest = "mardia")
## $multivariateNormality
##              Test          Statistic              p value Result
## 1 Mardia Skewness    67.430508778062 4.75799820400869e-07     NO
## 2 Mardia Kurtosis -0.230112114481001    0.818004651478012    YES
## 3             MVN               <NA>                 <NA>     NO
## 
## $univariateNormality
##               Test     Variable Statistic   p value Normality
## 1 Anderson-Darling Sepal.Length    0.8892  0.0225      NO    
## 2 Anderson-Darling Sepal.Width     0.9080  0.0202      NO    
## 3 Anderson-Darling Petal.Length    7.6785  <0.001      NO    
## 4 Anderson-Darling Petal.Width     5.1057  <0.001      NO    
## 
## $Descriptives
##                n     Mean   Std.Dev Median Min Max 25th 75th       Skew
## Sepal.Length 150 5.843333 0.8280661   5.80 4.3 7.9  5.1  6.4  0.3086407
## Sepal.Width  150 3.057333 0.4358663   3.00 2.0 4.4  2.8  3.3  0.3126147
## Petal.Length 150 3.758000 1.7652982   4.35 1.0 6.9  1.6  5.1 -0.2694109
## Petal.Width  150 1.199333 0.7622377   1.30 0.1 2.5  0.3  1.8 -0.1009166
##                Kurtosis
## Sepal.Length -0.6058125
## Sepal.Width   0.1387047
## Petal.Length -1.4168574
## Petal.Width  -1.3581792
mvn(data = X, mvnTest = "hz")
## $multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 2.336394       0  NO
## 
## $univariateNormality
##               Test     Variable Statistic   p value Normality
## 1 Anderson-Darling Sepal.Length    0.8892  0.0225      NO    
## 2 Anderson-Darling Sepal.Width     0.9080  0.0202      NO    
## 3 Anderson-Darling Petal.Length    7.6785  <0.001      NO    
## 4 Anderson-Darling Petal.Width     5.1057  <0.001      NO    
## 
## $Descriptives
##                n     Mean   Std.Dev Median Min Max 25th 75th       Skew
## Sepal.Length 150 5.843333 0.8280661   5.80 4.3 7.9  5.1  6.4  0.3086407
## Sepal.Width  150 3.057333 0.4358663   3.00 2.0 4.4  2.8  3.3  0.3126147
## Petal.Length 150 3.758000 1.7652982   4.35 1.0 6.9  1.6  5.1 -0.2694109
## Petal.Width  150 1.199333 0.7622377   1.30 0.1 2.5  0.3  1.8 -0.1009166
##                Kurtosis
## Sepal.Length -0.6058125
## Sepal.Width   0.1387047
## Petal.Length -1.4168574
## Petal.Width  -1.3581792
mvn(data = X, mvnTest = "royston")
## $multivariateNormality
##      Test        H      p value MVN
## 1 Royston 50.39667 3.098229e-11  NO
## 
## $univariateNormality
##               Test     Variable Statistic   p value Normality
## 1 Anderson-Darling Sepal.Length    0.8892  0.0225      NO    
## 2 Anderson-Darling Sepal.Width     0.9080  0.0202      NO    
## 3 Anderson-Darling Petal.Length    7.6785  <0.001      NO    
## 4 Anderson-Darling Petal.Width     5.1057  <0.001      NO    
## 
## $Descriptives
##                n     Mean   Std.Dev Median Min Max 25th 75th       Skew
## Sepal.Length 150 5.843333 0.8280661   5.80 4.3 7.9  5.1  6.4  0.3086407
## Sepal.Width  150 3.057333 0.4358663   3.00 2.0 4.4  2.8  3.3  0.3126147
## Petal.Length 150 3.758000 1.7652982   4.35 1.0 6.9  1.6  5.1 -0.2694109
## Petal.Width  150 1.199333 0.7622377   1.30 0.1 2.5  0.3  1.8 -0.1009166
##                Kurtosis
## Sepal.Length -0.6058125
## Sepal.Width   0.1387047
## Petal.Length -1.4168574
## Petal.Width  -1.3581792

Introducción

El análisis discriminante es una técnica de clasificación supervisada utilizada cuando las observaciones pertenecen a grupos conocidos y se desea construir una regla para clasificar nuevas observaciones.

Supongamos que existen \(K\) poblaciones

\[ \pi_1, \pi_2, \dots, \pi_K \]

y cada observación es un vector

\[ X = (X_1,\dots,X_p)^T \]

El objetivo es asignar una nueva observación al grupo más probable.

Supuestos del modelo

En el análisis discriminante lineal (LDA) se asume que

\[ X|\pi_k \sim N_p(\mu_k,\Sigma) \]

donde

Regla de clasificación

La probabilidad posterior es

\[ P(\pi_k|x) \]

y la observación se clasifica en el grupo que maximiza

\[ \delta_k(x) \]

donde

\[ \delta_k(x)=x^T\Sigma^{-1}\mu_k -\frac12\mu_k^T\Sigma^{-1}\mu_k +\log(\pi_k) \]

Esta función se llama función discriminante lineal.

Funciones discriminantes

Para dos grupos el problema se reduce a una función lineal

\[ w^Tx \]

donde

\[ w = \Sigma^{-1}(\mu_1-\mu_2) \]

La decisión depende del signo de esta función.

Ejemplo con el dataset Iris

library(MASS)
## Warning: package 'MASS' was built under R version 4.4.3
library(ggplot2)

data(iris)

head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
lda_model <- lda(Species ~ ., data = iris)

lda_model
## Call:
## lda(Species ~ ., data = iris)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## virginica         6.588       2.974        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8293776 -0.02410215
## Sepal.Width   1.5344731 -2.16452123
## Petal.Length -2.2012117  0.93192121
## Petal.Width  -2.8104603 -2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088
pred <- predict(lda_model)

head(pred$class)
## [1] setosa setosa setosa setosa setosa setosa
## Levels: setosa versicolor virginica
table(Predicted = pred$class, Real = iris$Species)
##             Real
## Predicted    setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         1
##   virginica       0          2        49
lda_values <- as.data.frame(pred$x)
lda_values$Species <- iris$Species

ggplot(lda_values, aes(LD1, LD2, color = Species)) +
  geom_point(size = 3) +
  labs(
    title = "Proyección discriminante LDA",
    x = "LD1",
    y = "LD2"
  ) +
  theme_minimal()