35  Análisis Exploratorio Multivariado: De la Estadística Descriptiva al Espacio Multidimensional

1 El problema: cuando mirar variable por variable no es suficiente

A lo largo de los capítulos anteriores de esta Parte III exploramos los pilares del análisis descriptivo: situamos el centro de una variable con medidas de localización (Sección 35.2), cuantificamos su dispersión (Sección 35.3), evaluamos la asimetría y el peso de sus colas (Sección 35.4), medimos la diversidad y la presencia de casos atípicos (Sección 35.5) y analizamos si dos características avanzan juntas mediante correlación (Sección 35.6).

Sin embargo, en la práctica de la minería de datos las observaciones rara vez viven en aislamiento. Una fila de una tabla de clientes, una transacción financiera o una medición de sensores industriales no está formada por una o dos variables desvinculadas, sino por un vector conjunto de \(p\) dimensiones:

\[ \mathbf{x}_i = (x_{i1}, x_{i2}, \dots, x_{ip})^T \in \mathbb{R}^p. \]

El error más común al enfrentarse a este escenario consiste en aplicar lo que podríamos llamar la falacia de las proyecciones marginales: descomponer el conjunto de datos en \(p\) columnas separadas, resumir cada una de forma univariada (media, desviación estándar, histograma) y asumir que la suma de esas descripciones captura la estructura global de los datos.

AdvertenciaLa falacia de las proyecciones marginales

Examinar cada variable por separado equivale a observar únicamente las sombras que un objeto tridimensional proyecta sobre las paredes de una habitación. Dos figuras con siluetas idénticas en dos planos ortogonales pueden ser, en el espacio tridimensional, objetos completamente distintos. En minería de datos, los fenómenos más críticos —agrupamientos, anomalías, cuellos de botella y dependencias— habitan en las relaciones conjuntas, no en las proyecciones aisladas.

Considere el siguiente caso intuitivo: un departamento de recursos humanos analiza a cinco candidatos mediante dos pruebas estandarizadas con escalas de 0 a 100: una de habilidad analítica (\(X_1\)) y otra de habilidad comunicativa (\(X_2\)). Históricamente, ambas habilidades presentan una correlación fuertemente positiva en la organización.

Candidato	Analítica (\(X_1\))	Comunicación (\(X_2\))	Diagnóstico univariado	Diagnóstico multivariado
A	50	52	Normal (cerca de la media)	Normal y coherente con la tendencia
B	85	88	Alto en ambas	Destacado balanceado
C	20	22	Bajo en ambas	Requiere capacitación
D	95	92	Sobresaliente	Sobresaliente balanceado
E	85	15	\(X_1\) normal-alto, \(X_2\) aceptable	¡Anomalía multivariada severa!

Si examinamos al candidato E variable por variable: 1. Su calificación en Analítica (85) es perfectamente plausible y coincide con la del candidato B. 2. Su calificación en Comunicación (15), aunque baja, se encuentra dentro de los rangos de la distribución marginal de la prueba.

Ningún filtro univariado (\(z\)-score individual o rango intercuartílico de Tukey) señalaría al candidato E como un caso extremo. No obstante, en el plano conjunto \((X_1, X_2)\), el candidato E rompe de manera radical la estructura de covarianza de los datos: posee un perfil extraordinariamente disonante.

NotaTres brújulas para el espacio multivariado

1. La localización no es solo un punto escalar: el centro de una nube en \(\mathbb{R}^p\) debe considerar cómo se compensan las coordenadas y si el estimador es invariante ante rotaciones del sistema de referencia.
2. La dispersión tiene orientación y forma: en \(\mathbb{R}^p\) la variabilidad no es un radio esférico, sino un elipsoide cuya orientación y excentricidad están determinadas por las covarianzas.
3. La distancia depende de la geometría estadística: medir separación en \(\mathbb{R}^p\) con la distancia euclidiana convencional ignora las correlaciones y produce juicios erróneos sobre la similitud y las anomalías.

---

2 Localización multivariada: el centro de una nube en \(\mathbb{R}^p\)

En el capítulo de localización univariada aprendimos que la media minimiza la suma de desviaciones cuadráticas, mientras que la mediana minimiza la suma de distancias absolutas y resiste valores extremos. ¿Cómo se extienden estas propiedades a una matriz de datos \(\mathbf{X} \in \mathbb{R}^{n \times p}\)?

2.1 El vector de medias muestral

El estimador clásico de localización multivariada es el vector de medias \(\bar{\mathbf{x}} \in \mathbb{R}^p\), obtenido calculando el promedio aritmético de cada columna:¹

¹ El vector de medias \(\bar{\mathbf{x}}\) minimiza la suma de distancias cuadráticas euclidianas a todos los puntos: \(\bar{\mathbf{x}} = \arg\min_{\mathbf{m} \in \mathbb{R}^p} \sum_{i=1}^n \|\mathbf{x}_i - \mathbf{m}\|_2^2\).

\[ \bar{\mathbf{x}} = \frac{1}{n} \sum_{i=1}^n \mathbf{x}_i = \begin{pmatrix} \bar{x}_1 \\ \bar{x}_2 \\ \vdots \\ \bar{x}_p \end{pmatrix}, \quad \text{donde } \bar{x}_j = \frac{1}{n}\sum_{i=1}^n x_{ij}. \]

Geométricamente, \(\bar{\mathbf{x}}\) es el centro de gravedad o centroide de la nube de puntos. Posee una propiedad fundamental: es un estimador lineal y conmuta con cualquier transformación afín de los datos.² Si definimos una transformación \(\mathbf{y}_i = \mathbf{A}\mathbf{x}_i + \mathbf{b}\), el nuevo centro de gravedad es simplemente \(\bar{\mathbf{y}} = \mathbf{A}\bar{\mathbf{x}} + \mathbf{b}\).

² Un estimador de localización \(\mathbf{T}(\mathbf{X})\) es equivariante afín si para cualquier matriz invertible \(\mathbf{A} \in \mathbb{R}^{p \times p}\) y cualquier vector de traslación \(\mathbf{b} \in \mathbb{R}^p\), se cumple \(\mathbf{T}(\mathbf{X}\mathbf{A}^T + \mathbf{1}\mathbf{b}^T) = \mathbf{A}\mathbf{T}(\mathbf{X}) + \mathbf{b}\).

Sin embargo, al igual que su contraparte univariada, el vector de medias hereda un punto de ruptura del 0% (Huber y Ronchetti 2009): una sola observación corrupta o errónea en cualquiera de las variables puede desplazar arbitrariamente todo el vector de medias en el espacio.

2.2 La falacia de la mediana por componentes

Ante la sensibilidad de la media, una tentación inmediata consiste en definir una “mediana multivariada” tomando la mediana univariada de cada coordenada de forma independiente:

\[ \mathbf{m}_{\text{marg}} = \begin{pmatrix} \operatorname{med}(x_{\cdot 1}) \\ \operatorname{med}(x_{\cdot 2}) \\ \vdots \\ \operatorname{med}(x_{\cdot p}) \end{pmatrix}. \]

Aunque sencilla de calcular, la mediana por componentes presenta dos deficiencias teóricas y prácticas graves (Small 1990):

1. No es equivariante por rotación: si rotamos la nube de puntos en el plano (por ejemplo, cambiando la orientación de los sensores o rotando los ejes), la mediana por componentes calculada sobre los nuevos ejes no coincide con la rotación de la mediana original.
2. Puede caer fuera de los datos y de su envolvente convexa: en configuraciones asimétricas o no convexas, el punto \(\mathbf{m}_{\text{marg}}\) puede ubicarse en una región donde la probabilidad de observar datos es nula.

# Demostración interactiva: Centroide vs Mediana por componentes vs Mediana geométrica
np.random.seed(42)

# Tres puntos en L
P = np.array([
    [1.0, 1.0],
    [1.0, 5.0],
    [5.0, 1.0]
])

# Vector de medias
mean_p = P.mean(axis=0)

# Mediana por componentes
med_marg = np.median(P, axis=0)

fig_loc = go.Figure()

# Conectar puntos formando el triángulo (casco convexo)
fig_loc.add_trace(go.Scatter(
    x=[P[0,0], P[1,0], P[2,0], P[0,0]],
    y=[P[0,1], P[1,1], P[2,1], P[0,1]],
    mode='lines',
    line=dict(color='rgba(150,150,150,0.5)', dash='dash'),
    name='Envolvente convexa',
    showlegend=True
))

# Puntos originales
fig_loc.add_trace(go.Scatter(
    x=P[:,0], y=P[:,1],
    mode='markers+text',
    text=['P1 (1,1)', 'P2 (1,5)', 'P3 (5,1)'],
    textposition=['bottom left', 'top left', 'bottom right'],
    marker=dict(size=14, color='royalblue'),
    name='Observaciones'
))

# Vector de medias
fig_loc.add_trace(go.Scatter(
    x=[mean_p[0]], y=[mean_p[1]],
    mode='markers+text',
    text=['Vector de Medias (2.33, 2.33)'],
    textposition='top right',
    marker=dict(size=14, color='crimson', symbol='star'),
    name='Centroide / Media'
))

# Mediana por componentes
fig_loc.add_trace(go.Scatter(
    x=[med_marg[0]], y=[med_marg[1]],
    mode='markers+text',
    text=['Mediana Marginal (1.0, 1.0) = P1'],
    textposition='top center',
    marker=dict(size=14, color='darkorange', symbol='diamond'),
    name='Mediana por componentes'
))

fig_loc.update_layout(
    title="Localización en ℝ²: Comparación de estimadores de centro",
    xaxis=dict(range=[0, 6], title="Variable 1"),
    yaxis=dict(range=[0, 6], title="Variable 2"),
    height=450,
    template="plotly_white"
)
fig_loc.show()

Como se observa en la figura anterior, la mediana por componentes colapsa arbitrariamente en el vértice \((1,1)\), mientras que el vector de medias se sitúa en el centro de gravedad interior del triángulo.

2.3 Mediana geométrica (Fermat–Weber) y profundidad estadística

Para obtener un centro verdaderamente multivariado y robusto, generalizamos el problema de Fermat-Weber: buscamos el punto \(\mathbf{m}_{\text{geom}} \in \mathbb{R}^p\) que minimice la suma de las distancias euclidianas directas a todas las observaciones (Small 1990):

\[ \mathbf{m}_{\text{geom}} = \arg\min_{\mathbf{m} \in \mathbb{R}^p} \sum_{i=1}^n \|\mathbf{x}_i - \mathbf{m}\|_2 = \arg\min_{\mathbf{m} \in \mathbb{R}^p} \sum_{i=1}^n \sqrt{\sum_{j=1}^p (x_{ij} - m_j)^2}. \]

A diferencia de la mediana por componentes: - La mediana geométrica (o mediana espacial) es equivariante ante rotaciones, traslaciones y cambios de escala isotrópicos. - Posee un punto de ruptura del 50%: se necesita corromper más de la mitad de las observaciones para moverla al infinito. - Aunque no tiene solución analítica cerrada en forma de ecuación simple, se aproxima de manera extraordinariamente rápida mediante el algoritmo iterativo de Weiszfeld.³

³ El algoritmo de Weiszfeld (1937) actualiza iterativamente la aproximación mediante \(\mathbf{m}^{(k+1)} = \left(\sum_{i=1}^n \frac{\mathbf{x}_i}{\|\mathbf{x}_i - \mathbf{m}^{(k)}\|}\right) / \left(\sum_{i=1}^n \frac{1}{\|\mathbf{x}_i - \mathbf{m}^{(k)}\|}\right)\).

Otra perspectiva fundamental proviene del concepto de profundidad estadística (halfspace depth o profundidad de Tukey) (Tukey 1975). La profundidad de un punto \(\mathbf{\theta}\) respecto a una nube de datos mide la fracción mínima de observaciones contenida en cualquier semiespacio cerrado que contenga a \(\mathbf{\theta}\). El punto de máxima profundidad de Tukey representa la mediana de máxima centralidad en \(\mathbb{R}^p\).

Medida	Definición operativa	Equivarianza por rotación	Punto de ruptura	Costo computacional
Vector de medias	Promedio aritmético por columna	Sí (afín)	0%	Muy bajo \(\mathcal{O}(np)\)
Mediana marginal	Mediana independiente por columna	No	50%	Bajo \(\mathcal{O}(np \log n)\)
Mediana geométrica	Minimiza suma de distancias euclidianas	Sí (ortogonal)	50%	Moderado (iterativo)
Medoid multivariado	Observación real que minimiza distancias	Depende de la distancia	Robusto	\(\mathcal{O}(n^2 p)\)

---

3 Variabilidad y dispersión multivariada: forma y orientación

En el análisis univariado, la dispersión se resume en un único número no negativo: la varianza \(s^2\) o la desviación estándar \(s\). Si una variable tiene \(s=10\), sabemos exactamente cuánto se alejan los datos de su promedio sobre esa recta.

En \(\mathbb{R}^p\), sin embargo, una nube de puntos puede tener la misma dispersión en cada coordenada por separado y aun así presentar estructuras geométricas radicalmente distintas. La dispersión multivariada no es un escalar; es un tensor geométrico que define tamaño, forma y orientación.

3.1 La matriz de varianza-covarianza

