python - Pythonの主成分分析（PCA）

Question

（26424 x 144）配列があり、Pythonを使用してその配列に対してPCAを実行したいと思います。ただし、このタスクを実行する方法を説明するWeb上の特定の場所はありません（独自にPCAを実行するサイトがいくつかあります。私が見つけることができるように実行する一般的な方法はありません）。どんな種類の助けもあれば誰でも素晴らしいことをするでしょう。

score 99 · Accepted Answer

別の回答がすでに受け入れられているにもかかわらず、私は自分の回答を投稿しました。受け入れられた回答は、非推奨の関数に依存しています。さらに、この非推奨の関数は、特異値分解（SVD）に基づいています。これは、（完全に有効ですが）PCAを計算するための2つの一般的な手法の中ではるかにメモリとプロセッサを消費します。OP内のデータ配列のサイズのため、これはここで特に関係があります。共分散ベースのPCAを使用すると、計算フローで使用される配列は、 26424 x 144 （元のデータ配列の次元）ではなく、 144x144になります。

これは、 SciPyのlinalgモジュールを使用したPCAの簡単な実装です。この実装は、最初に共分散行列を計算し、次にこの配列に対して後続のすべての計算を実行するため、SVDベースのPCAよりもはるかに少ないメモリを使用します。

（ NumPyのlinalgモジュールは、LAとしてのnumpy import linalgからのimportステートメントを除いて、以下のコードを変更せずに使用することもできます。）

このPCA実装の2つの重要なステップは次のとおりです。

共分散行列の計算; と
このcov行列の固有ベクトルと固有値を取得します

以下の関数では、パラメーターdims_rescaled_dataは、再スケーリングされたデータマトリックス内の必要な次元数を参照します。このパラメーターのデフォルト値は2次元ですが、以下のコードは2次元に制限されていませんが、元のデータ配列の列番号よりも小さい値にすることができます。

def PCA(data, dims_rescaled_data=2):
    """
    returns: data transformed in 2 dims/columns + regenerated original data
    pass in: data as 2D NumPy array
    """
    import numpy as NP
    from scipy import linalg as LA
    m, n = data.shape
    # mean center the data
    data -= data.mean(axis=0)
    # calculate the covariance matrix
    R = NP.cov(data, rowvar=False)
    # calculate eigenvectors & eigenvalues of the covariance matrix
    # use 'eigh' rather than 'eig' since R is symmetric, 
    # the performance gain is substantial
    evals, evecs = LA.eigh(R)
    # sort eigenvalue in decreasing order
    idx = NP.argsort(evals)[::-1]
    evecs = evecs[:,idx]
    # sort eigenvectors according to same index
    evals = evals[idx]
    # select the first n eigenvectors (n is desired dimension
    # of rescaled data array, or dims_rescaled_data)
    evecs = evecs[:, :dims_rescaled_data]
    # carry out the transformation on the data using eigenvectors
    # and return the re-scaled data, eigenvalues, and eigenvectors
    return NP.dot(evecs.T, data.T).T, evals, evecs

def test_PCA(data, dims_rescaled_data=2):
    '''
    test by attempting to recover original data array from
    the eigenvectors of its covariance matrix & comparing that
    'recovered' array with the original data
    '''
    _ , _ , eigenvectors = PCA(data, dim_rescaled_data=2)
    data_recovered = NP.dot(eigenvectors, m).T
    data_recovered += data_recovered.mean(axis=0)
    assert NP.allclose(data, data_recovered)
    

def plot_pca(data):
    from matplotlib import pyplot as MPL
    clr1 =  '#2026B2'
    fig = MPL.figure()
    ax1 = fig.add_subplot(111)
    data_resc, data_orig = PCA(data)
    ax1.plot(data_resc[:, 0], data_resc[:, 1], '.', mfc=clr1, mec=clr1)
    MPL.show()

>>> # iris, probably the most widely used reference data set in ML
>>> df = "~/iris.csv"
>>> data = NP.loadtxt(df, delimiter=',')
>>> # remove class labels
>>> data = data[:,:-1]
>>> plot_pca(data)

