Qué es PCA, y ya lo usaste

PCA no es un algoritmo nuevo 🌸

Si hiciste el capítulo 10, PCA ya está hecho. Lo único que falta es multiplicar y leer.

Y quiero empezar por lo que PCA no es, porque hay dos ideas sueltas por ahí que hacen daño. No es "quitar columnas": ninguna columna se elimina, todas participan en todas las componentes. Y no es "quedarse con lo importante": lo que se mide es dónde hay variación, y variar mucho no es lo mismo que importar.

Abrimos el archivo con más columnas

Este bloque necesita una tabla más ancha. Con las tres de siempre no hay nada que descomponer, porque están casi sin relación entre ellas. Así que meto dos más: el descuento y el monto final facturado.

import pandas as pd
import numpy as np

COLS = ['unidades', 'monto', 'descuento', 'satisfaccion',
        'monto_final_facturado']

df = pd.read_csv('ventas-miss-yera.csv')
for c in COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')
M = df[COLS].dropna().to_numpy()

Z = (M - M.mean(axis=0)) / M.std(axis=0)
print('M.shape:', M.shape)

M.shape: (1793, 5)

Esa Z es la tabla escalada: a cada columna le resté su promedio y la dividí por su desviación. Después de eso todas tienen media cero y desviación uno, o sea que ninguna manda por tener números más grandes. En el capítulo 7 viste por qué eso hace falta.

Las direcciones, otra vez

C = np.cov(Z, rowvar=False)
vals, vecs = np.linalg.eigh(C)
orden = np.argsort(vals)[::-1]
vals, vecs = vals[orden], vecs[:, orden]
print('autovalores:', np.round(vals, 4))

autovalores: [1.7482 1.0382 0.9982 0.9617 0.2565]

Y PCA es una multiplicación

P = Z @ vecs
print('P.shape:', P.shape)
print('varianza de cada componente:', np.round(P.var(axis=0, ddof=1), 4))

P.shape: (1793, 5)
varianza de cada componente: [1.7482 1.0382 0.9982 0.9617 0.2565]

Ahí está todo el truco. Las varianzas de las componentes son los autovalores, clavados. No se parecen: son los mismos números.

Eso es lo que quería decir "los autovalores reparten la varianza". Cada dirección se lleva un trozo, y el autovalor es el trozo.

Fíjate también en la forma: sigue siendo 1.793 por 5. No se perdió nada. Es la misma tabla girada, no una tabla recortada. El recorte viene después y es una decisión aparte.

Cuánto se lleva cada una

print('proporcion:', np.round(vals / vals.sum(), 4))
print('acumulada :', np.round(np.cumsum(vals / vals.sum()), 4))

proporcion: [0.3495 0.2075 0.1995 0.1922 0.0513]
acumulada : [0.3495 0.557  0.7565 0.9487 1.    ]

Y acá te tengo que dar una noticia mala sobre mis propios datos: este PCA no sirve de mucho.

Mira la primera componente: 34.95%. Y para llegar al 95% hacen falta cuatro de las cinco. O sea que puedo tirar una columna y poco más.

Cuando en un curso te enseñan PCA, el ejemplo siempre tiene la primera componente en 70% u 80%, porque eligieron un dataset donde funciona. En el mío las columnas son bastante independientes entre sí, así que no hay mucho que comprimir.

Eso también es un resultado. "PCA no encontró estructura" quiere decir "tus columnas ya estaban aportando cosas distintas", y eso es una buena noticia para un modelo, aunque sea una mala noticia para el gráfico bonito 🙂

Lo que sí encontró

COLS = ['unidades', 'monto', 'descuento', 'satisfaccion',
        'monto_final_facturado']
for nombre, peso in zip(COLS, vecs[:, 0]):
    print(f'  {nombre:24} {peso:8.4f}')

  unidades                  -0.0000
  monto                      0.6996
  descuento                 -0.0969
  satisfaccion               0.0631
  monto_final_facturado      0.7051

La componente 1 es plata. Monto y monto final facturado con casi el mismo peso, y lo demás en cero.