La estructura de dispersión conjunta se captura en la matriz de varianza-covarianza muestral \(\mathbf{S} \in \mathbb{R}^{p \times p}\):⁴

⁴ El divisor \(n-1\) (corrección de Bessel) asegura que \(\mathbf{S}\) sea un estimador insesgado de la matriz de covarianza poblacional \(\boldsymbol{\Sigma}\): \(\mathbb{E}[\mathbf{S}] = \boldsymbol{\Sigma}\).

\[ \mathbf{S} = \frac{1}{n-1} \sum_{i=1}^n (\mathbf{x}_i - \bar{\mathbf{x}})(\mathbf{x}_i - \bar{\mathbf{x}})^T = \begin{pmatrix} s_{11} & s_{12} & \cdots & s_{1p} \\ s_{21} & s_{22} & \cdots & s_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ s_{p1} & s_{p2} & \cdots & s_{pp} \end{pmatrix}. \]

Propiedades esenciales de \(\mathbf{S}\): 1. Diagonal principal: \(s_{jj} = s_j^2\) contiene la varianza muestral de cada variable individual. 2. Elementos fuera de la diagonal: \(s_{jk} = s_{kj}\) contiene la covarianza entre la variable \(j\) y la variable \(k\). 3. Simetría y positividad: \(\mathbf{S}\) es una matriz simétrica (\(s_{jk} = s_{kj}\)) y semidefinida positiva (\(\mathbf{a}^T \mathbf{S} \mathbf{a} \ge 0\) para cualquier vector \(\mathbf{a} \in \mathbb{R}^p\)). Si no existen variables linealmente redundantes y \(n > p\), \(\mathbf{S}\) es estrictamente definida positiva.

3.2 Resúmenes escalares: Varianza Total vs. Varianza Generalizada

Cuando necesitamos comparar la dispersión global entre dos grupos o monitorear la variabilidad total de un sistema, reducir la matriz \(\mathbf{S}\) a un único escalar resulta conveniente. Existen dos alternativas clásicas con significados geométricos muy diferentes:

1. Varianza Total (Traza de \(\mathbf{S}\)).

La varianza total se define como la traza de la matriz de covarianza:⁵

⁵ La traza de una matriz cuadrada es la suma de sus elementos diagonales y coincide con la suma de sus autovalores: \(\operatorname{tr}(\mathbf{S}) = \sum_{j=1}^p s_{jj} = \sum_{j=1}^p \lambda_j\).

\[ \operatorname{Var}_{\text{total}}(\mathbf{X}) = \operatorname{tr}(\mathbf{S}) = \sum_{j=1}^p s_{jj} = \sum_{j=1}^p \lambda_j, \]

donde \(\lambda_1, \dots, \lambda_p\) son los autovalores de \(\mathbf{S}\).

AdvertenciaLa ceguera de la Varianza Total

La varianza total es simplemente la suma de las varianzas univariadas. Ignora por completo las covarianzas \(s_{jk}\). Dos conjuntos de datos pueden tener exactamente la misma varianza total, pero uno ser una esfera isotrópica no correlacionada y el otro ser una línea delgada y altamente correlacionada.

2. Varianza Generalizada de Wilks (Determinante de \(\mathbf{S}\)).

Introducida por Samuel Wilks (Wilks 1932), la varianza generalizada se define como el determinante de la matriz de covarianza:

\[ \operatorname{Var}_{\text{gen}}(\mathbf{X}) = |\mathbf{S}| = \det(\mathbf{S}) = \prod_{j=1}^p \lambda_j. \]

Geométricamente, el determinante \(|\mathbf{S}|\) es proporcional al cuadrado del hipervolumen del elipsoide de dispersión generado por los datos (Johnson y Wichern 2007).

• Si las variables son independientes, \(|\mathbf{S}| = s_{11} s_{22} \cdots s_{pp}\).
• A medida que la correlación entre variables aumenta, el elipsoide se aplana y \(|\mathbf{S}|\) disminuye.
• Si existe una relación lineal exacta entre variables (colinealidad perfecta), el elipsoide colapsa a una dimensión inferior y \(|\mathbf{S}| = 0\), revelando que los datos no ocupan un volumen real en \(\mathbb{R}^p\).

# Comparación visual: Varianza Total vs Varianza Generalizada
np.random.seed(101)
n_samples = 400

# Caso A: Sin correlación (círculo)
cov_a = np.array([[2.0, 0.0], [0.0, 2.0]])
data_a = np.random.multivariate_normal([0, 0], cov_a, n_samples)

# Caso B: Fuerte correlación positiva (elipse alargada)
cov_b = np.array([[2.0, 1.8], [1.8, 2.0]])
data_b = np.random.multivariate_normal([0, 0], cov_b, n_samples)

tr_a, det_a = np.trace(cov_a), np.linalg.det(cov_a)
tr_b, det_b = np.trace(cov_b), np.linalg.det(cov_b)

fig_var = make_subplots(
    rows=1, cols=2,
    subplot_titles=[
        f"A: Sin correlación<br>Traza = {tr_a:.1f} | Det = {det_a:.2f}",
        f"B: Alta correlación (r = 0.9)<br>Traza = {tr_b:.1f} | Det = {det_b:.2f}"
    ]
)

fig_var.add_trace(go.Scatter(
    x=data_a[:,0], y=data_a[:,1], mode='markers',
    marker=dict(size=4, color='royalblue', opacity=0.6),
    name='Datos A'
), row=1, col=1)

fig_var.add_trace(go.Scatter(
    x=data_b[:,0], y=data_b[:,1], mode='markers',
    marker=dict(size=4, color='crimson', opacity=0.6),
    name='Datos B'
), row=1, col=2)

fig_var.update_layout(
    height=420,
    showlegend=False,
    template="plotly_white",
    title="Misma Varianza Total (Traza = 4.0), muy distinta Varianza Generalizada (Det)"
)
fig_var.update_xaxes(range=[-6, 6])
fig_var.update_yaxes(range=[-6, 6])
fig_var.show()

3.3 De la distancia euclidiana a la distancia de Mahalanobis

¿Qué tan lejos está una observación \(\mathbf{x}\) del centro de la población \(\bar{\mathbf{x}}\)?

La respuesta natural en geometría elemental es la distancia euclidiana:

\[ d_E(\mathbf{x}, \bar{\mathbf{x}}) = \sqrt{(\mathbf{x} - \bar{\mathbf{x}})^T (\mathbf{x} - \bar{\mathbf{x}})} = \sqrt{\sum_{j=1}^p (x_j - \bar{x}_j)^2}. \]

