Reducir dimensiones: PCA, t-SNE y UMAP

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n = 200
# Datos sintéticos: notas de Matemática y de Física de 200 estudiantes
mate = rng.normal(60, 12, size=n)
fisica = 60 + 0.8 * (mate - 60) + rng.normal(0, 6, size=n)
X = np.column_stack([mate, fisica])       # forma (200, 2)

# Centrar: restar la media de cada columna
Xc = X - X.mean(axis=0)
print("Forma:", X.shape)
print("Medias antes:", X.mean(axis=0).round(1))
print("Medias después:", Xc.mean(axis=0).round(1))

def sombra(Xc, angulo):
    """Varianza de la sombra y error de reconstrucción sobre una recta."""
    rad = np.deg2rad(angulo)
    u = np.array([np.cos(rad), np.sin(rad)])  # dirección, de largo 1
    t = Xc @ u                    # la coordenada de cada sombra
    reconstruido = np.outer(t, u)  # la sombra, de vuelta en el plano
    varianza = np.mean(t ** 2)
    error = np.mean(np.sum((Xc - reconstruido) ** 2, axis=1))
    return varianza, error

total = np.mean(np.sum(Xc ** 2, axis=1))   # varianza total de la nube
print(f"Varianza total: {total:.1f}")
for angulo in [0, 45, 90]:
    v, e = sombra(Xc, angulo)
    print(f"{angulo:3d}°:  varianza {v:6.1f}  + error {e:6.1f}"
          f"  = {v + e:.1f}")

angulos = np.arange(180)
resultados = np.array([sombra(Xc, a) for a in angulos])
varianzas, errores = resultados[:, 0], resultados[:, 1]
mejor = angulos[np.argmax(varianzas)]
print("Ángulo con más varianza:", mejor, "grados")
print("Ángulo con menos error: ", angulos[np.argmin(errores)], "grados")

plt.plot(angulos, varianzas, color="C0", label="Varianza")
plt.plot(angulos, errores, color="C2", linestyle="--",
         label="Error")
plt.scatter([mejor], [varianzas.max()], color="C1", s=120, zorder=3)
plt.ylim(0, 320)
plt.xlabel("Ángulo de la recta (grados)")
plt.ylabel("Puntos al cuadrado")
plt.title("Girar la recta reparte la varianza")
plt.legend(loc="upper center", ncol=2)
plt.show()

# Matriz de covarianza: (2, 200) @ (200, 2) da (2, 2)
C = Xc.T @ Xc / n
print("Covarianza:\n", C.round(1))

valores, vectores = np.linalg.eigh(C)   # de menor a mayor
print("Valores propios:", valores.round(1))
print("Vectores propios (columnas):\n", vectores.round(3))

u1 = vectores[:, 1]             # el del valor propio más grande
# El ángulo de la flecha u1, llevado al rango de 0° a 180°
angulo_u1 = np.rad2deg(np.arctan2(u1[1], u1[0])) % 180
print(f"Ángulo de la primera componente: {angulo_u1:.1f}°")

# 1. Un elemento de C a mano: el promedio de los productos
print(f"Promedio de mate·física: {np.mean(Xc[:, 0] * Xc[:, 1]):.1f}")

# 2. u·C·u frente a la varianza de la sombra, con la recta de 45°.
# Con un vector a la izquierda, @ hace lo mismo que la transpuesta uᵀ
u45 = np.array([np.cos(np.deg2rad(45)), np.sin(np.deg2rad(45))])
print(f"u·C·u a 45°: {u45 @ C @ u45:.1f}"
      f"   sombra a 45°: {sombra(Xc, 45)[0]:.1f}")

# 3. C empuja las direcciones hacia la nube; al vector propio, no
for angulo in [0, 90]:
    u = np.array([np.cos(np.deg2rad(angulo)), np.sin(np.deg2rad(angulo))])
    Cu = C @ u
    girado = np.rad2deg(np.arctan2(Cu[1], Cu[0]))
    print(f"u a {angulo:2d}°:  C·u apunta a {girado:.1f}°")
print("C·u1: ", (C @ u1).round(1))
print("λ·u1: ", (valores[1] * u1).round(1))
print(f"Máximo de los 180 ángulos: {varianzas.max():.1f}"
      f"   valor propio: {valores[1]:.1f}")

