Principal Component Analysis

Un poco de historia

Análisis de Componentes Principales
Análisis de Componentes Principales

Los componentes principales son combinaciones lineales de variables aleatorias que poseen propiedades especiales en términos de varianzas. Transformar la variable vectorial original en el vector de componentes principales equivale a una rotación de los ejes de coordenadas hacia un nuevo sistema de coordenadas que posee propiedades estadísticas inherentes. Esta elección del sistema de coordenadas contrasta con muchos de los problemas tratados anteriormente, en los cuales el sistema de coordenadas era irrelevante.

Los componentes principales resultan ser los vectores característicos (autovectores) de la matriz de covarianza. Así, el estudio de los componentes principales puede considerarse como una formulación estadística de los desarrollos habituales de raíces y vectores característicos (para matrices semidefinidas positivas) (Anderson, 2003).

Un problema central en el análisis de componentes principales es la reducción de dimensionalidad. Hablamos de reducción de dimensionalidad si es posible describir \(p\) variables por \(r\) variables tal que \(r < p\), ello a costa de una pequeña pérdida de información.

El objetivo de componentes principales es precisamente ello, dadas n observaciones de p variables (atributos), se analiza si es posible representar dicha información con un número menor de variables.

La técnica de componentes principales es debida a Hotelling(\(1933\)), pero los origenes pueden encontarse en ajuste ortogonales por mínimos cuadrados por Karl Pearson.

En general, la técnica de PCA (principal components analysis) tiene dos principales utilidades:

  • Permite representar observaciones de un espacio p-dimensional en uno de menor dimensión, permitiendo también, saber que variables afectan más la variabilidad de los datos.

  • Cuando existen variables altamente correlacionadas, el PCA permite transformar dichas variables en nuevas variables no correlacionadas.

Teoría de PCA

Supongamos que el vector aleatorio \(𝑋\), de \(𝑝\) componentes

\[\textbf{X}= \left( \begin{array}{cccc} X_{1,1} & X_{1,2} & ... & X_{1,p}\\ X_{2,1} & X_{2,2} & ... & X_{2,p}\\ \vdots & \vdots & \vdots & \vdots\\ X_{n,1} & X_{n,2} & ... & X_{n,p}\\ \end{array}\right)\]

tiene como matriz de covarianza a \(\Sigma\), donde cada renglón reprensenta a cada individuo, es decir, \(X_{i,1},...,X_{i,p}\)

\[Var(\textbf{X})=\Sigma = \left( \begin{array}{cccc} \sigma_1^2 & \sigma_{12} & ... & \sigma_{1p}\\ \sigma_{21}^2 & \sigma_{2}^2 & ... & \sigma_{2p}\\ \vdots & \vdots & ... & \vdots \\ \sigma_{n1} & \sigma_{2n} & ... & \sigma_{n}^2 \end{array}\right)\]

La distribución real de \(X\) es irrelevante, excepto por su matriz de covarianza; sin embargo, si \(X\) está distribuido normalmente, puede otorgarse un significado más profundo a los componentes principales.

El desarrollo incluirá los casos en los que \(\Sigma\) sea singular (es decir, semidefinida positiva) y aquellos en los que \(\Sigma\) tenga raíces múltiples.

Sea \(\mathbf{e}\) un vector columna de \(p\) componentes tal que \(\mathbf{e}'\mathbf{e} = 1\). La varianza de \(\mathbf{Y} = \mathbf{e}'\mathbf{X}\) es:

\[ \text{Var}(Y) = \text{Var}(\mathbf{e}'\mathbf{X}) \]

en forma explícita

\[\begin{array}{cccccc} Y_1 & = & e_{11} X_1 + & e_{12} X_2 + & ... + & e_{1p} X_p\\ Y_2 & = & e_{21} X_1 + & e_{22} X_2 + & ... + & e_{2p} X_p\\ & \vdots & & & & \\ Y_{n} & = & e_{p1} X_1 + & e_{p2} X_2 + & ... + & e_{pp} X_p \end{array}\]

Con varianza poblacional \[\begin{equation} Var(Y) = \textbf{e}^{'} \Sigma \textbf{e} \end{equation}\]

Para determinar la combinación lineal normalizada \(\mathbf{e}'\mathbf{X}\) con varianza máxima, debemos encontrar un vector \(\mathbf{e}\) que satisfaga \(\mathbf{e}'\mathbf{e} = 1\) y que maximice

\[ \phi = \mathbf{e}' \boldsymbol{\Sigma} \mathbf{e} - \lambda(\mathbf{e}'\mathbf{e} - 1) = \sum_i \sum_j e_i \sigma_{ij} e_j - \lambda\left(\sum_i e_i^2 - 1\right), \]

donde \(\lambda\) es un multiplicador de Lagrange.

El vector de derivadas parciales \((\partial \phi / \partial e_j)\) es

\[\begin{eqnarray*} \frac{\partial \phi} {\partial e_j} & = & 2 \Sigma \textbf{e} -2\lambda \textbf{e}\\ 0 & = & (\Sigma - \lambda )\textbf{e} \end{eqnarray*}\]

\[ (\boldsymbol{\Sigma} - \lambda \mathbf{I})\mathbf{e} = \mathbf{0}. \]

como \(\mathbf{e}'\boldsymbol{\Sigma}\mathbf{e}\) y \(\mathbf{e}'\mathbf{e}\) tienen derivadas en toda una región que contiene \(\mathbf{e}'\mathbf{e} = 1\), el vector \(\mathbf{e}\) que maximiza \(\mathbf{e}'\boldsymbol{\Sigma}\mathbf{e}\) debe satisfacer la ecuación anterior.

Para obtener una solución de \((\boldsymbol{\Sigma} - \lambda \mathbf{I})\mathbf{e} = \mathbf{0}\) con \(\mathbf{e}'\mathbf{e} = 1\), la matriz \(\boldsymbol{\Sigma} - \lambda \mathbf{I}\) debe ser singular; en otras palabras, \(\lambda\) debe satisfacer

\[ \det(\boldsymbol{\Sigma} - \lambda \mathbf{I}) = 0. \]

La función \(\det(\boldsymbol{\Sigma} - \lambda \mathbf{I})\) es un polinomio en \(\lambda\) de grado \(p\). Por lo tanto, la ecuación anterior tiene \(p\) raíces; sean éstas \[ \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p. \]

Por lo que

\[ \mathbf{e}'\boldsymbol{\Sigma}\mathbf{e} = \lambda \mathbf{e}'\mathbf{e}. \]

entonces, la varianza de \(\mathbf{e}'\mathbf{X}\) es \(\lambda\). Por lo tanto, para obtener la varianza máxima, basta con obtener la raíz mayor \(\lambda_1\).

Sea \(\mathbf{e}^{(1)}\) una solución normalizada de \[ (\boldsymbol{\Sigma} - \lambda_1 \mathbf{I})\mathbf{e} = \mathbf{0}. \] entonces \[ U_1 = \mathbf{e}^{(1)'}\mathbf{X} \] es una combinación lineal normalizada con varianza máxima.

Ahora, busquemos una combinación normalizada \(\mathbf{e}'\mathbf{X}\) que tenga la máxima varianza entre todas las combinaciones lineales no correlacionadas con \(U_1\). La falta de correlación implica que

\[ E(U_1 U) = 0, \]

ya que \(\boldsymbol{\Sigma}\mathbf{e}^{(1)} = \lambda_1 \mathbf{e}^{(1)}\). Por lo tanto, \(\mathbf{e}'\mathbf{X}\) es ortogonal a \(U_1\) tanto en el sentido estadístico (falta de correlación) como en el sentido geométrico (producto interno nulo entre los vectores \(\mathbf{e}\) y \(\mathbf{e}^{(1)}\)).

Fuente:

  • T. W. Anderson (2003). An Introduction to Multivariate Statistical Analysis. Third Ediion. Wiley-Interscience

Ejemplo en python:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from numpy.linalg import eig
from sklearn.decomposition import PCA
import random
sns.set_theme()

Consideramos las variables \(x_1\) y \(x_2\)


x1 = np.linspace(0,20, 50)
x2 = x1 + np.random.normal(0,3, len(x1))
# dataframe
X = pd.DataFrame({'x1':x1, 'x2':x2})
plt.figure(figsize=(16,8))
plt.plot(x1, x2, 'o')
plt.xlabel('x1')
plt.ylabel('x2')
plt.title('Variables originales')
plt.show()

Los siguientes gráficos muestran la aplicación del Análisis de Componentes Principales (PCA) para reducir la dimensionalidad de un conjunto de datos con dos variables reales, \(x_1\), y \(x_2\). Primero se grafican las variables originales (extrema derecha), mostrando cómo se distribuyen los puntos en el plano y revelando su posible correlación. Luego, tras aplicar el PCA, se representan las dos componentes principales (CP1 y CP2) (gráfica intermedia), que son nuevas variables no correlacionadas que capturan la máxima varianza del conjunto: CP1 explica la mayor parte de la variabilidad y CP2 la que queda de forma ortogonal. Finalmente, se muestra la proyección sobre la primera componente principal (gráfica extrema izquierda), reduciendo los datos a una sola dimensión donde se conserva la mayor parte de la información posible. En conjunto, las tres gráficas ilustran cómo el PCA transforma y simplifica los datos manteniendo su estructura esencial.

pca = PCA()
X_red = pca.fit_transform(X)

fig, axs = plt.subplots(1,3,figsize=(18,5))
axs[0].scatter(X['x1'],X['x2'],color='blue')
axs[0].set_xlabel('x1')
axs[0].set_ylabel('x2')
axs[0].set_title('Variables originales')
axs[1].scatter(X_red[:,0],X_red[:,1],alpha=0.7,color='blue',label="PCA")
axs[1].set_xlabel('CP1')
axs[1].set_ylabel('CP2')
axs[1].set_title('Componentes principales')
axs[1].plot([-15,20],[0,0],color='gray')
axs[1].plot([0,0],[-5,5],color='gray')
axs[1].legend(loc='best')
axs[2].scatter(X_red[:,0],[0 for x in X_red],color='blue',label="PCA (1st PC)")
axs[2].plot([-15,20],[0,0],color='gray')
axs[2].set_xlabel('CP1')
axs[2].set_title('1ra Componente Principal')
axs[2].legend(loc='best')
fig.tight_layout()
fig.show()

Caso: Notas de estudiantes

Supongamos que se tiene los datos de \(15\) estudiantes con sus respectivas calificaciones en cada asignatura(\(8\) asignaturas). Supongamos que se desea analizar todas las variables en una escala de dimensión menor, para mayor referencia Daniel Peña (2002).

Alumno Lengua Matemáticas Física Inglés Filosofía Historia Química Educ. Física
1 5 5 5 5 5 5 5 5
2 7 4 3 8 4 7 3 8
3 5 8 7 6 5 6 7 5
4 7 2 4 8 7 7 3 6
5 8 9 10 8 8 7 9 4
6 4 9 8 4 3 4 7 5
7 6 4 4 6 5 5 3 7
8 4 7 8 3 3 2 8 3
9 5 5 4 5 6 5 5 1
10 7 4 5 7 8 8 4 6
11 7 8 8 7 7 6 7 9
12 4 3 3 4 3 2 1 4
13 7 4 4 7 8 7 4 5
14 3 5 5 2 3 3 5 7
15 5 6 6 5 5 5 6 6 
lengua = [5, 7, 5, 7, 8, 4, 6, 4, 5, 7, 7, 4, 7, 3, 5]
matematicas = [5, 4, 8, 2, 9, 9, 4, 7, 5, 4, 8, 3, 4, 5, 6]
fisica = [5, 3, 7, 4, 10, 8, 4, 8, 4, 5, 8, 3, 4, 5, 6]
ingles = [5, 8, 6, 8, 8, 4, 6, 3, 5, 7, 7, 4, 7, 2, 5]
filosofia = [5, 4, 5, 7, 8, 3, 5, 3, 6, 8, 7, 3, 8, 3, 5]
historia = [5, 7, 6, 7, 7, 4, 5, 2, 5, 8, 6, 2, 7, 3, 5]
quimica = [5, 3, 7, 3, 9, 7, 3, 8, 5, 4, 7, 1, 4, 5, 6]
efisica = [5, 8, 5, 6, 4, 5, 7, 3, 1, 6, 9, 4, 5, 7, 6]

estudiantes = pd.DataFrame({'lengua':lengua,
                        'matemáticas':matematicas,
                        'física':fisica,
                        'inglés':ingles,
                        'filosofía':filosofia,
                        'historia':historia,
                        'química':quimica,
                        'efísica':efisica})
estudiantes
##     lengua  matemáticas  física  inglés  filosofía  historia  química  efísica
## 0        5            5       5       5          5         5        5        5
## 1        7            4       3       8          4         7        3        8
## 2        5            8       7       6          5         6        7        5
## 3        7            2       4       8          7         7        3        6
## 4        8            9      10       8          8         7        9        4
## 5        4            9       8       4          3         4        7        5
## 6        6            4       4       6          5         5        3        7
## 7        4            7       8       3          3         2        8        3
## 8        5            5       4       5          6         5        5        1
## 9        7            4       5       7          8         8        4        6
## 10       7            8       8       7          7         6        7        9
## 11       4            3       3       4          3         2        1        4
## 12       7            4       4       7          8         7        4        5
## 13       3            5       5       2          3         3        5        7
## 14       5            6       6       5          5         5        6        6

Como primer paso centramos los datos \[ \bar{x}_j = \frac{1}{n} \sum_{i=1}^{n} x_{ij}, \quad j = 1, 2, \dots, p \]

\[ x'_{ij} = x_{ij} - \bar{x}_j \]

\[ X_c = X - \mathbf{1}_n \bar{x}^\top \]

estudiantes_center = estudiantes - estudiantes.mean(axis=0)
estudiantes_center
##     lengua  matemáticas  física  ...  historia   química  efísica
## 0     -0.6    -0.533333    -0.6  ... -0.266667 -0.133333     -0.4
## 1      1.4    -1.533333    -2.6  ...  1.733333 -2.133333      2.6
## 2     -0.6     2.466667     1.4  ...  0.733333  1.866667     -0.4
## 3      1.4    -3.533333    -1.6  ...  1.733333 -2.133333      0.6
## 4      2.4     3.466667     4.4  ...  1.733333  3.866667     -1.4
## 5     -1.6     3.466667     2.4  ... -1.266667  1.866667     -0.4
## 6      0.4    -1.533333    -1.6  ... -0.266667 -2.133333      1.6
## 7     -1.6     1.466667     2.4  ... -3.266667  2.866667     -2.4
## 8     -0.6    -0.533333    -1.6  ... -0.266667 -0.133333     -4.4
## 9      1.4    -1.533333    -0.6  ...  2.733333 -1.133333      0.6
## 10     1.4     2.466667     2.4  ...  0.733333  1.866667      3.6
## 11    -1.6    -2.533333    -2.6  ... -3.266667 -4.133333     -1.4
## 12     1.4    -1.533333    -1.6  ...  1.733333 -1.133333     -0.4
## 13    -2.6    -0.533333    -0.6  ... -2.266667 -0.133333      1.6
## 14    -0.6     0.466667     0.4  ... -0.266667  0.866667      0.6
## 
## [15 rows x 8 columns]

Ahora calculamos la matriz de covarianza \[ S = \frac{1}{n-1} X_c^{\top} X_c \]

en forma expandida:

\[ s_{jk} = \frac{1}{n-1} \sum_{i=1}^{n} (x_{ij} - \bar{x}_j)(x_{ik} - \bar{x}_k), \quad j,k = 1,2,\dots,p \]

mat_cov = np.cov(estudiantes_center.T)
print(pd.DataFrame(mat_cov))
##           0         1         2  ...         5         6         7
## 0  2.257143 -0.271429  0.185714  ...  2.471429 -0.014286  0.885714
## 1 -0.271429  4.838095  4.300000  ... -0.366667  4.423810 -0.371429
## 2  0.185714  4.300000  4.542857  ... -0.100000  4.414286 -0.328571
## 3  2.714286 -0.523810 -0.142857  ...  3.166667 -0.380952  1.142857
## 4  2.428571 -0.261905  0.428571  ...  2.976190  0.380952  0.214286
## 5  2.471429 -0.366667 -0.100000  ...  3.495238 -0.038095  1.171429
## 6 -0.014286  4.423810  4.414286  ... -0.038095  4.838095 -0.771429
## 7  0.885714 -0.371429 -0.328571  ...  1.171429 -0.771429  3.971429
## 
## [8 rows x 8 columns]

Calculo de valores propios y vectores propios

Calculamos valores propios y vectores propios

# eigenvalues
values, vectors = eig(mat_cov)
print(f"Valores propios{values}")
## Valores propios[13.74809013 11.82720931  3.77290536  0.79666945  0.53112604  0.04103244
##   0.12790717  0.28839343]
print(f"Vectores propios{vectors}")
## Vectores propios[[-0.10537175  0.40711401  0.03314967 -0.15405155  0.39215124  0.71550906
##    0.35074162 -0.09781973]
##  [ 0.56246248  0.11900302 -0.13283045 -0.34720521 -0.22998894  0.06448064
##    0.20026394  0.65887825]
##  [ 0.53363953  0.18912541 -0.06515149  0.13366458  0.47189382  0.05397228
##   -0.65283708 -0.07454162]
##  [-0.16878408  0.4900049   0.01967903 -0.54659238  0.2593807  -0.59591898
##    0.03877109 -0.09140934]
##  [-0.08821408  0.48702046  0.30534414  0.66549573  0.04666778 -0.20503592
##    0.14686702  0.39130832]
##  [-0.14562687  0.49715281  0.02560583 -0.0972865  -0.66876246  0.24043353
##   -0.44582481 -0.13304252]
##  [ 0.55931708  0.17184302  0.04106607  0.16492072 -0.23217433 -0.14685022
##    0.42628131 -0.60956005]
##  [-0.13209482  0.17419349 -0.93864161  0.24400642  0.00585955 -0.06666033
##    0.08443805  0.00355905]]

Ordenamos de forma descendente los valores propios, y con base en ello también lo hacemos para los vectores propios.

# ordena valores propios de forma descendente y retorna el índice
sorted_index = np.argsort(values)[::-1]
# ordenamos los valores propios con respecto al índice
sorted_values = values[sorted_index]
# ordenamos los vectores propios con respecto al índice
sorted_vectors = vectors[:, sorted_index]

veamamos los vectores propios

# eigenvectores
pd.DataFrame(vectors)
##           0         1         2  ...         5         6         7
## 0 -0.105372  0.407114  0.033150  ...  0.715509  0.350742 -0.097820
## 1  0.562462  0.119003 -0.132830  ...  0.064481  0.200264  0.658878
## 2  0.533640  0.189125 -0.065151  ...  0.053972 -0.652837 -0.074542
## 3 -0.168784  0.490005  0.019679  ... -0.595919  0.038771 -0.091409
## 4 -0.088214  0.487020  0.305344  ... -0.205036  0.146867  0.391308
## 5 -0.145627  0.497153  0.025606  ...  0.240434 -0.445825 -0.133043
## 6  0.559317  0.171843  0.041066  ... -0.146850  0.426281 -0.609560
## 7 -0.132095  0.174193 -0.938642  ... -0.066660  0.084438  0.003559
## 
## [8 rows x 8 columns]

Número de componentes a retener

  • Retener suficientes componentes para garantizar un porcentaje predefinido de la varianza total, por ejemplo el \(80\%\).
  • Retener aquellas componentes cuyos autovalores superen el promedio de todos los autovalores, \(\sum_{i=1}^{p} \lambda_i/p\). Para la matriz de correlación este promedio es \(1.0\).
  • Utilizar una representación ggráfica de \(\lambda_i\) frente a \(i\), y determinar el punto en el que se produce la transición entre los autovalores grandes y los pequeños.
  • Utilizar tests de significancia.

Para probar \(H_{0k}\) se calcula el promedio de los últimos \(k\) autovalores

\[\begin{eqnarray*} \bar{\lambda} = \sum_{i=p-k+1}^{p} \frac{\lambda_i}{k} \end{eqnarray*} \] y se calcula el estadístico \[ \begin{eqnarray*} u = \left(n - \frac{2p+11}{6} \right)\left(k \ ln(\bar{\lambda}) - \sum_{i=p-k+1}^{p} ln(\lambda_i) \right) \end{eqnarray*} \] que sigue aproximadamente una distribución \(\chi_v^2\). Por lo que se rechaza \(H_0\) si \(u \leq \chi_v^2\)

Para el ejemplo de las calificaciones, considerando los siguientes datos, sobre la proporción de la varianza explicada por las componentes, se considerá \(2\) componentes como apropiado, ya que representan el \(82\%\) de la variabilidad. Para mayor descripción sobre el método puede consultarse Rencher (2003).

# desviación estándar
sd = np.sqrt(sorted_values)
# proporción de la varianza explicada
prop_var = sorted_values/np.sum(sorted_values)
# proporción de la varianza explicada acumulada
cum_prop = np.cumsum(prop_var)
# nombramos columnas
cp = ['CP1', 'CP2', 'CP3', 'CP4', 'CP5', 'CP6', 'CP7', 'CP8']
# data frame
tab_summary = pd.DataFrame({'CP':cp, 'sd':sd, 'prop_var':prop_var, 'cum_prop':cum_prop})
tab_summary = tab_summary.set_index('CP')
tab_summary
##            sd  prop_var  cum_prop
## CP                               
## CP1  3.707842  0.441587  0.441587
## CP2  3.439071  0.379889  0.821476
## CP3  1.942397  0.121185  0.942662
## CP4  0.892563  0.025589  0.968251
## CP5  0.728784  0.017060  0.985311
## CP6  0.537023  0.009263  0.994574
## CP7  0.357641  0.004108  0.998682
## CP8  0.202565  0.001318  1.000000

Graficamos la proporción de la varianza explicada

plt.figure(figsize=(16,8))
plt.bar(tab_summary.index, tab_summary['prop_var'])
plt.plot(tab_summary.index, tab_summary['cum_prop'])
plt.plot(tab_summary.index, tab_summary['cum_prop'], 'o')
plt.title('Proporción de la varianza explicada')
plt.xlabel('Componentes principales')
# etiquetas de la proporción
for i in range(len(prop_var)):
  plt.text(i,prop_var[i], str(np.round(prop_var[i],2)*100) + '%')
# etiquetas de la proporción acumulada
for i in range(len(cum_prop)):
  plt.text(i, cum_prop[i], str(np.round(cum_prop[i],2)*100) + '%')

Las componentes resultantes son: \[ \begin{eqnarray*} Y_1 & = & -0.1053 X_1 + 0.5624 X_2 + 0.5336 X_3 -0.1687 X_4 -0.0882X_5 -0.1456X_6 + 0.5593X_7 - 0.1320X_8\\ \\ Y_2 & = & 0.4071 X_1 + 0.1190 X_2 + 0.1891 X_3 + 0.4900 X_4 + 0.4870 X_5 + 0.4971 X_6 + 0.1718X_7 + 0.1781X_8\\ \\ & ... &\\ \\ Y_8 & = & 0.7155 X_1 + 0.0644 X_2 + 0.0539 X_3 + 0.5959 X_4 - 0.2050 X_5 + 0.2404 X_6 - 0.1468 X_7 -0.0666 X_8\\ \\ \end{eqnarray*} \] donde: \(X_1\) = Lengua, \(X_2\) = Matemáticas, \(X_3\) = Física, \(X_4\) = Inglés, \(X_5\) = Filosofía, \(X_6\) = Historia, \(X_7\) = Química, \(X_8\) = Educación física\

Componentes a retener

# número de componentes a retener
k = 8
# filtramos
vectors_subset = sorted_vectors[:,0:k]
cp_subset = cp[0:k]

Realizamos el producto punto del vector anterior con el de los datos centrados de las notas de estudiantes

Y_components = np.dot(vectors_subset.T , estudiantes_center.T ).T
# dataframe
CPs = pd.DataFrame(Y_components, columns=cp_subset)
CPs
##          CP1       CP2       CP3  ...       CP6       CP7       CP8
## 0  -0.397917 -1.135386  0.338296  ... -0.202152  0.027921 -0.048325
## 1  -4.462713  1.337783 -2.425423  ... -0.609411  0.313362  0.202133
## 2   3.160973  0.930718 -0.063081  ...  0.181827 -0.231453 -0.396125
## 3  -5.054451  2.401576  0.568402  ... -0.834902 -0.468278 -0.354643
## 4   5.511130  5.946156  1.709853  ... -0.040994 -0.097032  0.183248
## 5   5.167696 -2.116624 -0.995471  ...  0.530271 -0.514394  0.127862
## 6  -3.150999 -0.541694 -1.370309  ...  0.250520  0.186320  0.113192
## 7   5.326316 -4.015485  1.117649  ... -1.046670  0.195362  0.100423
## 8  -0.491391 -1.534265  4.463358  ...  0.249461  0.489873 -0.040692
## 9  -2.901626  3.494719  0.589926  ...  0.148428 -0.632034  0.312757
## 10  2.610277  4.094891 -3.186132  ...  0.617090  0.487451 -0.183766
## 11 -3.807830 -5.975827  0.768303  ...  0.869595 -0.202269 -0.061989
## 12 -3.157544  2.634247  1.568114  ...  0.352453  0.382189  0.085011
## 13  0.522671 -5.039588 -2.326223  ... -0.241698 -0.023084  0.104298
## 14  1.125407 -0.481221 -0.757261  ... -0.223817  0.086067 -0.143383
## 
## [15 rows x 8 columns]
plt.figure(figsize=(16,8))
plt.plot(CPs['CP1'], CPs['CP2'], 'o')
plt.xlabel('CP1')
plt.ylabel('CP2')
plt.title('Componentes principales')
for i in range(len(CPs['CP1'])):
  plt.text( CPs.iloc[i,0], CPs.iloc[i,1], CPs.index[i])
plt.show()

Módulo PCA() en python

from sklearn.decomposition import PCA
import random

PCA(n_components=None, svd_solver=‘auto’, n_oversamples=10, …)

  • n_components: número de componentes a retener.
  • svd_solver:
    • arpack : SVD truncado a n_components.
    • full : SVD completo.
    • auto : El solver es seleccionado por defecto.
from sklearn.decomposition import PCA
import random
pca = PCA()
Y = pca.fit_transform(X)
pd.DataFrame(Y)
##             0         1
## 0  -10.074595 -4.227641
## 1  -12.751586 -1.126515
## 2  -12.942233 -0.383212
## 3  -13.931173  1.117118
## 4  -11.649117 -0.484469
## 5  -13.180611  1.530371
## 6  -11.839432  0.821027
## 7   -6.564074 -3.619137
## 8   -7.157859 -2.493536
## 9   -1.154895 -7.623696
## 10 -12.111924  3.329471
## 11  -5.642393 -2.243140
## 12  -8.394051  0.928795
## 13  -5.589450 -1.168327
## 14  -7.366808  1.079669
## 15  -9.486265  3.652080
## 16  -1.949057 -2.933019
## 17  -2.365451 -1.975640
## 18  -2.017432 -1.743160
## 19  -3.585385  0.306254
## 20  -3.581099  0.864700
## 21  -0.813865 -1.196986
## 22   0.184911 -1.581626
## 23  -1.366371  0.451980
## 24   0.944648 -1.177074
## 25  -0.553235  0.805893
## 26  -2.555690  3.267350
## 27   5.316892 -3.635788
## 28   0.371833  1.616170
## 29   1.109063  1.479557
## 30   4.554457 -1.225235
## 31  -0.453496  4.086366
## 32   5.017911 -0.539714
## 33   3.464951  1.495482
## 34   2.806212  2.682680
## 35   4.298242  1.830283
## 36   8.386621 -1.484258
## 37   6.992924  0.399908
## 38   6.955247  0.998147
## 39   4.444644  3.941486
## 40   6.814263  2.256862
## 41  10.990578 -1.141070
## 42   6.194983  3.969150
## 43  12.878387 -1.806279
## 44   8.349280  3.051228
## 45  11.745463  0.393103
## 46  12.487291  0.252130
## 47  13.143232  0.192603
## 48  14.718367 -0.738603
## 49  16.907147 -2.251736

Los valores singulares

print(pca.singular_values_)
## [57.93653593 16.61593418]

La proporción acumulada de la varianza explicada

print(pca.explained_variance_ratio_.cumsum())
## [0.92399954 1.        ]

Correlación entre las componentes principales y variables originales

# número de componentes a retener
k = 2
# subconjunto
vectors_subset = sorted_vectors[:,0:k]
cp_subset = cp[0:k]
# componentes
Y_components = np.dot(vectors_subset.T , estudiantes_center.T ).T
CPs = pd.DataFrame(Y_components, columns=cp_subset)
CPs
##          CP1       CP2
## 0  -0.397917 -1.135386
## 1  -4.462713  1.337783
## 2   3.160973  0.930718
## 3  -5.054451  2.401576
## 4   5.511130  5.946156
## 5   5.167696 -2.116624
## 6  -3.150999 -0.541694
## 7   5.326316 -4.015485
## 8  -0.491391 -1.534265
## 9  -2.901626  3.494719
## 10  2.610277  4.094891
## 11 -3.807830 -5.975827
## 12 -3.157544  2.634247
## 13  0.522671 -5.039588
## 14  1.125407 -0.481221
cor1 = []
cor2 = []
for i in range(len(estudiantes_center.columns)):
  cor1.append(CPs['CP1'].corr(estudiantes.iloc[:,i]))
  cor2.append(CPs['CP2'].corr(estudiantes.iloc[:,i]))

df_cor = pd.DataFrame({'cor1':cor1, 'cor2':cor2})

Finalmente graficamos

plt.figure(figsize=(16,8))
plt.xlabel('CP1')
plt.ylabel('CP2')
plt.title('Correlación entre componentes principales y variables originales')
plt.plot(df_cor['cor1'], df_cor['cor2'], 'o', color='orange', markersize=30)
for i in range(len(df_cor['cor1'])):
  plt.text( df_cor.iloc[i,0], df_cor.iloc[i,1], estudiantes_center.columns[i])
plt.show()

Referencias:

  • T. W. Anderson (2003). An Introduction to Multivariate Statistical Analysis. Third Ediion. Wiley-Interscience
  • Alvin C. Rencher (2003). Methods of multivariate analysis. Jhon Wiley. Wiley series in probability and mathematical statistics.
  • Peña Daniel (2002). Análisis de datos multivariantes. Mc Graw Hill Interamericana.