Sin embargo, en estadística y minería de datos, la distancia euclidiana es profundamente engañosa por dos motivos: 1. Sensibilidad a las escalas de medida: si \(X_1\) se mide en metros y \(X_2\) en milímetros, \(X_2\) dominará numéricamente la distancia euclidiana por un factor de un millón. 2. Ceguera ante la correlación: asume que todas las direcciones del espacio son equiprobables (isocuantas esféricas). Si dos variables están correlacionadas positivamente, desplazarse a lo largo del eje principal de correlación es habitual, mientras que desplazarse perpendicularmente a la correlación es extraordinariamente improbable.

La solución matemática rigurosa es la distancia estadística o distancia de Mahalanobis (Mahalanobis 1936):

\[ D_M(\mathbf{x}, \bar{\mathbf{x}}) = \sqrt{(\mathbf{x} - \bar{\mathbf{x}})^T \mathbf{S}^{-1} (\mathbf{x} - \bar{\mathbf{x}})}. \]

Descomposición intuitiva en tres etapas.

Para entender por qué la distancia de Mahalanobis corrige todas las distorsiones, descompongamos \(\mathbf{S}^{-1}\) mediante su descomposición espectral \(\mathbf{S} = \mathbf{V} \mathbf{\Lambda} \mathbf{V}^T\) (donde \(\mathbf{V}\) son los autovectores y \(\mathbf{\Lambda}\) los autovalores):

1. Centrado: \((\mathbf{x} - \bar{\mathbf{x}})\) traslada el origen al centro de la nube.
2. Rotación ortogonal: \(\mathbf{V}^T(\mathbf{x} - \bar{\mathbf{x}})\) rota los ejes hacia las direcciones principales de variación natural de los datos, eliminando las correlaciones.
3. Estandarización por escala propia: multiplicar por \(\mathbf{\Lambda}^{-1/2}\) divide cada nuevo eje rotado entre su propia desviación estándar \(\sqrt{\lambda_j}\).

Por tanto, la distancia de Mahalanobis es simplemente la distancia euclidiana calculada después de descorrelacionar y estandarizar la nube de datos.⁶

⁶ Si \(\mathbf{S} = \mathbf{V}\mathbf{\Lambda}\mathbf{V}^T\), la transformación \(\mathbf{z}_i = \mathbf{\Lambda}^{-1/2}\mathbf{V}^T(\mathbf{x}_i - \bar{\mathbf{x}})\) genera variables incorrelacionadas con matriz de covarianza \(\mathbf{I}_p\), donde \(D_M(\mathbf{x}_i, \bar{\mathbf{x}}) = \|\mathbf{z}_i\|_2\).

# Visualización interactiva de curvas de nivel: Euclidiana vs Mahalanobis
np.random.seed(42)
mean_ex = np.array([0, 0])
cov_ex = np.array([[4.0, 3.2], [3.2, 4.0]]) # correlación r = 0.8

# Puntos de prueba a igual distancia euclidiana del centro
r_dist = 3.0
angles = np.array([np.pi/4, 3*np.pi/4]) # P1 a lo largo de la correlación, P2 ortogonal
test_pts = np.array([
    [r_dist * np.cos(angles[0]), r_dist * np.sin(angles[0])], # (2.12, 2.12)
    [r_dist * np.cos(angles[1]), r_dist * np.sin(angles[1])]  # (-2.12, 2.12)
])

inv_cov = np.linalg.inv(cov_ex)
d_m1 = np.sqrt(test_pts[0] @ inv_cov @ test_pts[0])
d_m2 = np.sqrt(test_pts[1] @ inv_cov @ test_pts[1])

# Malla para elipses de Mahalanobis
theta_grid = np.linspace(0, 2*np.pi, 200)
circle = np.array([np.cos(theta_grid), np.sin(theta_grid)])

# Autodescomposición de cov_ex
vals, vecs = np.linalg.eigh(cov_ex)
# Elipse D_M = 1 y D_M = 2
ellipse_1 = vecs @ np.diag(np.sqrt(vals) * 1.0) @ circle
ellipse_2 = vecs @ np.diag(np.sqrt(vals) * 2.0) @ circle

fig_mah = go.Figure()

# Nube de fondo simulada
bg_data = np.random.multivariate_normal(mean_ex, cov_ex, 300)
fig_mah.add_trace(go.Scatter(
    x=bg_data[:,0], y=bg_data[:,1], mode='markers',
    marker=dict(size=3, color='lightgray'),
    name='Distribución de fondo'
))

# Elipses de Mahalanobis
fig_mah.add_trace(go.Scatter(
    x=ellipse_1[0], y=ellipse_1[1], mode='lines',
    line=dict(color='crimson', dash='dash'), name='Contorno D_M = 1'
))
fig_mah.add_trace(go.Scatter(
    x=ellipse_2[0], y=ellipse_2[1], mode='lines',
    line=dict(color='crimson', width=2), name='Contorno D_M = 2'
))

# Círculo euclidiano de radio 3
fig_mah.add_trace(go.Scatter(
    x=r_dist * np.cos(theta_grid), y=r_dist * np.sin(theta_grid), mode='lines',
    line=dict(color='royalblue', width=1.5, dash='dot'), name='Círculo Euclidiano (d_E = 3)'
))

# Puntos de prueba
fig_mah.add_trace(go.Scatter(
    x=[test_pts[0,0]], y=[test_pts[0,1]], mode='markers+text',
    text=[f'P1: d_E=3.0, D_M={d_m1:.2f}'],
    textposition='top right',
    marker=dict(size=12, color='green', symbol='circle'),
    name='Punto 1 (esperado)'
))

fig_mah.add_trace(go.Scatter(
    x=[test_pts[1,0]], y=[test_pts[1,1]], mode='markers+text',
    text=[f'P2: d_E=3.0, D_M={d_m2:.2f}'],
    textposition='top left',
    marker=dict(size=12, color='purple', symbol='x'),
    name='Punto 2 (atípico)'
))

fig_mah.update_layout(
    title="Distancia Euclidiana vs. Distancia de Mahalanobis",
    xaxis=dict(range=[-5, 5], title="Variable 1"),
    yaxis=dict(range=[-5, 5], title="Variable 2"),
    height=480,
    template="plotly_white"
)
fig_mah.show()

Observe la conclusión: ambos puntos están a exactamente 3 unidades euclidianas del centro. Sin embargo, para la estructura estadística, el punto verde \(P_1\) tiene \(D_M = 1.12\) (es un dato completamente común y predecible), mientras que el punto morado \(P_2\) tiene \(D_M = 3.35\) (se ubica en una región de bajísima densidad y constituye una anomalía clara).

---

4 Forma y normalidad multivariada

En el capítulo de medidas de forma (Sección 35.4) estudiamos la asimetría (\(g_1\)) y la curtosis (\(g_2\)) para caracterizar cómo se desvía una distribución de la curva gaussiana univariada (Joanes y Gill 1998).