O sea que si tuvieras que resumir cada venta en un solo número, ese número sería "cuánto factura". Suena obvio dicho así, y el asunto es que nadie se lo dijo: lo encontró la matriz mirando 1.793 filas.

Esa es la utilidad honesta de PCA la mayoría de las veces: no comprimir, sino entender qué se mueve junto.

Recortar, que es lo que la gente quiere

Z2 = P[:, :2] @ vecs[:, :2].T
Z4 = P[:, :4] @ vecs[:, :4].T
print('error medio con 2 componentes:', round(float(np.abs(Z - Z2).mean()), 4))
print('error medio con 4 componentes:', round(float(np.abs(Z - Z4).mean()), 4))

error medio con 2 componentes: 0.4905
error medio con 4 componentes: 0.0966

Con dos componentes reconstruyo la tabla con un error medio de 0.49 sobre datos que tienen desviación 1. O sea que me equivoco casi media desviación en cada casilla, que es muchísimo.

Con cuatro bajo a 0.0966, que ya es razonable. Otra vez lo mismo: acá se puede tirar una columna, no cuatro.

Tres cosas que hay que decir y casi nunca se dicen

Uno: escalar no es opcional. Si no escalas, la columna con los números más grandes se lleva la primera componente entera, y no porque importe más. Es el mismo problema del capítulo 7.

Dos: las componentes no se explican. "Componente 1" no es una variable de negocio. Puedes leer sus pesos y ponerle nombre, como hice arriba con "plata", pero si le entregas a alguien un modelo entrenado sobre componentes, no vas a poder decirle qué significa ninguna.

Tres: PCA no sabe qué quieres predecir. Busca dónde hay variación, no dónde hay respuesta. Puede tirar tranquilamente la componente que menos varía y resulta que era justo la que separaba a tus clientes buenos de los malos. Pasa, y cuando pasa no avisa 🐣

Lo que te llevas

• PCA es multiplicar tus datos escalados por los autovectores.
• La varianza de cada componente es su autovalor.
• PCA gira la tabla, no la recorta. Recortar es una decisión aparte.
• Si la varianza se reparte parejo, no hay nada que comprimir, y eso es un resultado.
• Sin escalar, PCA solo encuentra la columna más grande.
• PCA no sabe qué quieres predecir.

Comprueba que se entendió

Ejercicios

C = np.cov(Z, rowvar=False)
vals, vecs = np.linalg.eigh(C)
P = Z @ vecs
print(np.round(np.corrcoef(P, rowvar=False), 6))

[[ 1.  0.  0.  0. -0.]
 [ 0.  1.  0. -0.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0. -0.  0.  1.  0.]
 [-0.  0.  0.  0.  1.]]

COLS = ['unidades', 'monto', 'descuento', 'satisfaccion',
        'monto_final_facturado']
Cm = np.cov(M - M.mean(axis=0), rowvar=False)
vm, wm = np.linalg.eigh(Cm)
o = np.argsort(vm)[::-1]
vm, wm = vm[o], wm[:, o]
print('proporcion sin escalar:', np.round(vm / vm.sum(), 4))
for nombre, peso in zip(COLS, wm[:, 0]):
    print(f'  {nombre:24} {peso:8.4f}')

proporcion sin escalar: [0.8715 0.1285 0.     0.     0.    ]
  unidades                  -0.0000
  monto                     -0.7696
  descuento                  0.0000
  satisfaccion              -0.0000
  monto_final_facturado     -0.6385

C = np.cov(Z, rowvar=False)
vals, vecs = np.linalg.eigh(C)
orden = np.argsort(vals)[::-1]
vals, vecs = vals[orden], vecs[:, orden]
P = Z @ vecs

for k in range(1, 6):
    Zk = P[:, :k] @ vecs[:, :k].T
    print(k, 'componentes -> error', round(float(np.abs(Z - Zk).mean()), 4))

1 componentes -> error 0.5975
2 componentes -> error 0.4905
3 componentes -> error 0.3578
4 componentes -> error 0.0966
5 componentes -> error 0.0