以下のプロットは、虹彩データでのこのPCA関数の視覚的表現です。ご覧のとおり、2D変換は、クラスIをクラスIIおよびクラスIIIから明確に分離します（ただし、クラスIIをクラスIIIから分離することはできません。実際、別の次元が必要です）。

ここに画像の説明を入力してください

score 62 · Accepted Answer

PCA関数はmatplotlibモジュールにあります。

import numpy as np
from matplotlib.mlab import PCA

data = np.array(np.random.randint(10,size=(10,3)))
results = PCA(data)

結果には、PCAのさまざまなパラメーターが格納されます。これは、MATLAB構文との互換性レイヤーであるmatplotlibのmlab部分からのものです。

編集：ブログnextgeneticsで、matplotlib mlabモジュールを使用してPCAを実行および表示する方法のすばらしいデモンストレーションを見つけました。楽しんで、そのブログを確認してください。

score 30 · Accepted Answer

numpyを使用する別のPythonPCA。@dougと同じアイデアですが、実行されませんでした。

from numpy import array, dot, mean, std, empty, argsort
from numpy.linalg import eigh, solve
from numpy.random import randn
from matplotlib.pyplot import subplots, show

def cov(X):
    """
    Covariance matrix
    note: specifically for mean-centered data
    note: numpy's `cov` uses N-1 as normalization
    """
    return dot(X.T, X) / X.shape[0]
    # N = data.shape[1]
    # C = empty((N, N))
    # for j in range(N):
    #   C[j, j] = mean(data[:, j] * data[:, j])
    #   for k in range(j + 1, N):
    #       C[j, k] = C[k, j] = mean(data[:, j] * data[:, k])
    # return C

def pca(data, pc_count = None):
    """
    Principal component analysis using eigenvalues
    note: this mean-centers and auto-scales the data (in-place)
    """
    data -= mean(data, 0)
    data /= std(data, 0)
    C = cov(data)
    E, V = eigh(C)
    key = argsort(E)[::-1][:pc_count]
    E, V = E[key], V[:, key]
    U = dot(data, V)  # used to be dot(V.T, data.T).T
    return U, E, V

""" test data """
data = array([randn(8) for k in range(150)])
data[:50, 2:4] += 5
data[50:, 2:5] += 5

""" visualize """
trans = pca(data, 3)[0]
fig, (ax1, ax2) = subplots(1, 2)
ax1.scatter(data[:50, 0], data[:50, 1], c = 'r')
ax1.scatter(data[50:, 0], data[50:, 1], c = 'b')
ax2.scatter(trans[:50, 0], trans[:50, 1], c = 'r')
ax2.scatter(trans[50:, 0], trans[50:, 1], c = 'b')
show()

これは、はるかに短いものと同じものをもたらします

from sklearn.decomposition import PCA

def pca2(data, pc_count = None):
    return PCA(n_components = 4).fit_transform(data)

私が理解しているように、固有値（最初の方法）を使用すると、高次元のデータとサンプルが少ない場合に適していますが、特異値分解を使用する場合は、次元よりもサンプルが多い場合に適しています。

score 15 · Accepted Answer

これはの仕事ですnumpy。

numpyそして、これは、のような組み込みモジュールを使用して、主要なコンポーネント分析を実行する方法を示すチュートリアルmean,cov,double,cumsum,dot,linalg,array,rankです。

http://glowingpython.blogspot.sg/2011/07/principal-component-analysis-with-numpy.html

scipyここにも長い説明があることに注意してください-https://github.com/scikit-learn/scikit-learn/blob/babe4a5d0637ca172d47e1dfdd2f6f3c3ecb28db/scikits/learn/utils/extmath.py#L105

scikit-learnより多くのコード例があるライブラリで-https://github.com/scikit-learn/scikit-learn/blob/babe4a5d0637ca172d47e1dfdd2f6f3c3ecb28db/scikits/learn/utils/extmath.py#L105

score 7 · Accepted Answer

scikit-learnオプションは次のとおりです。どちらの方法でも、 PCAはスケールの影響を受けるため、StandardScalerが使用されました

方法1：分散の少なくともx％（以下の例では90％）が保持されるように、scikit-learnに主成分の最小数を選択させます。

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

iris = load_iris()

# mean-centers and auto-scales the data
standardizedData = StandardScaler().fit_transform(iris.data)

