Анализ основных компонентов (PCA) в Python

Ответ 1

Вы можете найти функцию PCA в модуле matplotlib:

import numpy as np
from matplotlib.mlab import PCA

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

будут сохранены различные параметры СПС. Это из mlab-части matplotlib, которая является слоем совместимости с синтаксисом MATLAB

EDIT: в блоге nextgenetics Я нашел замечательную демонстрацию того, как выполнять и отображать PCA с модулем matplotlib mlab, получать удовольствие и проверять этот блог

Ответ 2

Я отправил свой ответ, хотя еще один ответ уже принят; принятый ответ основан на устаревшей функции; кроме того, эта устаревшая функция основана на сингулярном декомпозиции значений (SVD), который (хотя и совершенно корректный) является гораздо более интенсивным с точки зрения памяти и процессора двух общих методов расчета PCA. Это особенно актуально здесь из-за размера массива данных в OP. Используя ковариационный PCA, массив, используемый в потоке вычислений, составляет всего 144 x 144, а не 26424 x 144 (размеры исходного массива данных).

Вот простая рабочая реализация PCA с использованием модуля linalg от SciPy. Поскольку эта реализация сначала вычисляет матрицу ковариации, а затем выполняет все последующие вычисления в этом массиве, она использует гораздо меньше памяти, чем SVA на основе SVD.

(модуль linalg в NumPy также может использоваться без изменения кода ниже, кроме инструкции import, которая будет из numpy import linalg как LA.)

Два ключевых шага в этой реализации PCA:

• вычисление ковариационной матрицы; и

• взяв eivenvectors и собственные значения этой матрицы cov

В приведенной ниже функции параметр dims_rescaled_data​​strong > относится к требуемому количеству измерений в масштабируемой матрице данных; этот параметр имеет значение по умолчанию только двух измерений, но приведенный ниже код не ограничен двумя, но может быть любым значением, меньшим, чем номер столбца исходного массива данных.

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-преобразование чисто разделяет класс я из класса II и класса III (но не класса II из класса III, который на самом деле требует другого измерения).

Ответ 3

Другой Python PCA с использованием numpy. Та же идея, что и @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(data):
    """
    Covariance matrix
    note: specifically for mean-centered data
    note: numpy `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)

Как я понимаю, использование собственных значений (первый способ) лучше для высокоразмерных данных и меньшего количества выборок, тогда как использование сингулярного разложения значений лучше, если у вас больше образцов, чем размеры.

Ответ 4

Это задание для numpy.

И вот учебник, демонстрирующий, как анализ pincipal компонентов может быть выполнен с использованием 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

Ответ 5

Вот варианты scikit-learn. При обоих методах использовался StandardScaler, поскольку PCA зависит от масштаба

Метод 1: Попросите scikit-learn выбрать минимальное количество основных компонентов, чтобы сохранить не менее x% (90% в приведенном ниже примере) дисперсии.

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

Ответ 6

UPDATE: matplotlib.mlab.PCA с момента выпуска 2.2 (2018-03-06) действительно устарел.

Библиотека matplotlib.mlab.PCA (используется в этом ответе) не устарела. Поэтому для всех людей, прибывающих сюда через Google, я отправлю полный рабочий пример, протестированный с помощью Python 2.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.pyplot чтобы нарисовать эти данные, чтобы убедиться, что PCA дает "хорошие" результаты. Список names используется только для аннотирования наших пяти векторов.

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

Посмотрев на наши исходные векторы, мы увидим, что данные [0] ("A") и данные [3] ("D") довольно схожи, как и данные [1] ("B") и данные [2] (" C "). Это отражено в 2D-графике данных, преобразованных в PCA.

Ответ 7

В дополнение ко всем остальным ответам, здесь приведен код для построения графика biplot с использованием 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()

Ответ 8

Я сделал небольшой скрипт для сравнения различных PCAs, представленных в качестве ответа здесь:

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 - это ответ в ответе Марка.

pca_doug - это ответ на вопрос о даге.

Вот пример вывода (но результат во многом зависит от размера данных и параметров 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

РЕДАКТИРОВАТЬ:

Функция diffsnorm из fbpca вычисляет ошибку спектрально-нормой разложения Шура.

Ответ 9

Ради 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)

Ответ 10

этот пример кода загружает японскую кривую доходности и создает компоненты 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()

Ответ 11

введите описание изображения здесь

Как рассчитать PCA для скорости ветра (рис. В качестве образца) с данными этого типа в формате .csv.