Al pasar a \(\mathbb{R}^p\), la distribución de referencia por excelencia es la distribución normal multivariada \(\mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), cuya función de densidad es:⁷

⁷ La distribución normal multivariada se denota \(\mathbf{X} \sim \mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\). Si \(\boldsymbol{\Sigma}\) no es invertible, la distribución reside en un subespacio afín de dimensión inferior y se denomina normal singular.

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

Nótese que el exponente contiene precisamente el cuadrado de la distancia de Mahalanobis: \(-\frac{1}{2} D_M^2(\mathbf{x}, \boldsymbol{\mu})\). Por esta razón, las curvas o superficies de nivel donde la densidad \(f(\mathbf{x})\) es constante corresponden exactamente a hiperelipsoides de Mahalanobis.

4.1 La falacia: ¿marginales normales implican normalidad conjunta?

Uno de los mitos más arraigados en el análisis exploratorio es asumir que si cada variable por separado supera una prueba de normalidad univariada (como Shapiro-Wilk o Kolmogorov-Smirnov), entonces el conjunto completo es normal multivariado.

PrecauciónAdvertencia teórica fundamental

La normalidad de todas las distribuciones marginales es una condición necesaria, pero NO suficiente para la normalidad multivariada. Es perfectamente posible construir distribuciones donde cada coordenada individual sea una normal estándar \(\mathcal{N}(0,1)\), pero cuya distribución conjunta tenga huecos, formas no elípticas o dependencias no lineales severas (Johnson y Wichern 2007).

Considere el siguiente contraejemplo clásico: sea \(X \sim \mathcal{N}(0, 1)\), y definamos \(Y\) como:

\[ Y = \begin{cases} X & \text{si } |X| < c \\ -X & \text{si } |X| \ge c \end{cases} \]

para una constante \(c > 0\). Marginalmente, tanto \(X\) como \(Y\) tienen distribuciones exactamente normales estándar \(\mathcal{N}(0,1)\). Sin embargo, su distribución conjunta dista radicalmente de ser una normal bivariada: no tiene curvas de nivel elípticas y contiene franjas vacías.

# Contraejemplo interactivo: marginales normales pero conjunto no normal
np.random.seed(42)
n_pts = 1500
x_norm = np.random.normal(0, 1, n_pts)
c_val = 1.0
y_non_norm = np.where(np.abs(x_norm) < c_val, x_norm, -x_norm)

fig_contra = make_subplots(
    rows=1, cols=3,
    subplot_titles=[
        "Marginal X (Normal)",
        "Marginal Y (Normal)",
        "Conjunta (X, Y) — ¡NO es normal!"
    ]
)

fig_contra.add_trace(go.Histogram(x=x_norm, nbinsx=30, marker_color='royalblue', name='X'), row=1, col=1)
fig_contra.add_trace(go.Histogram(x=y_non_norm, nbinsx=30, marker_color='darkorange', name='Y'), row=1, col=2)
fig_contra.add_trace(go.Scatter(x=x_norm, y=y_non_norm, mode='markers', marker=dict(size=4, color='crimson', opacity=0.5), name='(X,Y)'), row=1, col=3)

fig_contra.update_layout(
    height=380,
    showlegend=False,
    template="plotly_white",
    title="Contraejemplo: Marginales gaussianas que no forman una normal multivariada"
)
fig_contra.show()

4.2 Asimetría y curtosis multivariadas de Mardia

Para contrastar rigurosamente si una muestra \(p\)-dimensional proviene de una población normal multivariada, K. V. Mardia (Mardia 1970) introdujo extensiones algebraicas basadas en la distancia de Mahalanobis cruzada entre pares de observaciones:

Sean \(g_{ij} = (\mathbf{x}_i - \bar{\mathbf{x}})^T \mathbf{S}^{-1} (\mathbf{x}_j - \bar{\mathbf{x}})\).

1. Asimetría multivariada de Mardia (\(b_{1,p}\)): \[ b_{1,p} = \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n g_{ij}^3. \] Bajo la hipótesis nula de normalidad multivariada, el estadístico \(\frac{n}{6} b_{1,p}\) sigue asintóticamente una distribución Chi-cuadrado con \(\frac{p(p+1)(p+2)}{6}\) grados de libertad.

2. Curtosis multivariada de Mardia (\(b_{2,p}\)): \[ b_{2,p} = \frac{1}{n} \sum_{i=1}^n g_{ii}^2 = \frac{1}{n} \sum_{i=1}^n \left[ D_M^2(\mathbf{x}_i, \bar{\mathbf{x}}) \right]^2. \] Bajo normalidad multivariada, la esperanza teórica es \(\mathbb{E}[b_{2,p}] = p(p+2)\). Valores significativamente mayores indican colas pesadas multivariadas (leptocurtosis conjunta).

---

5 Heterogeneidad y anomalías multivariadas

En el capítulo de heterogeneidad (Sección 35.5) estudiamos cómo identificar valores atípicos mediante el rango intercuartílico (\(1.5 \times \text{IQR}\)) o umbrales de \(z\)-score (\(\pm 3\sigma\)). En el espacio multivariado, estas herramientas univariadas fallan por dos razones fundamentales:

1. El fenómeno del outlier oculto: un punto puede estar dentro del rango normal de cada variable individual y aun así ser un valor atípico extremo al considerar la correlación entre ellas.
2. La maldición de la dimensionalidad: a medida que \(p\) crece, el volumen del hipercubo que encierra a los datos crece exponencialmente respecto a la región ocupada por la distribución central, haciendo que las observaciones se vuelvan dispersas.

5.1 La distribución teórica de \(D_M^2\) y umbrales Chi-cuadrado

Si una población de \(p\) variables sigue una distribución normal multivariada \(\mathcal{N}_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})\), el cuadrado de la distancia de Mahalanobis sigue exactamente una distribución Chi-cuadrado con \(p\) grados de libertad:⁸

⁸ Esta propiedad se deriva de que \(D_M^2 = \mathbf{Z}^T \mathbf{Z} = \sum_{j=1}^p Z_j^2\), donde cada \(Z_j \sim \mathcal{N}(0, 1)\) de forma independiente tras la diagonalización ortogonal de \(\boldsymbol{\Sigma}\).

\[ D_M^2(\mathbf{X}, \boldsymbol{\mu}) = (\mathbf{X} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{X} - \boldsymbol{\mu}) \sim \chi^2_p. \]

Esto proporciona una regla de decisión probabilística inmediata para el análisis exploratorio: una observación \(\mathbf{x}_i\) se clasifica como un outlier multivariado potencial con un nivel de significancia \(\alpha\) (típicamente \(\alpha = 0.01\) o \(\alpha = 0.001\)) si su distancia supera el valor crítico:

\[ D_M^2(\mathbf{x}_i, \bar{\mathbf{x}}) > \chi^2_{p, 1-\alpha}. \]

Dimensión (\(p\))	Umbral al 95% (\(\chi^2_{p, 0.95}\))	Umbral al 97.5% (\(\chi^2_{p, 0.975}\))	Umbral al 99% (\(\chi^2_{p, 0.99}\))
2	5.99	7.38	9.21
3	7.81	9.35	11.34
5	11.07	12.83	15.09
10	18.31	20.48	23.21

5.2 Efecto máscara (masking) y pantallazo (swamping)

La regla del umbral \(\chi^2_p\) presenta una trampa sutil cuando se utiliza con los estimadores clásicos \((\bar{\mathbf{x}}, \mathbf{S})\). Como los propios outliers forman parte del cálculo de la media y la covarianza:

• Efecto Máscara (Masking): un grupo de outliers severos puede arrastrar a \(\bar{\mathbf{x}}\) hacia ellos e inflar artificialmente la matriz \(\mathbf{S}\). Al expandirse \(\mathbf{S}\), las distancias de Mahalanobis de los propios outliers se reducen, haciendo que pasen desapercibidos como si fueran puntos normales.
• Efecto Pantallazo (Swamping): la inflación de \(\mathbf{S}\) y el desplazamiento de \(\bar{\mathbf{x}}\) provocan que puntos perfectamente normales y limpios queden fuera de la nube central y sean erróneamente clasificados como anomalías.

5.3 Estimación robusta: Determinante de Covarianza Mínima (MCD)

Para evitar el enmascaramiento, Peter Rousseeuw (Rousseeuw 1984; Rousseeuw y Van Driessen 1999) introdujo el estimador del Determinante de Covarianza Mínima (Minimum Covariance Determinant o MCD).

La idea intuitiva es magistral: de las \(n\) observaciones totales, buscamos el subconjunto de tamaño \(h\) (donde \(\frac{n+p+1}{2} \le h < n\), comúnmente \(h \approx 0.75n\)) cuya matriz de covarianza tenga el menor determinante posible. 1. Al minimizar el determinante \(|\mathbf{S}_h|\), el algoritmo selecciona automáticamente el subconjunto más denso y libre de valores atípicos. 2. La media y covarianza se estiman exclusivamente sobre ese subconjunto limpio: \((\bar{\mathbf{x}}_{\text{MCD}}, \mathbf{S}_{\text{MCD}})\). 3. Las distancias de Mahalanobis robustas se calculan respecto a este centro no contaminado: \[ D_{\text{MCD}}(\mathbf{x}_i) = \sqrt{(\mathbf{x}_i - \bar{\mathbf{x}}_{\text{MCD}})^T \mathbf{S}_{\text{MCD}}^{-1} (\mathbf{x}_i - \bar{\mathbf{x}}_{\text{MCD}})}. \]

---

6 Dependencia multivariada y correlación parcial

En el capítulo de correlación (Sección 35.6) medimos la fuerza de la asociación bivariada con el coeficiente \(r\) de Pearson. Sin embargo, en un sistema con múltiples variables, una correlación observada entre \(X\) y \(Y\) puede ser espuria, producto de una tercera variable confusora \(Z\) que afecta a ambas.

6.1 La Paradoja de Simpson

El caso más extremo de distorsión al ignorar variables adicionales es la Paradoja de Simpson (Simpson 1951): un fenómeno estadístico donde una tendencia o correlación observada en varios subgrupos desaparece o se invierte completamente cuando los grupos se combinan en un análisis global.

# Visualización interactiva de la Paradoja de Simpson
np.random.seed(42)
n_grp = 100

# Grupo 1: Jóvenes (menor salario base, correlación positiva horas-productividad)
x1 = np.random.uniform(20, 35, n_grp)
y1 = 15 + 0.8 * x1 + np.random.normal(0, 2, n_grp)

# Grupo 2: Intermedios
x2 = np.random.uniform(30, 45, n_grp)
y2 = 30 + 0.8 * x2 + np.random.normal(0, 2, n_grp)

# Grupo 3: Seniors
x3 = np.random.uniform(40, 55, n_grp)
y3 = 45 + 0.8 * x3 + np.random.normal(0, 2, n_grp)

# Agregados
x_all = np.concatenate([x1, x2, x3])
y_all = np.concatenate([y1, y2, y3])
grp_labels = ['Grupo A (Jr)'] * n_grp + ['Grupo B (Mid)'] * n_grp + ['Grupo C (Sr)'] * n_grp

df_simp = pd.DataFrame({'Horas': x_all, 'Productividad': y_all, 'Nivel': grp_labels})

fig_simp = px.scatter(
    df_simp, x='Horas', y='Productividad', color='Nivel',
    color_discrete_sequence=['#1f77b4', '#ff7f0e', '#2ca02c'],
    trendline='ols',
    title="Paradoja de Simpson: Dentro de cada grupo r > 0, pero el análisis global r_total puede distorsionarse"
)
fig_simp.update_layout(height=420, template="plotly_white")
fig_simp.show()

6.2 Correlación Parcial

Para discernir si la relación entre dos variables \(X_j\) y \(X_k\) es genuina o está mediada por otras variables \(\mathbf{Z} = \{X_1, \dots, X_p\} \setminus \{X_j, X_k\}\), calculamos la correlación parcial \(\rho_{jk \cdot \mathbf{Z}}\).

La correlación parcial mide la correlación de Pearson entre los residuos resultantes de haber ajustado una regresión lineal de \(X_j\) sobre \(\mathbf{Z}\) y de \(X_k\) sobre \(\mathbf{Z}\). Para el caso de tres variables (\(X_1, X_2\) condicionado a \(X_3\)):

\[ r_{12 \cdot 3} = \frac{r_{12} - r_{13} r_{23}}{\sqrt{(1 - r_{13}^2)(1 - r_{23}^2)}}. \]

6.3 La Matriz de Precisión \(\boldsymbol{\Theta} = \mathbf{S}^{-1}\)

En dimensiones mayores (\(p > 3\)), calcular todas las regresiones residuales por separado es ineficiente. La teoría multivariada proporciona una conexión matemática directa a través de la inversa de la matriz de covarianza, conocida como la matriz de precisión o matriz de concentración:⁹

⁹ La matriz de precisión \(\boldsymbol{\Theta} = \boldsymbol{\Sigma}^{-1}\) codifica la estructura de independencia condicional en modelos gráficos gaussianos: \(X_j \perp X_k \mid \mathbf{X}_{\setminus \{j,k\}} \iff \theta_{jk} = 0\).

\[ \boldsymbol{\Theta} = \mathbf{S}^{-1} = (\theta_{jk})_{p \times p}. \]