centro = X.mean(axis=0)
plt.scatter(X[:, 0], X[:, 1], s=14, color="C2", alpha=0.5)
for k, color in [(1, "C1"), (0, "C0")]:
    # Flecha de largo 2 desviaciones estándar en cada dirección;
    # el signo de un vector propio es arbitrario: la giramos hacia +x
    flecha = 2 * np.sqrt(valores[k]) * vectores[:, k]
    flecha = flecha * np.sign(flecha[0])
    plt.annotate("", xy=centro + flecha, xytext=centro,
                 arrowprops=dict(arrowstyle="-|>", color=color, lw=3,
                                 mutation_scale=22))
plt.axis("equal")
plt.xlabel("Matemática")
plt.ylabel("Física")
plt.title("Las dos componentes principales")
plt.show()

from sklearn.datasets import load_digits

digitos = load_digits()
Xd, yd = digitos.data, digitos.target
print("Forma de los dígitos:", Xd.shape)

def pca(X, k):
    """PCA desde cero: centrar, covarianza y vectores propios."""
    media = X.mean(axis=0)
    Xc = X - media
    C = Xc.T @ Xc / len(X)
    valores, vectores = np.linalg.eigh(C)
    orden = np.argsort(valores)[::-1]     # de mayor a menor
    valores, vectores = valores[orden], vectores[:, orden]
    componentes = vectores[:, :k].T        # forma (k, columnas)
    Z = Xc @ componentes.T                 # coordenadas nuevas
    return Z, componentes, valores / valores.sum()

Z, comps, proporcion = pca(Xd, 2)
print("Forma de Z:", Z.shape)
print("Proporción de varianza, 5 primeras:", proporcion[:5].round(3))

from sklearn.decomposition import PCA

modelo = PCA(n_components=2).fit(Xd)
Z_sk = modelo.transform(Xd)
print("scikit-learn:", modelo.explained_variance_ratio_.round(3))
# Cada componente puede salir con el signo contrario: da igual
print("Mismas direcciones:",
      np.allclose(np.abs(modelo.components_), np.abs(comps)))
print("Mismas coordenadas:", np.allclose(np.abs(Z_sk), np.abs(Z)))

completo = PCA().fit(Xd)                # las 64 componentes
acumulada = np.cumsum(completo.explained_variance_ratio_)
for meta in [0.5, 0.8, 0.9, 0.95]:
    k = np.argmax(acumulada >= meta) + 1  # primera que llega a la meta
    print(f"Para {meta:.0%} de la varianza: {k:2d} componentes")

k90 = np.argmax(acumulada >= 0.9) + 1
plt.plot(np.arange(1, 65), acumulada, color="C0")
plt.axhline(0.9, color="C2", linestyle="--", linewidth=1.5)
plt.scatter([k90], [acumulada[k90 - 1]], color="C1", s=120, zorder=3)
plt.xlabel("Número de componentes")
plt.ylabel("Varianza acumulada")
plt.title("¿Cuántas componentes hacen falta?")
plt.show()

from sklearn.manifold import TSNE
from sklearn.neighbors import NearestNeighbors
import umap

tsne = TSNE(n_components=2, perplexity=30, random_state=0)
Z_tsne = tsne.fit_transform(Xd)
reductor = umap.UMAP(n_neighbors=15, min_dist=0.1,
                     random_state=0, n_jobs=1)
Z_umap = reductor.fit_transform(Xd)

def vecino_igual(E, y):
    """Fracción de puntos cuyo vecino más cercano es el mismo dígito."""
    vecinos = NearestNeighbors(n_neighbors=2).fit(E)
    indices = vecinos.kneighbors(E, return_distance=False)[:, 1]
    return np.mean(y[indices] == y)

for nombre, E in [("64 dimensiones", Xd), ("PCA 2D", Z_sk),
                  ("t-SNE 2D", Z_tsne), ("UMAP 2D", Z_umap)]:
    print(f"{nombre:15s} {vecino_igual(E, yd):.1%}")