pca = PCA(.90)

principalComponents = pca.fit_transform(X = standardizedData)

# To get how many principal components was chosen
print(pca.n_components_)

方法2：主成分の数を選択します（この場合、2つを選択しました）

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

iris = load_iris()

standardizedData = StandardScaler().fit_transform(iris.data)

pca = PCA(n_components=2)

principalComponents = pca.fit_transform(X = standardizedData)

# to get how much variance was retained
print(pca.explained_variance_ratio_.sum())

ソース：https ：//towardsdatascience.com/pca-using-python-scikit-learn-e653f8989e60

score 4 · Accepted Answer

更新： matplotlib.mlab.PCAリリース2.2（2018-03-06）以降、実際に非推奨になりました。

ライブラリmatplotlib.mlab.PCA（この回答で使用）は非推奨ではありません。そこで、Google経由でここに到着するすべての人のために、Python2.7でテストされた完全な実例を投稿します。

次のコードは非推奨のライブラリを使用しているため、注意して使用してください。

from matplotlib.mlab import PCA
import numpy
data = numpy.array( [[3,2,5], [-2,1,6], [-1,0,4], [4,3,4], [10,-5,-6]] )
pca = PCA(data)

現在、 `pca.Y'には、主成分基底ベクトルに関する元のデータ行列があります。PCAオブジェクトの詳細については、こちらをご覧ください。

>>> pca.Y
array([[ 0.67629162, -0.49384752,  0.14489202],
   [ 1.26314784,  0.60164795,  0.02858026],
   [ 0.64937611,  0.69057287, -0.06833576],
   [ 0.60697227, -0.90088738, -0.11194732],
   [-3.19578784,  0.10251408,  0.00681079]])

matplotlib.pyplotPCAが「良好な」結果をもたらすことを自分自身に納得させるために、このデータを描画するために使用できます。このnamesリストは、5つのベクトルに注釈を付けるために使用されています。

import matplotlib.pyplot
names = [ "A", "B", "C", "D", "E" ]
matplotlib.pyplot.scatter(pca.Y[:,0], pca.Y[:,1])
for label, x, y in zip(names, pca.Y[:,0], pca.Y[:,1]):
    matplotlib.pyplot.annotate( label, xy=(x, y), xytext=(-2, 2), textcoords='offset points', ha='right', va='bottom' )
matplotlib.pyplot.show()

元のベクトルを見ると、data [0]（ "A"）とdata [3]（ "D"）は、data [1]（ "B"）とdata [2]（ " C "）。これは、PCA変換されたデータの2Dプロットに反映されます。

score 4 · Accepted Answer

他のすべての回答に加えて、`biplot`using`sklearn`とをプロットするためのコードを次に示し`matplotlib`ます。

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.decomposition import PCA
import pandas as pd
from sklearn.preprocessing import StandardScaler

iris = datasets.load_iris()
X = iris.data
y = iris.target
#In general a good idea is to scale the data
scaler = StandardScaler()
scaler.fit(X)
X=scaler.transform(X)    

pca = PCA()
x_new = pca.fit_transform(X)

def myplot(score,coeff,labels=None):
    xs = score[:,0]
    ys = score[:,1]
    n = coeff.shape[0]
    scalex = 1.0/(xs.max() - xs.min())
    scaley = 1.0/(ys.max() - ys.min())
    plt.scatter(xs * scalex,ys * scaley, c = y)
    for i in range(n):
        plt.arrow(0, 0, coeff[i,0], coeff[i,1],color = 'r',alpha = 0.5)
        if labels is None:
            plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, "Var"+str(i+1), color = 'g', ha = 'center', va = 'center')
        else:
            plt.text(coeff[i,0]* 1.15, coeff[i,1] * 1.15, labels[i], color = 'g', ha = 'center', va = 'center')
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.xlabel("PC{}".format(1))
plt.ylabel("PC{}".format(2))
plt.grid()

#Call the function. Use only the 2 PCs.
myplot(x_new[:,0:2],np.transpose(pca.components_[0:2, :]))
plt.show()

score 2 · Accepted Answer

さまざまなPCAを比較するための小さなスクリプトを作成しました。ここに回答として表示されます。

import numpy as np
from scipy.linalg import svd

shape = (26424, 144)
repeat = 20
pca_components = 2

data = np.array(np.random.randint(255, size=shape)).astype('float64')

# data normalization
# data.dot(data.T)
# (U, s, Va) = svd(data, full_matrices=False)
# data = data / s[0]

from fbpca import diffsnorm
from timeit import default_timer as timer

from scipy.linalg import svd
start = timer()
for i in range(repeat):
    (U, s, Va) = svd(data, full_matrices=False)
time = timer() - start
err = diffsnorm(data, U, s, Va)
print('svd time: %.3fms, error: %E' % (time*1000/repeat, err))


from matplotlib.mlab import PCA
start = timer()
_pca = PCA(data)
for i in range(repeat):
    U = _pca.project(data)
time = timer() - start
err = diffsnorm(data, U, _pca.fracs, _pca.Wt)
print('matplotlib PCA time: %.3fms, error: %E' % (time*1000/repeat, err))

from fbpca import pca
start = timer()
for i in range(repeat):
    (U, s, Va) = pca(data, pca_components, True)
time = timer() - start
err = diffsnorm(data, U, s, Va)
print('facebook pca time: %.3fms, error: %E' % (time*1000/repeat, err))


from sklearn.decomposition import PCA
start = timer()
_pca = PCA(n_components = pca_components)
_pca.fit(data)
for i in range(repeat):
    U = _pca.transform(data)
time = timer() - start
err = diffsnorm(data, U, _pca.explained_variance_, _pca.components_)
print('sklearn PCA time: %.3fms, error: %E' % (time*1000/repeat, err))

start = timer()
for i in range(repeat):
    (U, s, Va) = pca_mark(data, pca_components)
time = timer() - start
err = diffsnorm(data, U, s, Va.T)
print('pca by Mark time: %.3fms, error: %E' % (time*1000/repeat, err))

start = timer()
for i in range(repeat):
    (U, s, Va) = pca_doug(data, pca_components)
time = timer() - start
err = diffsnorm(data, U, s[:pca_components], Va.T)
print('pca by doug time: %.3fms, error: %E' % (time*1000/repeat, err))

pca_markは、Markの回答のpcaです。

pca_dougは、ダグの回答のpcaです。

出力例を次に示します（ただし、結果はデータサイズとpca_componentsに大きく依存するため、独自のデータを使用して独自のテストを実行することをお勧めします。また、Facebookのpcaは正規化されたデータ用に最適化されているため、より高速になり、その場合はより正確になります）：

svd time: 3212.228ms, error: 1.907320E-10
matplotlib PCA time: 879.210ms, error: 2.478853E+05
facebook pca time: 485.483ms, error: 1.260335E+04
sklearn PCA time: 169.832ms, error: 7.469847E+07
pca by Mark time: 293.758ms, error: 1.713129E+02
pca by doug time: 300.326ms, error: 1.707492E+02

編集：

fbpcaのdiffsnorm関数は、シュール分解のスペクトルノルム誤差を計算します。

score 0 · Accepted Answer

動作するためdef plot_pca(data):に、行を交換する必要があります

data_resc, data_orig = PCA(data)
ax1.plot(data_resc[:, 0], data_resc[:, 1], '.', mfc=clr1, mec=clr1)

線で

newData, data_resc, data_orig = PCA(data)
ax1.plot(newData[:, 0], newData[:, 1], '.', mfc=clr1, mec=clr1)

score 0 · Accepted Answer

このサンプルコードは、日本のイールドカーブをロードし、PCAコンポーネントを作成します。次に、PCAを使用して特定の日付の移動を推定し、実際の移動と比較します。

%matplotlib inline

import numpy as np
import scipy as sc
from scipy import stats
from IPython.display import display, HTML
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import datetime
from datetime import timedelta

import quandl as ql

start = "2016-10-04"
end = "2019-10-04"

ql_data = ql.get("MOFJ/INTEREST_RATE_JAPAN", start_date = start, end_date = end).sort_index(ascending= False)

eigVal_, eigVec_ = np.linalg.eig(((ql_data[:300]).diff(-1)*100).cov()) # take latest 300 data-rows and normalize to bp
print('number of PCA are', len(eigVal_))

loc_ = 10
plt.plot(eigVec_[:,0], label = 'PCA1')
plt.plot(eigVec_[:,1], label = 'PCA2')
plt.plot(eigVec_[:,2], label = 'PCA3')
plt.xticks(range(len(eigVec_[:,0])), ql_data.columns)
plt.legend()
plt.show()

x = ql_data.diff(-1).iloc[loc_].values * 100 # set the differences
x_ = x[:,np.newaxis]
a1, _, _, _ = np.linalg.lstsq(eigVec_[:,0][:, np.newaxis], x_) # linear regression without intercept
a2, _, _, _ = np.linalg.lstsq(eigVec_[:,1][:, np.newaxis], x_)
a3, _, _, _ = np.linalg.lstsq(eigVec_[:,2][:, np.newaxis], x_)

pca_mv = m1 * eigVec_[:,0] + m2 * eigVec_[:,1] + m3 * eigVec_[:,2] + c1 + c2 + c3
pca_MV = a1[0][0] * eigVec_[:,0] + a2[0][0] * eigVec_[:,1] + a3[0][0] * eigVec_[:,2]
pca_mV = b1 * eigVec_[:,0] + b2 * eigVec_[:,1] + b3 * eigVec_[:,2]

display(pd.DataFrame([eigVec_[:,0], eigVec_[:,1], eigVec_[:,2], x, pca_MV]))
print('PCA1 regression is', a1, a2, a3)


plt.plot(pca_MV)
plt.title('this is with regression and no intercept')
plt.plot(ql_data.diff(-1).iloc[loc_].values * 100, )
plt.title('this is with actual moves')
plt.show()

score 0 · Accepted Answer

これは、簡単に理解できる手順を含め、PCAで見つけることができる最も簡単な答えかもしれません。最大の情報を提供する144から2つの主要な次元を保持したいとします。

まず、2次元配列をデータフレームに変換します。

import pandas as pd

# Here X is your array of size (26424 x 144)
data = pd.DataFrame(X)

次に、2つの方法があります。

方法1：手動計算

ステップ1：Xに列の標準化を適用する

from sklearn import preprocessing

scalar = preprocessing.StandardScaler()
standardized_data = scalar.fit_transform(data)

ステップ2：元の行列Xの共分散行列Sを見つける

sample_data = standardized_data
covar_matrix = np.cov(sample_data)

ステップ3：Sの固有値と固有ベクトルを見つけます（ここでは2Dなので、それぞれ2つ）

from scipy.linalg import eigh

# eigh() function will provide eigen-values and eigen-vectors for a given matrix.
# eigvals=(low value, high value) takes eigen value numbers in ascending order
values, vectors = eigh(covar_matrix, eigvals=(142,143))

# Converting the eigen vectors into (2,d) shape for easyness of further computations
vectors = vectors.T

ステップ4：データを変換する

# Projecting the original data sample on the plane formed by two principal eigen vectors by vector-vector multiplication.

new_coordinates = np.matmul(vectors, sample_data.T)
print(new_coordinates.T)

これnew_coordinates.Tは、2つの主成分を持つサイズ（26424 x 2）になります。

方法2：Scikit-Learnを使用する

ステップ1：Xに列の標準化を適用する

from sklearn import preprocessing

scalar = preprocessing.StandardScaler()
standardized_data = scalar.fit_transform(data)

ステップ2：pcaを初期化する

from sklearn import decomposition

# n_components = numbers of dimenstions you want to retain
pca = decomposition.PCA(n_components=2)

ステップ3：pcaを使用してデータを適合させる

# This line takes care of calculating co-variance matrix, eigen values, eigen vectors and multiplying top 2 eigen vectors with data-matrix X.
pca_data = pca.fit_transform(sample_data)

これpca_dataは、2つの主成分を持つサイズ（26424 x 2）になります。

python - Pythonの主成分分析（PCA）

11 に答える 11

他のすべての回答に加えて、biplotusingsklearnとをプロットするためのコードを次に示しmatplotlibます。

Related

Reference

他のすべての回答に加えて、`biplot`using`sklearn`とをプロットするためのコードを次に示し`matplotlib`ます。