La correlación parcial entre \(X_j\) y \(X_k\) condicionada a todas las demás variables del conjunto de datos viene dada por:

\[ r_{jk \cdot \text{resto}} = -\frac{\theta_{jk}}{\sqrt{\theta_{jj} \theta_{kk}}}. \]

NotaInterpretación en redes gaussianas y minería de datos

Si la entrada \(\theta_{jk} = 0\), las variables \(X_j\) y \(X_k\) son condicionalmente independientes dado el resto de las variables. Esto significa que toda la correlación aparente entre \(X_j\) y \(X_k\) puede explicarse completamente a través de caminos intermedios formados por otras variables del sistema. Esta propiedad es el fundamento de los modelos gráficos probabilísticos y el descubrimiento causal en minería de datos.

---

7 Técnicas de visualización multivariada

A diferencia del plano bidimensional o el espacio tridimensional, los espacios de \(p \ge 4\) dimensiones no admiten una proyección directa natural en la retina humana. En minería de datos recurrimos a tres técnicas complementarias para inspeccionar la geometría conjunta:

7.1 Matriz de dispersión interactiva (SPLOM) con elipsoides de confianza

Permite inspeccionar simultáneamente todas las proyecciones bidimensionales \(\binom{p}{2}\), complementadas con diagramas de densidad univariada en la diagonal principal y elipsoides de Mahalanobis para evaluar colinealidades locales.

7.2 Coordenadas Paralelas (Parallel Coordinates)

En lugar de disponer los ejes de forma ortogonal (lo cual restringe la visualización a 3 dimensiones), los \(p\) ejes se disponen en paralelo verticalmente. Cada observación \(\mathbf{x}_i \in \mathbb{R}^p\) se convierte en una línea poligonal continua que cruza cada eje en su valor correspondiente. - Permite detectar agrupamientos (clusters) como haces densos de líneas que viajan juntas. - Correlaciones fuertemente positivas aparecen como líneas paralelas entre dos ejes. - Correlaciones negativas aparecen como líneas que se cruzan en forma de tijera entre dos ejes contiguos.

7.3 Diagramas de Radar (Spider Charts)

Útiles para comparar el perfil multivariado de centroides o de observaciones individuales específicas respecto a las medias globales.

---

8 Laboratorio práctico interactivo

A continuación, implementaremos un flujo de análisis exploratorio multivariado completo sobre un conjunto de datos realista de monitoreo ambiental y de sensores industriales, evaluando localización, matrices de covarianza, correlaciones parciales y detección robusta de anomalías.

# 1. Simulación de un dataset multivariado con estructura de dependencia y outliers
np.random.seed(2026)
n_inliers = 200
n_outliers = 10

# 4 variables: Temperatura, Presión, Vibración, Consumo eléctrico
mean_real = np.array([50.0, 100.0, 1.5, 30.0])
cov_real = np.array([
    [16.0,  9.6,  1.2,  8.0],   # Temp
    [ 9.6, 25.0,  2.0, 15.0],   # Presión (r_12 = 9.6 / (4*5) = 0.48)
    [ 1.2,  2.0,  0.36, 0.9],   # Vibración
    [ 8.0, 15.0,  0.9, 20.0]    # Consumo
])

# Generar datos normales correlacionados
inliers = np.random.multivariate_normal(mean_real, cov_real, n_inliers)

# Generar outliers multivariados patológicos:
# Puntos normales en sus marginales pero disonantes con la correlación conjunta
outliers = np.array([
    [52.0, 120.0, 0.4, 45.0],   # Alta presión pero vibración casi nula
    [40.0,  80.0, 3.2, 15.0],   # Baja temperatura pero vibración extrema
    [62.0, 115.0, 0.6, 18.0],   # Alta temp pero consumo anormalmente bajo
    [48.0,  95.0, 2.8, 48.0],
    [55.0, 102.0, 0.2, 52.0],
    [42.0, 118.0, 1.4, 20.0],
    [58.0,  82.0, 2.5, 42.0],
    [46.0, 105.0, 0.5, 16.0],
    [64.0,  88.0, 2.9, 22.0],
    [38.0, 112.0, 0.3, 38.0]
])

X_raw = np.vstack([inliers, outliers])
cols = ['Temperatura', 'Presión', 'Vibración', 'Consumo']
df_lab = pd.DataFrame(X_raw, columns=cols)
df_lab['Tipo_Real'] = ['Normal'] * n_inliers + ['Outlier'] * n_outliers

print(f"Dataset generado: {df_lab.shape[0]} observaciones, {df_lab.shape[1]-1} variables.")

Dataset generado: 210 observaciones, 4 variables.

8.1 Paso 1: Localización y dispersión multivariada

# Vector de medias
mean_vec = df_lab[cols].mean()

# Mediana geométrica aproximada mediante el algoritmo de Weiszfeld
def geometric_median(X, eps=1e-5, max_iter=200):
    y = np.mean(X, axis=0)
    for _ in range(max_iter):
        distances = np.linalg.norm(X - y, axis=1)
        distances = np.where(distances == 0, eps, distances)
        weights = 1.0 / distances
        y_next = np.sum(X * weights[:, np.newaxis], axis=0) / np.sum(weights)
        if np.linalg.norm(y_next - y) < eps:
            break
        y = y_next
    return y

geom_med = geometric_median(df_lab[cols].values)

# Matriz de covarianza y resúmenes escalares
S_mat = df_lab[cols].cov().values
var_total = np.trace(S_mat)
var_gen = np.linalg.det(S_mat)

df_resumen_loc = pd.DataFrame({
    'Variable': cols,
    'Media Muestral': mean_vec.values.round(3),
    'Mediana Geométrica': geom_med.round(3)
})

print(df_resumen_loc.to_string(index=False))
print(f"\nVarianza Total (Traza de S): {var_total:.3f}")
print(f"Varianza Generalizada de Wilks (|S|): {var_gen:.3f}")

   Variable  Media Muestral  Mediana Geométrica
Temperatura          49.837              49.834
    Presión          99.362              99.128
  Vibración           1.476               1.447
    Consumo          29.575              29.409

Varianza Total (Traza de S): 83.527
Varianza Generalizada de Wilks (|S|): 3840.341

8.2 Paso 2: Detección de anomalías con Distancia de Mahalanobis vs. Filtros Univariados

from sklearn.covariance import MinCovDet

X_mat = df_lab[cols].values
p = len(cols)

# Distancia de Mahalanobis clásica
inv_S = np.linalg.inv(S_mat)
diff_classic = X_mat - mean_vec.values
dm2_classic = np.sum(diff_classic @ inv_S * diff_classic, axis=1)

# Distancia de Mahalanobis robusta (MCD)
mcd = MinCovDet(random_state=42).fit(X_mat)
dm2_robust = mcd.mahalanobis(X_mat)

