Como desenhar um scree plot em python? [fechadas]

11

Estou usando decomposição de vetor singular em uma matriz e obtendo as matrizes U, S e Vt. Neste ponto, estou tentando escolher um limite para o número de dimensões a serem mantidas. Foi-me sugerido que olhasse para uma trama de seixos, mas estou imaginando como proceder para traçá-la numpy. Atualmente, estou fazendo o seguinte usando bibliotecas numpy e scipy em python:

U, S, Vt = svd(A)

Alguma sugestão?

lenda
fonte
1
pegue a diagonal de S, se já não é uma diagonal, quadrie-a, classifique-a em ordem decrescente, pegue a soma acumulada, divida pelo último valor e plote-a.
Shabbychef
@shabbychef: Você quer dizer, pegue a soma acumulada e divida pela soma de todos os valores, certo?
Legend
sim. No matlab, seria[U,S,V] = svd(X);S = cumsum(sort(diag(S).^2,1,'descend'));S = S ./ S(end);plot(S);
shabbychef

Respostas:

13

Aqui está um exemplo que pode ser colado em um prompt do IPython e gerar uma imagem como abaixo (ele usa dados aleatórios):

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()

insira a descrição da imagem aqui

Josh Hemann
fonte
Hermann: +1 Obrigado pelo seu tempo! Eu sei que tem sido um longo tempo, mas, no entanto, isso é muito bom ter :)
Legend
o que é num_vars? parece não estar definido no seu script.
TheChymera
@TheChymera - Obrigado por capturar isso, eu atualizei minha resposta.
Josh Hemann 12/08
@ Josh Hemann Sim, eu também descobri isso no médio prazo - mas eu acho que poderia ser melhor para calculá-lo a partir da forma de A
TheChymera
1
@ JoshHemann, você pode explicar isso: eigvals = S ** 2 / np.cumsum (S) [- 1] ?? I ter visto a partir de alguns documentos de eigvals = S ** 2 / (n-1), onde n é o número de recursos
MAKIS