Comment dessiner un tracé éboulis en python? [fermé]

11

J'utilise la décomposition vectorielle singulière sur une matrice et j'obtiens les matrices U, S et Vt. À ce stade, j'essaie de choisir un seuil pour le nombre de dimensions à conserver. On m'a suggéré de regarder un tracé d'éboulis, mais je me demande comment s'y prendre pour le tracer en numpy. Actuellement, je fais ce qui suit en utilisant les bibliothèques numpy et scipy en python:

U, S, Vt = svd(A)

Aucune suggestion?

Légende
la source
1
prendre la diagonale de S, si ce n'est pas déjà une diagonale, la mettre au carré, la trier par ordre décroissant, prendre la somme cumulée, diviser par la dernière valeur, puis la tracer.
shabbychef
@shabbychef: Vous voulez dire, prenez la somme cumulée et divisez par la somme de toutes les valeurs, n'est-ce pas?
Legend
Oui. En matlab, ce serait[U,S,V] = svd(X);S = cumsum(sort(diag(S).^2,1,'descend'));S = S ./ S(end);plot(S);
shabbychef

Réponses:

13

Voici un exemple qui peut être collé dans une invite IPython et générer une image comme ci-dessous (il utilise des données aléatoires):

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

#Make a random array and then make it positive-definite
num_vars = 6
num_obs = 9
A = np.random.randn(num_obs, num_vars)
A = np.asmatrix(A.T) * np.asmatrix(A)
U, S, V = np.linalg.svd(A) 
eigvals = S**2 / np.sum(S**2)  # NOTE (@amoeba): These are not PCA eigenvalues. 
                               # This question is about SVD.

fig = plt.figure(figsize=(8,5))
sing_vals = np.arange(num_vars) + 1
plt.plot(sing_vals, eigvals, 'ro-', linewidth=2)
plt.title('Scree Plot')
plt.xlabel('Principal Component')
plt.ylabel('Eigenvalue')
#I don't like the default legend so I typically make mine like below, e.g.
#with smaller fonts and a bit transparent so I do not cover up data, and make
#it moveable by the viewer in case upper-right is a bad place for it 
leg = plt.legend(['Eigenvalues from SVD'], loc='best', borderpad=0.3, 
                 shadow=False, prop=matplotlib.font_manager.FontProperties(size='small'),
                 markerscale=0.4)
leg.get_frame().set_alpha(0.4)
leg.draggable(state=True)
plt.show()

entrez la description de l'image ici

Josh Hemann
la source
Hermann: +1 Merci pour votre temps! Je sais que ça fait longtemps mais néanmoins c'est vraiment bien d'avoir :)
Legend
c'est quoi num_vars? il ne semble pas être défini dans votre script.
TheChymera
@TheChymera - Merci d'avoir attrapé cela, j'ai mis à jour ma réponse.
Josh Hemann
@Josh Hemann oui, j'ai aussi compris cela entre-temps - mais je pense qu'il serait préférable de le calculer à partir de la forme de A
TheChymera
1
@JoshHemann pouvez-vous expliquer cela: eigvals = S ** 2 / np.cumsum (S) [- 1] ?? J'ai vu sur la base de certains papiers des eigvals = S ** 2 / (n-1) où n est le nombre de fonctionnalités
makis