# Umbral teórico Chi-cuadrado al 99%
cutoff_chi2 = stats.chi2.ppf(0.99, df=p)

df_lab['DM2_Clasica'] = dm2_classic
df_lab['DM2_Robusta'] = dm2_robust
df_lab['Outlier_Clasico'] = df_lab['DM2_Clasica'] > cutoff_chi2
df_lab['Outlier_Robusto'] = df_lab['DM2_Robusta'] > cutoff_chi2

# Comparación de detección
print(f"Umbral Chi² (p={p}, alpha=0.01): {cutoff_chi2:.2f}")
print("\nMatriz de confusión detección Clásica vs Real:")
print(pd.crosstab(df_lab['Tipo_Real'], df_lab['Outlier_Clasico'], rownames=['Real'], colnames=['Clasico (> cutoff)']))

print("\nMatriz de confusión detección Robusta (MCD) vs Real:")
print(pd.crosstab(df_lab['Tipo_Real'], df_lab['Outlier_Robusto'], rownames=['Real'], colnames=['Robusto (> cutoff)']))

Umbral Chi² (p=4, alpha=0.01): 13.28

Matriz de confusión detección Clásica vs Real:
Clasico (> cutoff)  False  True 
Real                            
Normal                199      1
Outlier                 0     10

Matriz de confusión detección Robusta (MCD) vs Real:
Robusto (> cutoff)  False  True 
Real                            
Normal                194      6
Outlier                 0     10

8.3 Paso 3: Correlación parcial y matriz de precisión

# Inversa de la covarianza = Matriz de precisión
theta_mat = np.linalg.inv(S_mat)

# Cálculo matricial de correlaciones parciales
d_diag = np.sqrt(np.diag(theta_mat))
partial_corr = -theta_mat / np.outer(d_diag, d_diag)
np.fill_diagonal(partial_corr, 1.0)

df_pcorr = pd.DataFrame(partial_corr, index=cols, columns=cols).round(3)
print("Matriz de Correlaciones Parciales (condicionada al resto de variables):")
print(df_pcorr)

Matriz de Correlaciones Parciales (condicionada al resto de variables):
             Temperatura  Presión  Vibración  Consumo
Temperatura        1.000    0.078      0.391    0.191
Presión            0.078    1.000      0.198    0.459
Vibración          0.391    0.198      1.000   -0.042
Consumo            0.191    0.459     -0.042    1.000

8.4 Paso 4: Visualización interactiva con Coordenadas Paralelas

# Visualización en coordenadas paralelas distinguiendo anomalías robustas
df_lab['Diagnostico'] = np.where(df_lab['Outlier_Robusto'], 'Anomalía Multivariada', 'Comportamiento Normal')

fig_par = px.parallel_coordinates(
    df_lab,
    dimensions=cols,
    color='DM2_Robusta',
    color_continuous_scale=px.colors.diverging.Tealrose,
    title="Coordenadas Paralelas del Sistema: Detección de Perfiles Anómalos vía Mahalanobis Robusto"
)
fig_par.update_layout(height=480, template="plotly_white")
fig_par.show()

---

9 Conclusiones y síntesis

El análisis exploratorio de datos multivariado transforma radicalmente la forma en que entendemos un conjunto de datos en minería y ciencia de datos:

1. La localización es vectorial y geométrica: el vector de medias es el centro de masa afín, pero ante valores atípicos la mediana geométrica de Fermat-Weber y la profundidad de Tukey ofrecen estimadores consistentes e invariantes ante rotaciones.
2. La dispersión es dimensional: la varianza total mide la magnitud global de las varianzas pero es ciega a la dependencia, mientras que la varianza generalizada de Wilks captura el hipervolumen real de la nube elipsoidal.
3. La distancia estadística vence a la distancia física: la distancia de Mahalanobis incorpora la covarianza para juzgar qué tan anómalo es un punto, constituyendo la herramienta primaria para detectar valores atípicos multivariados que escapan por completo a los filtros univariados.
4. La dependencia requiere condicionamiento: los coeficientes de correlación bivariados pueden engañar debido a variables confusoras o a la Paradoja de Simpson; la matriz de precisión y la correlación parcial aíslan los enlaces directos de dependencia en el sistema.

Estas bases conceptuales sientan el andamiaje necesario para abordar con éxito las etapas subsiguientes del ciclo CRISP-DM: el preprocesamiento de datos (escalado, imputación multivariada y reducción de dimensionalidad con PCA en el Capítulo de Reducción de Dimensionalidad con PCA) y los algoritmos de aprendizaje no supervisado y clasificación.

---

Referencias

Huber, Peter J., y Elvezio M. Ronchetti. 2009. Robust Statistics. 2.ª ed. Wiley. https://doi.org/10.1002/9780470434697.

Joanes, D. N., y C. A. Gill. 1998. «Comparing Measures of Sample Skewness and Kurtosis». Journal of the Royal Statistical Society: Series D (The Statistician) 47 (1): 183-89. https://doi.org/10.1111/1467-9884.00122.

Johnson, Richard A., y Dean W. Wichern. 2007. Applied Multivariate Statistical Analysis. 6.ª ed. Prentice Hall.

Mahalanobis, Prasanta Chandra. 1936. «On the generalized distance in statistics». Proceedings of the National Institute of Sciences of India 2 (1): 49-55.

Mardia, K. V. 1970. «Measures of multivariate skewness and kurtosis with applications». Biometrika 57 (3): 519-30. https://doi.org/10.1093/biomet/57.3.519.

Rousseeuw, Peter J. 1984. «Least Median of Squares Regression». Journal of the American Statistical Association 79 (388): 871-80. https://doi.org/10.1080/01621459.1984.10477105.

Rousseeuw, Peter J., y Katrien Van Driessen. 1999. «A Fast Algorithm for the Minimum Covariance Determinant Estimator». Technometrics 41 (3): 212-23. https://doi.org/10.1080/00401706.1999.10485670.

Simpson, Edward H. 1951. «The Interpretation of Interaction in Contingency Tables». Journal of the Royal Statistical Society: Series B (Methodological) 13 (2): 238-41. https://doi.org/10.1111/j.2517-6161.1951.tb00088.x.

Small, Christopher G. 1990. «A Survey of Multidimensional Medians». International Statistical Review 58 (3): 263-77. https://doi.org/10.2307/1403809.

Tukey, John W. 1975. «Mathematics and the picturing of data». Proceedings of the International Congress of Mathematicians 2: 523-31.

Wilks, Samuel S. 1932. «Certain generalizations in the analysis of variance». Biometrika 24 (3/4): 471-94. https://doi.org/10.1093/biomet/24.3-4.471.