주성분 분석에서는 파이썬
https://stackoverflow.com/questions/1730600
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian문제
을 사용하고 싶 PCA(principal component analysis)에 대한 차원 감소입니다.가 numpy 또는 scipy 이미 그것이 있거나,나는 나 자신을 사용하여 numpy.linalg.eigh?
나지 않아 그냥 사용하려는 단일 값이 분해(SVD)기 때문에 입력된 데이터는 매우 높은차원(~460 크기),그래서 내가 생각하 SVD 느린 것보다는 컴퓨팅의 고유 벡터의 공분산 행렬을 선택할 수 있습니다.
을 찾아 초록색,파란색,빨간색,디버깅 구현는 이미 대해 올바른 결정을 내릴 때 사용하는 방법,그리고 어쩌면 다른 최적화하는지 모르겠다.
해결책
당신은 볼 수 있습니다 MDP.
직접 테스트 할 기회가 없었지만 PCA 기능을 위해 정확히 북마크했습니다.
다른 팁
몇 달 후, 여기에 작은 클래스 PCA와 사진이 있습니다.
#!/usr/bin/env python
""" a small class for Principal Component Analysis
Usage:
p = PCA( A, fraction=0.90 )
In:
A: an array of e.g. 1000 observations x 20 variables, 1000 rows x 20 columns
fraction: use principal components that account for e.g.
90 % of the total variance
Out:
p.U, p.d, p.Vt: from numpy.linalg.svd, A = U . d . Vt
p.dinv: 1/d or 0, see NR
p.eigen: the eigenvalues of A*A, in decreasing order (p.d**2).
eigen[j] / eigen.sum() is variable j's fraction of the total variance;
look at the first few eigen[] to see how many PCs get to 90 %, 95 % ...
p.npc: number of principal components,
e.g. 2 if the top 2 eigenvalues are >= `fraction` of the total.
It's ok to change this; methods use the current value.
Methods:
The methods of class PCA transform vectors or arrays of e.g.
20 variables, 2 principal components and 1000 observations,
using partial matrices U' d' Vt', parts of the full U d Vt:
A ~ U' . d' . Vt' where e.g.
U' is 1000 x 2
d' is diag([ d0, d1 ]), the 2 largest singular values
Vt' is 2 x 20. Dropping the primes,
d . Vt 2 principal vars = p.vars_pc( 20 vars )
U 1000 obs = p.pc_obs( 2 principal vars )
U . d . Vt 1000 obs, p.obs( 20 vars ) = pc_obs( vars_pc( vars ))
fast approximate A . vars, using the `npc` principal components
Ut 2 pcs = p.obs_pc( 1000 obs )
V . dinv 20 vars = p.pc_vars( 2 principal vars )
V . dinv . Ut 20 vars, p.vars( 1000 obs ) = pc_vars( obs_pc( obs )),
fast approximate Ainverse . obs: vars that give ~ those obs.
Notes:
PCA does not center or scale A; you usually want to first
A -= A.mean(A, axis=0)
A /= A.std(A, axis=0)
with the little class Center or the like, below.
See also:
http://en.wikipedia.org/wiki/Principal_component_analysis
http://en.wikipedia.org/wiki/Singular_value_decomposition
Press et al., Numerical Recipes (2 or 3 ed), SVD
PCA micro-tutorial
iris-pca .py .png
"""
from __future__ import division
import numpy as np
dot = np.dot
# import bz.numpyutil as nu
# dot = nu.pdot
__version__ = "2010-04-14 apr"
__author_email__ = "denis-bz-py at t-online dot de"
#...............................................................................
class PCA:
def __init__( self, A, fraction=0.90 ):
assert 0 <= fraction <= 1
# A = U . diag(d) . Vt, O( m n^2 ), lapack_lite --
self.U, self.d, self.Vt = np.linalg.svd( A, full_matrices=False )
assert np.all( self.d[:-1] >= self.d[1:] ) # sorted
self.eigen = self.d**2
self.sumvariance = np.cumsum(self.eigen)
self.sumvariance /= self.sumvariance[-1]
self.npc = np.searchsorted( self.sumvariance, fraction ) + 1
self.dinv = np.array([ 1/d if d > self.d[0] * 1e-6 else 0
for d in self.d ])
def pc( self ):
""" e.g. 1000 x 2 U[:, :npc] * d[:npc], to plot etc. """
n = self.npc
return self.U[:, :n] * self.d[:n]
# These 1-line methods may not be worth the bother;
# then use U d Vt directly --
def vars_pc( self, x ):
n = self.npc
return self.d[:n] * dot( self.Vt[:n], x.T ).T # 20 vars -> 2 principal
def pc_vars( self, p ):
n = self.npc
return dot( self.Vt[:n].T, (self.dinv[:n] * p).T ) .T # 2 PC -> 20 vars
def pc_obs( self, p ):
n = self.npc
return dot( self.U[:, :n], p.T ) # 2 principal -> 1000 obs
def obs_pc( self, obs ):
n = self.npc
return dot( self.U[:, :n].T, obs ) .T # 1000 obs -> 2 principal
def obs( self, x ):
return self.pc_obs( self.vars_pc(x) ) # 20 vars -> 2 principal -> 1000 obs
def vars( self, obs ):
return self.pc_vars( self.obs_pc(obs) ) # 1000 obs -> 2 principal -> 20 vars
class Center:
""" A -= A.mean() /= A.std(), inplace -- use A.copy() if need be
uncenter(x) == original A . x
"""
# mttiw
def __init__( self, A, axis=0, scale=True, verbose=1 ):
self.mean = A.mean(axis=axis)
if verbose:
print "Center -= A.mean:", self.mean
A -= self.mean
if scale:
std = A.std(axis=axis)
self.std = np.where( std, std, 1. )
if verbose:
print "Center /= A.std:", self.std
A /= self.std
else:
self.std = np.ones( A.shape[-1] )
self.A = A
def uncenter( self, x ):
return np.dot( self.A, x * self.std ) + np.dot( x, self.mean )
#...............................................................................
if __name__ == "__main__":
import sys
csv = "iris4.csv" # wikipedia Iris_flower_data_set
# 5.1,3.5,1.4,0.2 # ,Iris-setosa ...
N = 1000
K = 20
fraction = .90
seed = 1
exec "\n".join( sys.argv[1:] ) # N= ...
np.random.seed(seed)
np.set_printoptions( 1, threshold=100, suppress=True ) # .1f
try:
A = np.genfromtxt( csv, delimiter="," )
N, K = A.shape
except IOError:
A = np.random.normal( size=(N, K) ) # gen correlated ?
print "csv: %s N: %d K: %d fraction: %.2g" % (csv, N, K, fraction)
Center(A)
print "A:", A
print "PCA ..." ,
p = PCA( A, fraction=fraction )
print "npc:", p.npc
print "% variance:", p.sumvariance * 100
print "Vt[0], weights that give PC 0:", p.Vt[0]
print "A . Vt[0]:", dot( A, p.Vt[0] )
print "pc:", p.pc()
print "\nobs <-> pc <-> x: with fraction=1, diffs should be ~ 0"
x = np.ones(K)
# x = np.ones(( 3, K ))
print "x:", x
pc = p.vars_pc(x) # d' Vt' x
print "vars_pc(x):", pc
print "back to ~ x:", p.pc_vars(pc)
Ax = dot( A, x.T )
pcx = p.obs(x) # U' d' Vt' x
print "Ax:", Ax
print "A'x:", pcx
print "max |Ax - A'x|: %.2g" % np.linalg.norm( Ax - pcx, np.inf )
b = Ax # ~ back to original x, Ainv A x
back = p.vars(b)
print "~ back again:", back
print "max |back - x|: %.2g" % np.linalg.norm( back - x, np.inf )
# end pca.py
PCA 사용 numpy.linalg.svd 매우 쉽습니다. 간단한 데모는 다음과 같습니다.
import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import lena
# the underlying signal is a sinusoidally modulated image
img = lena()
t = np.arange(100)
time = np.sin(0.1*t)
real = time[:,np.newaxis,np.newaxis] * img[np.newaxis,...]
# we add some noise
noisy = real + np.random.randn(*real.shape)*255
# (observations, features) matrix
M = noisy.reshape(noisy.shape[0],-1)
# singular value decomposition factorises your data matrix such that:
#
# M = U*S*V.T (where '*' is matrix multiplication)
#
# * U and V are the singular matrices, containing orthogonal vectors of
# unit length in their rows and columns respectively.
#
# * S is a diagonal matrix containing the singular values of M - these
# values squared divided by the number of observations will give the
# variance explained by each PC.
#
# * if M is considered to be an (observations, features) matrix, the PCs
# themselves would correspond to the rows of S^(1/2)*V.T. if M is
# (features, observations) then the PCs would be the columns of
# U*S^(1/2).
#
# * since U and V both contain orthonormal vectors, U*V.T is equivalent
# to a whitened version of M.
U, s, Vt = np.linalg.svd(M, full_matrices=False)
V = Vt.T
# PCs are already sorted by descending order
# of the singular values (i.e. by the
# proportion of total variance they explain)
# if we use all of the PCs we can reconstruct the noisy signal perfectly
S = np.diag(s)
Mhat = np.dot(U, np.dot(S, V.T))
print "Using all PCs, MSE = %.6G" %(np.mean((M - Mhat)**2))
# if we use only the first 20 PCs the reconstruction is less accurate
Mhat2 = np.dot(U[:, :20], np.dot(S[:20, :20], V[:,:20].T))
print "Using first 20 PCs, MSE = %.6G" %(np.mean((M - Mhat2)**2))
fig, [ax1, ax2, ax3] = plt.subplots(1, 3)
ax1.imshow(img)
ax1.set_title('true image')
ax2.imshow(noisy.mean(0))
ax2.set_title('mean of noisy images')
ax3.imshow((s[0]**(1./2) * V[:,0]).reshape(img.shape))
ax3.set_title('first spatial PC')
plt.show()
Sklearn을 사용할 수 있습니다.
import sklearn.decomposition as deco
import numpy as np
x = (x - np.mean(x, 0)) / np.std(x, 0) # You need to normalize your data first
pca = deco.PCA(n_components) # n_components is the components number after reduction
x_r = pca.fit(x).transform(x)
print ('explained variance (first %d components): %.2f'%(n_components, sum(pca.explained_variance_ratio_)))
SVD는 460 차원에서 정상적으로 작동해야합니다. Atom Netbook에서 약 7 초가 걸립니다. EIG () 메소드가 사용됩니다 더 시간 (더 많은 부동 소수점 조작을 사용해야 함)은 거의 항상 덜 정확합니다.
460 가지 예제가 미만인 경우 원하는 것은 산란 매트릭스 (x -datamean)^t (x -mean)를 대각선화하고, 데이터 포인트가 열이라고 가정하고 (x -datamean)에 의해 왼쪽 다중화됩니다. 저것 ~할 것 같다 데이터보다 차원이 더 많은 경우 더 빨리하십시오.
당신은 자신의 사용을 쉽게 "롤"할 수 있습니다. scipy.linalg (사전 중심 데이터 세트를 가정합니다 data):
covmat = data.dot(data.T)
evs, evmat = scipy.linalg.eig(covmat)
그 다음에 evs 당신의 고유 값입니다 evmat 투영 행렬입니다.
당신이 유지하고 싶다면 d 치수, 첫 번째를 사용하십시오 d 고유 값과 먼저 d 고유 벡터.
을 고려하면 scipy.linalg 분해와 Numpy가 매트릭스 곱셈이 있습니까?
나는 그냥 읽은 책 Machine Learning:는 알고리즘적 관점.모든 코드 예제에서는 이 책에 의해 작성되었 Python(및 거 Numpy).코드 조각의 chatper10.2 주요 구성 요소 분석 어쩌면 가치가 읽고 있습니다.그것은 사용 numpy.linalg.eig.
방법으로,생각 SVD 처리할 수 있습 460*460 차원습니다.내 계산 6500*6500SVD 와 numpy/scipy.linalg.svd 에서 아주 오래된 PC:펜티엄 III733mHz.솔직히 말해서,스크립트가 필요 많은 메모리(에 대해 1.xG)고 많은 시간(30 분)을 얻을 수 SVD 결과입니다.그러나 내가 생각하 460*460 에 현대 PC 하지 않는 것은 큰 문제가 될지도가 필요로 하 SVD 의 거대한 숫자 번입니다.
IT에서 완전한 단수 값 분해 (SVD)가 필요하지 않습니다. 모든 고유 값과 고유 벡터를 계산하고 큰 행렬에 적합 할 수 있습니다.Scipy 그리고 스파 스 모듈은 스파 스와 밀도가 높은 행렬 모두에서 작동하는 일반적인 선형 알그레브라 함수를 제공하며, 그중에는 EIG* 기능이 있습니다.
http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#matrix-factorizations
Scikit-Learn a Python PCA 구현 지금은 조밀 한 행렬 만 지원합니다.
타이밍 :
In [1]: A = np.random.randn(1000, 1000)
In [2]: %timeit scipy.sparse.linalg.eigsh(A)
1 loops, best of 3: 802 ms per loop
In [3]: %timeit np.linalg.svd(A)
1 loops, best of 3: 5.91 s per loop
여기 Numpy, Scipy 및 C-Extensions를 사용하여 Python 용 PCA 모듈의 또 다른 구현입니다. 이 모듈은 SVD 또는 Nipals (비선형 반복 부분 최소 제곱) 알고리즘을 사용하여 PCA를 운반합니다.
3D 벡터로 작업하는 경우 도구 벨트를 사용하여 SVD를 간결하게 적용 할 수 있습니다. VG. Numpy 위에있는 가벼운 층입니다.
import numpy as np
import vg
vg.principal_components(data)
첫 번째 주요 구성 요소 만 원한다면 편리한 별칭도 있습니다.
vg.major_axis(data)
나는 마지막 스타트 업에서 도서관을 만들었습니다. 여기서는 다음과 같은 용도로 동기를 부여했습니다.
제휴하지 않습니다 StackOverflow