def dibujar_mapa(eje, E, titulo, destacado=1):
    """Los puntos en gris, uno de los dígitos en coral y las etiquetas."""
    es = yd == destacado
    eje.scatter(E[~es, 0], E[~es, 1], s=3, color="C2", alpha=0.35,
                edgecolors="none")
    eje.scatter(E[es, 0], E[es, 1], s=4, color="C1", edgecolors="none")
    for d in range(10):
        # La etiqueta va un poco arriba del centro de cada dígito
        centro = np.median(E[yd == d], axis=0)
        eje.annotate(str(d), centro, xytext=(0, 10),
                     textcoords="offset points", ha="center",
                     fontsize=14, fontweight="semibold", color="C0")
    eje.set_title(titulo)
    eje.set_xticks([])
    eje.set_yticks([])

fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
dibujar_mapa(ejes[0, 0], Z_sk, "PCA")
dibujar_mapa(ejes[0, 1], Z_tsne, "t-SNE")
dibujar_mapa(ejes[1, 0], Z_umap, "UMAP")
# t-SNE otra vez, ahora mirando solo 5 vecinos
Z_tsne5 = TSNE(n_components=2, perplexity=5,
               random_state=0).fit_transform(Xd)
dibujar_mapa(ejes[1, 1], Z_tsne5, "t-SNE, perplejidad 5")
plt.show()

rng = np.random.default_rng(1)
# Datos sintéticos en 10 dimensiones: un grupo apretado y uno disperso
apretado = rng.normal(0, 1, size=(200, 10))
disperso = rng.normal(0, 10, size=(200, 10)) + 100
A = np.vstack([apretado, disperso])
E = TSNE(n_components=2, random_state=0).fit_transform(A)

def ancho(P):
    """Distancia promedio de los puntos a su centro."""
    return np.mean(np.linalg.norm(P - P.mean(axis=0), axis=1))

print(f"En 10 dimensiones: el disperso es "
      f"{ancho(A[200:]) / ancho(A[:200]):.1f} veces más ancho")
print(f"En el mapa t-SNE:  el disperso es "
      f"{ancho(E[200:]) / ancho(E[:200]):.1f} veces más ancho")

from sklearn.datasets import fetch_olivetti_faces
from matplotlib.colors import LinearSegmentedColormap

caras = fetch_olivetti_faces()       # se descarga la primera vez
Xf = caras.data                      # 400 caras de 64 × 64 = 4096
print("Forma:", Xf.shape)

# Guardamos una foto de cada persona para probar al final
prueba = np.arange(len(Xf)) % 10 == 9
Xf_entr, Xf_prueba = Xf[~prueba], Xf[prueba]
print("Para ajustar:", Xf_entr.shape, "  para probar:", Xf_prueba.shape)

pca_caras = PCA(n_components=200, random_state=0).fit(Xf_entr)
acum = np.cumsum(pca_caras.explained_variance_ratio_)
for k in [5, 20, 100, 200]:
    print(f"{k:3d} componentes: {acum[k - 1]:.1%} de la varianza")

# Un mapa de color de la tinta al blanco (en Colab: "gray")
gris = LinearSegmentedColormap.from_list("gris", ["C0", "white"])
fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
imagenes = [pca_caras.mean_] + list(pca_caras.components_[:3])
titulos = ["Cara promedio", "Eigencara 1", "Eigencara 2", "Eigencara 3"]
for eje, img, titulo in zip(ejes.flat, imagenes, titulos):
    eje.imshow(img.reshape(64, 64), cmap=gris)
    eje.set_title(titulo, fontsize=14)
    eje.axis("off")
plt.show()

cara = Xf_prueba[[0]]            # una foto que el PCA no vio, forma (1, 4096)
fig, ejes = plt.subplots(2, 2, figsize=(4.8, 4.8))
ejes[0, 0].imshow(cara.reshape(64, 64), cmap=gris, vmin=0, vmax=1)
ejes[0, 0].set_title("Original: 4096", fontsize=14)
for eje, k in zip([ejes[0, 1], ejes[1, 0], ejes[1, 1]], [5, 20, 100]):
    W = pca_caras.components_[:k]       # las k primeras, forma (k, 4096)
    z = (cara - pca_caras.mean_) @ W.T  # k números por cara
    reconstruida = pca_caras.mean_ + z @ W
    error = np.mean((cara - reconstruida) ** 2)
    print(f"{k:3d} números: error cuadrático medio {error:.4f}")
    eje.imshow(reconstruida.reshape(64, 64), cmap=gris, vmin=0, vmax=1)
    eje.set_title(f"Con {k} números", fontsize=14)
for eje in ejes.flat:
    eje.axis("off")
plt.show()