note_2
Устранение избыточности в данных.
Шокуров Антон В.
shokurov.anton.v@yandex.ru
http://машинноезрение.рф
Версия 0.13

Анотация

В данной заметке подйет речь о преобразования над самими признаками, т.е. над матрицей X. Значения каждого из прецендента нам не важно. Постараемся из X выжать как можно больше данных.

Это предварительная версия! Любые замечания приветсвуются.

In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Метод Главных Компонент

Синтезированные примеры

Формирование данных на приямой, которая под углом.

In [2]:
# Сгенерируем данные
x = np.random.uniform(-10, 10, 30) # Раномерная выборка на отрезке [-10, 10] в количестве 30 элементов.
x.shape
Out[2]:
(30,)
In [3]:
X = np.row_stack( (x, np.zeros_like(x) ) ) # Две строчки -- к каждому числу добавлена вторая координата (равная 0).
X # Выводим матрицу.
Out[3]:
array([[ 3.52083406,  2.13637234,  5.39226446, -3.42488832,  6.98779202,
         2.72217559, -1.07567581, -1.97679152,  5.73396276, -5.68561776,
         4.87352042, -7.35158786,  5.8722167 ,  1.43974041, -6.58981135,
        -4.02080892,  0.22950724, -5.36340312,  8.22779209, -4.0727585 ,
        -4.5386593 , -7.23263026, -6.96286994, -2.59496404, -8.99609292,
        -8.59966218, -2.40104642,  4.0441422 ,  9.87070958, -0.12108042],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ,  0.        ,  0.        ,  0.        ,  0.        ]])
In [4]:
plt.plot(X[0,:], X[1,:], ".") # Нарисуем их
# Они находятся на прямой.
Out[4]:
[<matplotlib.lines.Line2D at 0x7faa219a7be0>]
In [5]:
# Теперь подкрутим данные. Разместим их на под уголом.
theta = np.radians( 30 ) # Градусы в радианы.
c, s = np.cos(theta), np.sin(theta) # Поэлементвное присвоение,
# Фактически эквивалентно:
# c = np.cos(theta)
# s = np.sin(theta)
In [6]:
R = np.array([[c,-s], [s, c]]) # Формируем матрицу поворота,
T = np.array([[2.5], [-1.3]]) # а также вектор смещения.
In [7]:
XX = R.dot(X) + T # Применяем к нашим данным вращение и трансляцию.
XXt = XX.transpose()
XXt
Out[7]:
array([[ 5.54913173,  0.46041703],
       [ 4.35015272, -0.23181383],
       [ 7.16983801,  1.39613223],
       [-0.46604029, -3.01244416],
       [ 8.55160541,  2.19389601],
       [ 4.85747322,  0.0610878 ],
       [ 1.56843742, -1.8378379 ],
       [ 0.78804832, -2.28839576],
       [ 7.46575742,  1.56698138],
       [-2.42388941, -4.14280888],
       [ 6.72059249,  1.13676021],
       [-3.86666185, -4.97579393],
       [ 7.58548884,  1.63610835],
       [ 3.74685177, -0.5801298 ],
       [-3.20694403, -4.59490567],
       [-0.98212267, -3.31040446],
       [ 2.6987591 , -1.18524638],
       [-2.14484335, -3.98170156],
       [ 9.62547697,  2.81389605],
       [-1.02711232, -3.33637925],
       [-1.43059425, -3.56932965],
       [-3.76364154, -4.91631513],
       [-3.53002225, -4.78143497],
       [ 0.25269522, -2.59748202],
       [-5.290845  , -5.79804646],
       [-4.94752591, -5.59983109],
       [ 0.4206328 , -2.50052321],
       [ 6.00232988,  0.7220711 ],
       [11.04828525,  3.63535479],
       [ 2.39514128, -1.36054021]])
In [8]:
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, ".")
# Данные теперь под углом. Ордината (красная точка) находится сбоку.
Out[8]:
[<matplotlib.lines.Line2D at 0x7faa218c4160>,
 <matplotlib.lines.Line2D at 0x7faa218c4358>]

Метод главных компонент (PCA), данные на прямой.

Как восстановить их? Как устранить факт того, что они под уголом? Для этого есть метод главных компонент.

Он позволяет снизить размерность данных, т.е. найти подпространство, которое содержит данные. Самое пространство и преобразование линейны.

In [9]:
# Загружаем метод главных компонент (PCA -- Principal component analysis):
from sklearn.decomposition import PCA
In [10]:
# PCA это конструктор объектов. Для выполнения вычислений нужно объект создать.
# Самое главное задать требуемое размер пространства.
pca = PCA( n_components = 1 ) # Хотим преобзовать данные к прямой, т.е. сделать их одномерными (как были исходные).
In [11]:
pca.fit( XXt ) # Запускаем метод pca на данных XXt, т.е. ишим подплоскость (смещение и поворот) содержащую данные.
# Выводятся используемые значения параметров метода.
Out[11]:
PCA(copy=True, iterated_power='auto', n_components=1, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [12]:
XXtd = pca.transform( XXt ) # Применяем найденное (сохраненное в объекте pca) преобразование к данным.
XXtd
Out[12]:
array([[ 4.18607801],
       [ 2.8016163 ],
       [ 6.05750842],
       [-2.75964436],
       [ 7.65303598],
       [ 3.38741955],
       [-0.41043185],
       [-1.31154756],
       [ 6.39920672],
       [-5.0203738 ],
       [ 5.53876437],
       [-6.6863439 ],
       [ 6.53746066],
       [ 2.10498437],
       [-5.92456739],
       [-3.35556496],
       [ 0.8947512 ],
       [-4.69815916],
       [ 8.89303605],
       [-3.40751454],
       [-3.87341534],
       [-6.5673863 ],
       [-6.29762598],
       [-1.92972009],
       [-8.33084896],
       [-7.93441822],
       [-1.73580246],
       [ 4.70938616],
       [10.53595354],
       [ 0.54416353]])
In [13]:
pca.explained_variance_ratio_ # Хорошо-ли данные коррелируют с вычисленным преобзованием?
# Фактически говорит о том насколько хорошо удалось снизить размерность,
# много ли данных оказалось искаженными (спроецированными с сильным смещением)
# 1 -- 100% по данных по первой оси (она же единственная).
Out[13]:
array([1.])
In [14]:
# Теперь нужно сравнить восстановленные данные с исходными. Для этого их нужно привести к нормальному виду:
from sklearn import preprocessing # preprocessing -- первичные преобразования данных.
# Будем применять preprocessing.scale -- делает выборочное среднее равно нулю, а дисперсию равной 1.
# Иногда, для корретной работы примера, требуется:
XXtd = -XXtd # см ниже.
XXtds = preprocessing.scale( XXtd[:, 0] ) # Нормируем. Используем функию scale из модуля sklearn.preprocessing.
XXtds
Out[14]:
array([-0.78155872, -0.52307378, -1.1309628 ,  0.51523744, -1.42885464,
       -0.63244576,  0.07662939,  0.24487155, -1.19475934,  0.93732532,
       -1.03411106,  1.24836908, -1.22057194, -0.39300961,  1.10614214,
        0.62649837, -0.16705389,  0.87716646, -1.66036797,  0.63619758,
        0.72318326,  1.22615919,  1.17579378,  0.36028702,  1.55540523,
        1.48138991,  0.32408176, -0.87926259, -1.96710771, -0.10159767])
In [15]:
xs = preprocessing.scale( x ) # x -- это исходные даные. Их тоже нормируем.
xs
Out[15]:
array([ 0.78155872,  0.52307378,  1.1309628 , -0.51523744,  1.42885464,
        0.63244576, -0.07662939, -0.24487155,  1.19475934, -0.93732532,
        1.03411106, -1.24836908,  1.22057194,  0.39300961, -1.10614214,
       -0.62649837,  0.16705389, -0.87716646,  1.66036797, -0.63619758,
       -0.72318326, -1.22615919, -1.17579378, -0.36028702, -1.55540523,
       -1.48138991, -0.32408176,  0.87926259,  1.96710771,  0.10159767])
In [16]:
# Вычислим расстояние между векторами, т.е. отнормированными данными.
np.linalg.norm((xs - XXtds)) # Вызываем функцию из модуля, без from...
# Длина фактически равна нулю. Значит данные были восстановлены полностью.
# Если норма не равна нулю, то требуется массив отразить. см выше: XXtd = -XXtd. Нужно убрать комментарий.
Out[16]:
10.954451150103322
In [17]:
# Что будет если на две оси спроецировать?
pca = PCA( n_components = 2 ) # Настраиваем объект на двумерную плоскость.
XXtd = pca.fit_transform( XXt ) # Настроим параметры модели и сразу преобразуем данные к двум компонентам.
XXtd # Видим, что даные состоят из двумерных точек. Вторая компонента фактически равна нулю (машинное эпсилон).
Out[17]:
array([[ 4.18607801e+00,  1.64309362e-15],
       [ 2.80161630e+00, -2.45545890e-16],
       [ 6.05750842e+00, -1.54576389e-16],
       [-2.75964436e+00, -2.08490654e-16],
       [ 7.65303598e+00,  2.87782376e-16],
       [ 3.38741955e+00, -1.29735034e-16],
       [-4.10431851e-01,  1.35620592e-16],
       [-1.31154756e+00,  1.05830470e-16],
       [ 6.39920672e+00,  1.58513331e-16],
       [-5.02037380e+00, -1.48894016e-16],
       [ 5.53876437e+00, -9.89519877e-17],
       [-6.68634390e+00, -1.93668576e-16],
       [ 6.53746066e+00, -3.01124555e-16],
       [ 2.10498437e+00, -1.54045925e-16],
       [-5.92456739e+00, -1.38602972e-16],
       [-3.35556496e+00,  1.66139144e-16],
       [ 8.94751200e-01,  9.19300560e-17],
       [-4.69815916e+00, -2.30279004e-16],
       [ 8.89303605e+00,  4.35889764e-16],
       [-3.40751454e+00, -1.67018406e-16],
       [-3.87341534e+00,  3.87034397e-16],
       [-6.56738630e+00, -3.21898668e-16],
       [-6.29762598e+00, -3.08676439e-16],
       [-1.92972009e+00,  1.56342470e-18],
       [-8.33084896e+00,  1.13003583e-15],
       [-7.93441822e+00,  3.80281705e-16],
       [-1.73580246e+00,  1.07216379e-16],
       [ 4.70938616e+00, -1.53763245e-16],
       [ 1.05359535e+01, -2.52768153e-16],
       [ 5.44163535e-01, -6.94761054e-17]])
In [18]:
# Построим преобращованные (с целью снижения размерность) точки.
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график будет обманчевым. Ось по y отмасштабирована (разные масштабы)!
Out[18]:
[<matplotlib.lines.Line2D at 0x7faa15daabe0>]
In [19]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" )
plt.axis('equal') # Деламем оси одинакового масштаба.
# Правильнее так. Теперь видно, что данные расположены вдоль оси x.
Out[19]:
(-9.274189085506404,
 11.47929366739617,
 -4.2014828272745842e-16,
 1.7413432304822776e-15)
In [20]:
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось все содержит. Для второй оси коеффициент равен 0.
Out[20]:
array([1.00000000e+00, 6.23501456e-33])

Видно, что увеличиение количества осей мы не теряем информацию о том, что данные на самом деле одномерные.

Метод главных компонент (PCA), данные на прямой плюс шум

In [21]:
# Добавим шума (к y, раньше он был раен 0)
X[1,:] = np.random.randn( (30) )/4 # вторая компонента содержит номальный шум, раньше там был 0.
XX = R.dot(X) + T # Как и раньше применяем вращение и трансляцию.
XXt = XX.transpose() # Для удобства. Например, при вызове plot, fit_transform.
In [22]:
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, ".")
Out[22]:
[<matplotlib.lines.Line2D at 0x7faa15cfb160>,
 <matplotlib.lines.Line2D at 0x7faa15cfb2e8>]
In [23]:
pca = PCA( n_components = 1 ) # Настраиваем объект на одну компоненту.
XXtd = pca.fit_transform( XXt ) # Ищим параметры преобразования и преобразуем.
In [24]:
pca.explained_variance_ratio_ # Теперь важность первой оси немного упала. Но совсем немного!
Out[24]:
array([0.9984875])
In [25]:
# Теперь спроецируем на две оси.
pca = PCA( n_components = 2 )
XXtd = pca.fit_transform( XXt )
XXtd # Видно, что вторая компонента фактически равна шуму.
Out[25]:
array([[ 4.18713846e+00, -1.22973516e-01],
       [ 2.80169926e+00, -1.78193069e-04],
       [ 6.05496916e+00,  3.55826164e-01],
       [-2.76074351e+00,  1.33487484e-01],
       [ 7.65625210e+00, -3.92194437e-01],
       [ 3.38956849e+00, -2.68642538e-01],
       [-4.11067365e-01,  8.17032074e-02],
       [-1.31355222e+00,  2.57659201e-01],
       [ 6.39999356e+00, -7.86781996e-02],
       [-5.02103670e+00,  6.77003980e-02],
       [ 5.53672679e+00,  2.88113629e-01],
       [-6.68733878e+00,  1.04840243e-01],
       [ 6.53573445e+00,  2.51125714e-01],
       [ 2.10658760e+00, -2.02034683e-01],
       [-5.92539410e+00,  8.57126464e-02],
       [-3.35646125e+00,  1.04632631e-01],
       [ 8.93249944e-01,  2.00119130e-01],
       [-4.69665830e+00, -2.14580675e-01],
       [ 8.89463522e+00, -1.75598150e-01],
       [-3.40508985e+00, -3.30701936e-01],
       [-3.87339223e+00, -1.78093860e-02],
       [-6.56658829e+00, -1.29621698e-01],
       [-6.29689449e+00, -1.19876616e-01],
       [-1.93379661e+00,  5.26769994e-01],
       [-8.33124043e+00,  1.95024213e-02],
       [-7.93316554e+00, -1.94412514e-01],
       [-1.73464421e+00, -1.58386059e-01],
       [ 4.71052695e+00, -1.31503717e-01],
       [ 1.05361590e+01,  1.32779500e-02],
       [ 5.43822803e-01,  4.67215044e-02]])
In [26]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график как и раньше будет обманный. Ось по y отмасштабирована!
Out[26]:
[<matplotlib.lines.Line2D at 0x7faa15d20518>]
In [27]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Правильнее так. Теперь видно, что данные расположены вдоль оси x.
plt.axis('equal') # Деламем оси одинакового масштаба. Вдоль y расположен шум.
Out[27]:
(-9.274610400447406,
 11.47952902202962,
 -0.4381426585664787,
 0.5727182154961763)
In [28]:
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось содержит почти все. Для второй оси коеффициент почти равен 0, но не 0.
Out[28]:
array([0.9984875, 0.0015125])
In [29]:
XXtd[:,1] # Это шум.
Out[29]:
array([-1.22973516e-01, -1.78193069e-04,  3.55826164e-01,  1.33487484e-01,
       -3.92194437e-01, -2.68642538e-01,  8.17032074e-02,  2.57659201e-01,
       -7.86781996e-02,  6.77003980e-02,  2.88113629e-01,  1.04840243e-01,
        2.51125714e-01, -2.02034683e-01,  8.57126464e-02,  1.04632631e-01,
        2.00119130e-01, -2.14580675e-01, -1.75598150e-01, -3.30701936e-01,
       -1.78093860e-02, -1.29621698e-01, -1.19876616e-01,  5.26769994e-01,
        1.95024213e-02, -1.94412514e-01, -1.58386059e-01, -1.31503717e-01,
        1.32779500e-02,  4.67215044e-02])
In [30]:
import scipy.stats as models
models.norm.fit( XXtd[:,1] )[1], models.norm.fit( X[1,:] )[1]
# Почти как исходное 1/4 = 0.25. У исходного и преобразованного оценка почти совпадает.
Out[30]:
(0.2084655146248336, 0.21242956176112457)

Аналогичный пример в 3d

данные находятся в плоскости.

In [31]:
from sklearn.decomposition import PCA

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
In [32]:
e = np.exp(1)
np.random.seed(4)


def pdf(x):
    return 0.5 * (stats.norm(scale=0.25 / e).pdf(x)
                  + stats.norm(scale=4 / e).pdf(x))

y = np.random.normal(scale=0.5, size=(30000))
x = np.random.normal(scale=0.5, size=(30000))
z = np.random.normal(scale=0.1, size=len(x))

density = pdf(x) * pdf(y)
pdf_z = pdf(5 * z)

density *= pdf_z

a = x + y
b = 2 * y
c = a - b + z

norm = np.sqrt(a.var() + b.var())
a /= norm
b /= norm
In [33]:
def plot_figs(fig_num, elev, azim):
    fig = plt.figure(fig_num, figsize=(8, 6))
    plt.clf()
    ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)

    ax.scatter(a[::10], b[::10], c[::10], c=density[::10], marker='+', alpha=.4)
    Y = np.c_[a, b, c]

    # Using SciPy's SVD, this would be:
    # _, pca_score, V = scipy.linalg.svd(Y, full_matrices=False)

    pca = PCA(n_components=3)
    pca.fit(Y)
    pca_score = pca.explained_variance_ratio_
    V = pca.components_

    x_pca_axis, y_pca_axis, z_pca_axis = V.T * pca_score / pca_score.min()

    x_pca_axis, y_pca_axis, z_pca_axis = 3 * V.T
    x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]
    y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]
    z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]
    x_pca_plane.shape = (2, 2)
    y_pca_plane.shape = (2, 2)
    z_pca_plane.shape = (2, 2)
    ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)
    ax.w_xaxis.set_ticklabels([])
    ax.w_yaxis.set_ticklabels([])
    ax.w_zaxis.set_ticklabels([])
In [34]:
elev = -40
azim = -80
plot_figs(1, elev, azim)
In [35]:
elev = 30
azim = 20
plot_figs(2, elev, azim)

Другие примеры с PCA

In [36]:
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
X
Out[36]:
array([[-1, -1],
       [-2, -1],
       [-3, -2],
       [ 1,  1],
       [ 2,  1],
       [ 3,  2]])
In [37]:
plt.plot(X[:,0], X[:,1], "*")
plt.axis('equal')
Out[37]:
(-3.3, 3.3, -2.2, 2.2)
In [38]:
pca = PCA(n_components=2)
In [39]:
pca.fit(X)
Out[39]:
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='auto', tol=0.0, whiten=False)
In [40]:
print(pca.explained_variance_ratio_) 
[0.99244289 0.00755711]
In [41]:
print(pca.singular_values_)
[6.30061232 0.54980396]
In [42]:
X2 = pca.transform(X)
X2
Out[42]:
array([[ 1.38340578,  0.2935787 ],
       [ 2.22189802, -0.25133484],
       [ 3.6053038 ,  0.04224385],
       [-1.38340578, -0.2935787 ],
       [-2.22189802,  0.25133484],
       [-3.6053038 , -0.04224385]])
In [43]:
plt.plot(X2[:,0], X2[:,1], "*")
#plt.axis('equal')
Out[43]:
[<matplotlib.lines.Line2D at 0x7faa15792d30>]
In [44]:
pca_full = PCA(n_components=2, svd_solver='full')
pca_full.fit( X )
Out[44]:
PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
  svd_solver='full', tol=0.0, whiten=False)
In [45]:
pca_arpack = PCA(n_components=1, svd_solver='arpack') #долюно быть n_components < n_features
pca_arpack.fit( X )
Out[45]:
PCA(copy=True, iterated_power='auto', n_components=1, random_state=None,
  svd_solver='arpack', tol=0.0, whiten=False)
In [46]:
X_arpack=pca_arpack.transform(X)
X_arpack
Out[46]:
array([[ 1.38340578],
       [ 2.22189802],
       [ 3.6053038 ],
       [-1.38340578],
       [-2.22189802],
       [-3.6053038 ]])
In [47]:
plt.plot(np.sort(X_arpack[:,0]), "*")
Out[47]:
[<matplotlib.lines.Line2D at 0x7faa1577b898>]

Ядровый Метод Главных Компонент

Kernel PCA
Выделяемые компоненты совсем не обязательно должны быть "прямыми". Они могут быть и кривыми засчет подмены понятия скалярного произведения. Оно заменяется на так называемое ядро.

In [48]:
from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_circles
In [49]:
# Для обеспечения повтаряемости эксперимента
np.random.seed(0)
In [50]:
# Создадим "кривые" данные, две концентрические окружности с центром в точке (0, 0)
X, y = make_circles(n_samples=400, factor=.3, noise=.05)
X.shape, y.shape
Out[50]:
((400, 2), (400,))
In [51]:
np.unique(y)
Out[51]:
array([0, 1])
In [52]:
plt.title("Исходное пространство")
reds = y == 0
blues = y == 1

plt.scatter(X[reds, 0], X[reds, 1], c="red",
            s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
            s=20, edgecolor='k')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[52]:
Text(0,0.5,'$x_2$')
In [53]:
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
In [54]:
plt.scatter(X_pca[reds, 0], X_pca[reds, 1], c="red",
            s=20, edgecolor='k')
plt.scatter(X_pca[blues, 0], X_pca[blues, 1], c="blue",
            s=20, edgecolor='k')
plt.title("После МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
Out[54]:
Text(0,0.5,'2ая главная компонента')

В общем ничего не поменялось. Теперь пробуем ядровый метод главных компонент.

In [55]:
kpca = KernelPCA(n_components=2, kernel="rbf", fit_inverse_transform=True, gamma=10) # Ядро rbf, radial basis function
X_kpca = kpca.fit_transform(X)
X_back = kpca.inverse_transform(X_kpca)
In [56]:
plt.scatter(X_kpca[reds, 0], X_kpca[reds, 1], c="red",
            s=20, edgecolor='orange')
plt.scatter(X_kpca[blues, 0], X_kpca[blues, 1], c="blue",
            s=20, edgecolor='c')
plt.title("После Ядровый МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
Out[56]:
Text(0,0.5,'2ая главная компонента')

Исходные два множества четко разделились.

Начертим кривые координаты.

In [57]:
l = np.linspace(-1.5, 1.5, 5)
l
Out[57]:
array([-1.5 , -0.75,  0.  ,  0.75,  1.5 ])
In [58]:
grid = np.meshgrid( l, l )
grid
Out[58]:
[array([[-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
        [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
        [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
        [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
        [-1.5 , -0.75,  0.  ,  0.75,  1.5 ]]),
 array([[-1.5 , -1.5 , -1.5 , -1.5 , -1.5 ],
        [-0.75, -0.75, -0.75, -0.75, -0.75],
        [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
        [ 0.75,  0.75,  0.75,  0.75,  0.75],
        [ 1.5 ,  1.5 ,  1.5 ,  1.5 ,  1.5 ]])]
In [59]:
grid[0][1,3],grid[1][1,3]
Out[59]:
(0.75, -0.75)
In [60]:
type(grid), len(grid) #! Это список, а не массив !
Out[60]:
(list, 2)
In [61]:
grid[0]
Out[61]:
array([[-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
       [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
       [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
       [-1.5 , -0.75,  0.  ,  0.75,  1.5 ],
       [-1.5 , -0.75,  0.  ,  0.75,  1.5 ]])
In [62]:
grid[1]
Out[62]:
array([[-1.5 , -1.5 , -1.5 , -1.5 , -1.5 ],
       [-0.75, -0.75, -0.75, -0.75, -0.75],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.75,  0.75,  0.75,  0.75,  0.75],
       [ 1.5 ,  1.5 ,  1.5 ,  1.5 ,  1.5 ]])
In [63]:
np.ravel([[1, 2],[3,4]]) # выпрямляет массив в одномерный, схоже с flatten
Out[63]:
array([1, 2, 3, 4])
In [64]:
plt.title("Кривые координаты в исходном пространстве")
reds = y == 0
blues = y == 1

plt.scatter(X[reds, 0], X[reds, 1], c="red",
            s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
            s=20, edgecolor='k')

X1, X2 = np.meshgrid(np.linspace(-1.0, 1.0, 50), np.linspace(-1.0, 1.0, 50))
X_grid = np.array([np.ravel(X1), np.ravel(X2)]).T
# projection on the first principal component (in the phi space)
Z_grid = kpca.transform(X_grid)[:, 0].reshape(X1.shape)
plt.contour(X1, X2, Z_grid, colors='grey', linewidths=1, origin='lower')

plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[64]:
Text(0,0.5,'$x_2$')
In [65]:
plt.scatter(X_back[reds, 0], X_back[reds, 1], c="red",
            s=20, edgecolor='k')
plt.scatter(X_back[blues, 0], X_back[blues, 1], c="blue",
            s=20, edgecolor='k')
plt.title("Обратно в исходном пространстве")
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[65]:
Text(0,0.5,'$x_2$')

Метод главных компонент на примере Iris dataset

In [66]:
# Загрузка данных из самой бибилотеки. 
from sklearn import datasets # В dataset есть нексолько традиционых наборов.
In [67]:
iris = datasets.load_iris()

X = iris.data # Сами данные
y = iris.target # Класс, в данной заметке он нужен нам только для визуализации.
target_names = iris.target_names # Имена классов.
X.shape, y.shape, target_names.shape # Покажем чему раны размеры.
Out[67]:
((150, 4), (150,), (3,))

Инача говоря данные содержат 150 примеров. Для каждого примера заданы 4 признака, а также разметку (чем пример является). Для числовых обозначений классов даны соответствующие текстовые названия.

In [68]:
X # Данные, четыре признака. 150 примера (прецендента).
Out[68]:
array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       [4.6, 3.1, 1.5, 0.2],
       [5. , 3.6, 1.4, 0.2],
       [5.4, 3.9, 1.7, 0.4],
       [4.6, 3.4, 1.4, 0.3],
       [5. , 3.4, 1.5, 0.2],
       [4.4, 2.9, 1.4, 0.2],
       [4.9, 3.1, 1.5, 0.1],
       [5.4, 3.7, 1.5, 0.2],
       [4.8, 3.4, 1.6, 0.2],
       [4.8, 3. , 1.4, 0.1],
       [4.3, 3. , 1.1, 0.1],
       [5.8, 4. , 1.2, 0.2],
       [5.7, 4.4, 1.5, 0.4],
       [5.4, 3.9, 1.3, 0.4],
       [5.1, 3.5, 1.4, 0.3],
       [5.7, 3.8, 1.7, 0.3],
       [5.1, 3.8, 1.5, 0.3],
       [5.4, 3.4, 1.7, 0.2],
       [5.1, 3.7, 1.5, 0.4],
       [4.6, 3.6, 1. , 0.2],
       [5.1, 3.3, 1.7, 0.5],
       [4.8, 3.4, 1.9, 0.2],
       [5. , 3. , 1.6, 0.2],
       [5. , 3.4, 1.6, 0.4],
       [5.2, 3.5, 1.5, 0.2],
       [5.2, 3.4, 1.4, 0.2],
       [4.7, 3.2, 1.6, 0.2],
       [4.8, 3.1, 1.6, 0.2],
       [5.4, 3.4, 1.5, 0.4],
       [5.2, 4.1, 1.5, 0.1],
       [5.5, 4.2, 1.4, 0.2],
       [4.9, 3.1, 1.5, 0.1],
       [5. , 3.2, 1.2, 0.2],
       [5.5, 3.5, 1.3, 0.2],
       [4.9, 3.1, 1.5, 0.1],
       [4.4, 3. , 1.3, 0.2],
       [5.1, 3.4, 1.5, 0.2],
       [5. , 3.5, 1.3, 0.3],
       [4.5, 2.3, 1.3, 0.3],
       [4.4, 3.2, 1.3, 0.2],
       [5. , 3.5, 1.6, 0.6],
       [5.1, 3.8, 1.9, 0.4],
       [4.8, 3. , 1.4, 0.3],
       [5.1, 3.8, 1.6, 0.2],
       [4.6, 3.2, 1.4, 0.2],
       [5.3, 3.7, 1.5, 0.2],
       [5. , 3.3, 1.4, 0.2],
       [7. , 3.2, 4.7, 1.4],
       [6.4, 3.2, 4.5, 1.5],
       [6.9, 3.1, 4.9, 1.5],
       [5.5, 2.3, 4. , 1.3],
       [6.5, 2.8, 4.6, 1.5],
       [5.7, 2.8, 4.5, 1.3],
       [6.3, 3.3, 4.7, 1.6],
       [4.9, 2.4, 3.3, 1. ],
       [6.6, 2.9, 4.6, 1.3],
       [5.2, 2.7, 3.9, 1.4],
       [5. , 2. , 3.5, 1. ],
       [5.9, 3. , 4.2, 1.5],
       [6. , 2.2, 4. , 1. ],
       [6.1, 2.9, 4.7, 1.4],
       [5.6, 2.9, 3.6, 1.3],
       [6.7, 3.1, 4.4, 1.4],
       [5.6, 3. , 4.5, 1.5],
       [5.8, 2.7, 4.1, 1. ],
       [6.2, 2.2, 4.5, 1.5],
       [5.6, 2.5, 3.9, 1.1],
       [5.9, 3.2, 4.8, 1.8],
       [6.1, 2.8, 4. , 1.3],
       [6.3, 2.5, 4.9, 1.5],
       [6.1, 2.8, 4.7, 1.2],
       [6.4, 2.9, 4.3, 1.3],
       [6.6, 3. , 4.4, 1.4],
       [6.8, 2.8, 4.8, 1.4],
       [6.7, 3. , 5. , 1.7],
       [6. , 2.9, 4.5, 1.5],
       [5.7, 2.6, 3.5, 1. ],
       [5.5, 2.4, 3.8, 1.1],
       [5.5, 2.4, 3.7, 1. ],
       [5.8, 2.7, 3.9, 1.2],
       [6. , 2.7, 5.1, 1.6],
       [5.4, 3. , 4.5, 1.5],
       [6. , 3.4, 4.5, 1.6],
       [6.7, 3.1, 4.7, 1.5],
       [6.3, 2.3, 4.4, 1.3],
       [5.6, 3. , 4.1, 1.3],
       [5.5, 2.5, 4. , 1.3],
       [5.5, 2.6, 4.4, 1.2],
       [6.1, 3. , 4.6, 1.4],
       [5.8, 2.6, 4. , 1.2],
       [5. , 2.3, 3.3, 1. ],
       [5.6, 2.7, 4.2, 1.3],
       [5.7, 3. , 4.2, 1.2],
       [5.7, 2.9, 4.2, 1.3],
       [6.2, 2.9, 4.3, 1.3],
       [5.1, 2.5, 3. , 1.1],
       [5.7, 2.8, 4.1, 1.3],
       [6.3, 3.3, 6. , 2.5],
       [5.8, 2.7, 5.1, 1.9],
       [7.1, 3. , 5.9, 2.1],
       [6.3, 2.9, 5.6, 1.8],
       [6.5, 3. , 5.8, 2.2],
       [7.6, 3. , 6.6, 2.1],
       [4.9, 2.5, 4.5, 1.7],
       [7.3, 2.9, 6.3, 1.8],
       [6.7, 2.5, 5.8, 1.8],
       [7.2, 3.6, 6.1, 2.5],
       [6.5, 3.2, 5.1, 2. ],
       [6.4, 2.7, 5.3, 1.9],
       [6.8, 3. , 5.5, 2.1],
       [5.7, 2.5, 5. , 2. ],
       [5.8, 2.8, 5.1, 2.4],
       [6.4, 3.2, 5.3, 2.3],
       [6.5, 3. , 5.5, 1.8],
       [7.7, 3.8, 6.7, 2.2],
       [7.7, 2.6, 6.9, 2.3],
       [6. , 2.2, 5. , 1.5],
       [6.9, 3.2, 5.7, 2.3],
       [5.6, 2.8, 4.9, 2. ],
       [7.7, 2.8, 6.7, 2. ],
       [6.3, 2.7, 4.9, 1.8],
       [6.7, 3.3, 5.7, 2.1],
       [7.2, 3.2, 6. , 1.8],
       [6.2, 2.8, 4.8, 1.8],
       [6.1, 3. , 4.9, 1.8],
       [6.4, 2.8, 5.6, 2.1],
       [7.2, 3. , 5.8, 1.6],
       [7.4, 2.8, 6.1, 1.9],
       [7.9, 3.8, 6.4, 2. ],
       [6.4, 2.8, 5.6, 2.2],
       [6.3, 2.8, 5.1, 1.5],
       [6.1, 2.6, 5.6, 1.4],
       [7.7, 3. , 6.1, 2.3],
       [6.3, 3.4, 5.6, 2.4],
       [6.4, 3.1, 5.5, 1.8],
       [6. , 3. , 4.8, 1.8],
       [6.9, 3.1, 5.4, 2.1],
       [6.7, 3.1, 5.6, 2.4],
       [6.9, 3.1, 5.1, 2.3],
       [5.8, 2.7, 5.1, 1.9],
       [6.8, 3.2, 5.9, 2.3],
       [6.7, 3.3, 5.7, 2.5],
       [6.7, 3. , 5.2, 2.3],
       [6.3, 2.5, 5. , 1.9],
       [6.5, 3. , 5.2, 2. ],
       [6.2, 3.4, 5.4, 2.3],
       [5.9, 3. , 5.1, 1.8]])
In [69]:
y # Класса. Фактически сначала идут преценденты относящиеся к первому классу, потом второму, а уже потом третьему.
Out[69]:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])
In [70]:
target_names # Имена классов, т.е. 0 классу соответствует имя setosa, ..., 2 -- virginica
# Да, тип элементов может быть текстом!
Out[70]:
array(['setosa', 'versicolor', 'virginica'], dtype='<U10')
In [71]:
np.unique(y) # Найдем номера различающихся классов. Функция ищет уникальные элементы массива.
Out[71]:
array([0, 1, 2])
In [72]:
# Напомню.
y * 2 + 1 # Выполняет данное вычисление поэлементно. 
Out[72]:
array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
       3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
       3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
       5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
       5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5])
In [73]:
# По аналогии.
y==1 # Применили проверку на равенство 1 для каждого из элементов массива. Истина -- True, ложь -- False.
Out[73]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True,  True,  True,  True,  True,  True,  True,  True,  True,
        True, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False])
In [74]:
# Напечатам размер каждого из классов, т.е. количество прецендентов в каждом классе.
# Далее идет так называемый цикл, конкретно цикл for. Он выполняет некий код для каждого из элементов по одному.
for i in np.unique(y): # i пробегает список элементов [0,1,2], т.е. принимает каждое из данных значений по отдельности.
    # Тело начинается с отступа. Замечу, в предыдущей строчке в конце был символ :!
    print( "размер класса ", i, sum( y==i ) ) # Для каждого из жлементов выполянется данный код. Идет печать.
    # В арифметических операциях истина (True) равна 1, а лож (False).
# Видим, что классы имеют один и тот же размер, равный 50.
размер класса  0 50
размер класса  1 50
размер класса  2 50
In [75]:
nums = np.unique(y) # Фактически это [0, 1, 2]
list( zip( nums, target_names ) ) # Списки сшиваются по-элементно вместе. Как молния с бегунком у куртки.
# nums[0] и target_names[0] в первый элемент нового списка... (0, 'setosa') ... и так далее
# по аналогии все далается при большем количестве списоков. Рзамеры их должны сопадать!
Out[75]:
[(0, 'setosa'), (1, 'versicolor'), (2, 'virginica')]
In [76]:
X_1 = X[ y == 1 ] # Выбираем только те элементы для которых условие истинно, т.е. относящиеся к классу 1.
X_1
Out[76]:
array([[7. , 3.2, 4.7, 1.4],
       [6.4, 3.2, 4.5, 1.5],
       [6.9, 3.1, 4.9, 1.5],
       [5.5, 2.3, 4. , 1.3],
       [6.5, 2.8, 4.6, 1.5],
       [5.7, 2.8, 4.5, 1.3],
       [6.3, 3.3, 4.7, 1.6],
       [4.9, 2.4, 3.3, 1. ],
       [6.6, 2.9, 4.6, 1.3],
       [5.2, 2.7, 3.9, 1.4],
       [5. , 2. , 3.5, 1. ],
       [5.9, 3. , 4.2, 1.5],
       [6. , 2.2, 4. , 1. ],
       [6.1, 2.9, 4.7, 1.4],
       [5.6, 2.9, 3.6, 1.3],
       [6.7, 3.1, 4.4, 1.4],
       [5.6, 3. , 4.5, 1.5],
       [5.8, 2.7, 4.1, 1. ],
       [6.2, 2.2, 4.5, 1.5],
       [5.6, 2.5, 3.9, 1.1],
       [5.9, 3.2, 4.8, 1.8],
       [6.1, 2.8, 4. , 1.3],
       [6.3, 2.5, 4.9, 1.5],
       [6.1, 2.8, 4.7, 1.2],
       [6.4, 2.9, 4.3, 1.3],
       [6.6, 3. , 4.4, 1.4],
       [6.8, 2.8, 4.8, 1.4],
       [6.7, 3. , 5. , 1.7],
       [6. , 2.9, 4.5, 1.5],
       [5.7, 2.6, 3.5, 1. ],
       [5.5, 2.4, 3.8, 1.1],
       [5.5, 2.4, 3.7, 1. ],
       [5.8, 2.7, 3.9, 1.2],
       [6. , 2.7, 5.1, 1.6],
       [5.4, 3. , 4.5, 1.5],
       [6. , 3.4, 4.5, 1.6],
       [6.7, 3.1, 4.7, 1.5],
       [6.3, 2.3, 4.4, 1.3],
       [5.6, 3. , 4.1, 1.3],
       [5.5, 2.5, 4. , 1.3],
       [5.5, 2.6, 4.4, 1.2],
       [6.1, 3. , 4.6, 1.4],
       [5.8, 2.6, 4. , 1.2],
       [5. , 2.3, 3.3, 1. ],
       [5.6, 2.7, 4.2, 1.3],
       [5.7, 3. , 4.2, 1.2],
       [5.7, 2.9, 4.2, 1.3],
       [6.2, 2.9, 4.3, 1.3],
       [5.1, 2.5, 3. , 1.1],
       [5.7, 2.8, 4.1, 1.3]])
In [77]:
X_1.shape, sum(y==1) # Перепроверяем размеры.
Out[77]:
((50, 4), 50)
In [78]:
X[y == 1, 0] # Берем массив нулевых признаков отсносящихся к классу 1. Фактически 0ой вектор-столбец из X_1.
# Иначе говоря, идет проекция всего признакового пространства (там было 4 признака) на 0 ось.
Out[78]:
array([7. , 6.4, 6.9, 5.5, 6.5, 5.7, 6.3, 4.9, 6.6, 5.2, 5. , 5.9, 6. ,
       6.1, 5.6, 6.7, 5.6, 5.8, 6.2, 5.6, 5.9, 6.1, 6.3, 6.1, 6.4, 6.6,
       6.8, 6.7, 6. , 5.7, 5.5, 5.5, 5.8, 6. , 5.4, 6. , 6.7, 6.3, 5.6,
       5.5, 5.5, 6.1, 5.8, 5. , 5.6, 5.7, 5.7, 6.2, 5.1, 5.7])
In [79]:
colors = ['navy', 'turquoise', 'darkorange']
lw = 2

for color, i, target_name in zip(colors, [0, 1, 2], target_names): # Хитрым образом записаный цикл for.
    plt.scatter(X[y == i, 0], X[y == i, 1], color=color, alpha=.8, lw=lw, # Отрисовываем каждый класс по отдельности.
                label=target_name) # Названия наборов данных.
plt.legend(loc='best', shadow=False, scatterpoints = 1 ) # Отображение легенды.
plt.title('Исходные IRIS данные') # Название графика.
Out[79]:
Text(0.5,1,'Исходные IRIS данные')

Синий класс можно почти отделить, а вот оранжевый и берюзовый скорее нет. Как быть? Приобразуем систему координат так, что бы было видно. Обращаю внимание: исходные даные 4х мерные, а картинка 2х мерная.

Упр. Изучить другие проекции.

In [80]:
pca = PCA(n_components=2, svd_solver='full')
X_r = pca.fit_transform(X) # Сразу строим модель и преобразуем.
In [81]:
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
    plt.scatter(X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=.8, lw=lw,
                label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('PCA of IRIS dataset')
Out[81]:
Text(0.5,1,'PCA of IRIS dataset')

Теперь данные гораздо проще отделить данные друг отдруга. По крайней мере можно провести четкую границу/ломанную между классами.

Посмотрим на все это дело в 3d

In [82]:
from mpl_toolkits.mplot3d import Axes3D
from sklearn import decomposition
from sklearn import datasets
np.random.seed(5)
In [83]:
iris = datasets.load_iris()
X = iris.data
y = iris.target
In [84]:
np.choose([0, 2, 1, 2], ['a','b','c']) # Один список используется для индексации элементов из второго списка
Out[84]:
array(['a', 'c', 'b', 'c'], dtype='<U1')
In [85]:
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=60)

for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
    ax.text3D(X[y == label, 0].mean(),
              X[y == label, 1].mean() + 1.0,
              X[y == label, 2].mean(), name,
              horizontalalignment='center')
              
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=yy, cmap=plt.cm.spectral, edgecolor='k')
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-85-72de35e5ea5e> in <module>
     10 # Отрисовываем данные.
     11 yy = np.choose(y, [1, 2, 0]).astype(np.float)
---> 12 ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=yy, cmap=plt.cm.spectral, edgecolor='k')

AttributeError: module 'matplotlib.cm' has no attribute 'spectral'
In [ ]:
pca = decomposition.PCA(n_components=3)
pca.fit(X)
Xpca = pca.transform(X)

fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=90)

for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
    ax.text3D(Xpca[y == label, 0].mean(),
              Xpca[y == label, 1].mean() + 1.0,
              Xpca[y == label, 2].mean(), name,
              horizontalalignment='center')
              
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(Xpca[:, 0], Xpca[:, 1], Xpca[:, 2], c=yy, cmap=plt.cm.spectral, edgecolor='k')

Ослабим часть условий на метод главных компонент

In [ ]:
# В Методе главных компонент считалось, что разброс энергии (шума) постоянен по всем осям.
# В данном методе каждая из ортогональных осей независима.
from sklearn.decomposition import FactorAnalysis
In [86]:
fa = FactorAnalysis( n_components=2 )
X_fa = fa.fit_transform(X)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-86-dbd505457ce6> in <module>
----> 1 fa = FactorAnalysis( n_components=2 )
      2 X_fa = fa.fit_transform(X)

NameError: name 'FactorAnalysis' is not defined
In [87]:
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
    plt.scatter(X_fa[y == i, 0], X_fa[y == i, 1], color=color, alpha=.8, lw=lw,
                label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('FA of IRIS dataset')
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-87-a62d916e4f87> in <module>
      1 for color, i, target_name in zip(colors, [0, 1, 2], target_names):
----> 2     plt.scatter(X_fa[y == i, 0], X_fa[y == i, 1], color=color, alpha=.8, lw=lw,
      3                 label=target_name)
      4 plt.legend(loc='best', shadow=False, scatterpoints=1)
      5 plt.title('FA of IRIS dataset')

NameError: name 'X_fa' is not defined

Когда будет разница видна?

В методе главных компонентов оси ортонормированы. Это жесткое требование. от него можно отказатся применив следующий подход.

Еще более общий случай метода главных компонент

In [88]:
# Теперь оси необязательно ортогональны
from sklearn.decomposition import FastICA, FactorAnalysis
In [89]:
# ?rng
In [90]:
# Создадим данные
rng = np.random.RandomState(42)
S = rng.standard_t(1.5,  size=(20000, 2))  #standard_t 1.5,  size=(20000, 2) randn
S[:, 0] *= 2.
S, S.shape
Out[90]:
(array([[  0.84805472,  -0.41316453],
        [ 36.41730631,   0.53773481],
        [ -1.76643788,  -0.67843317],
        ...,
        [ -1.35307067,   1.32228283],
        [ -3.47661776, -29.51830934],
        [ -5.41834779,   0.1860394 ]]), (20000, 2))
In [91]:
SS = S / S.std()
plt.scatter(SS[:, 0], SS[:, 1], s=2, marker='o', zorder=10,
                color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
Out[91]:
(-3, 3)
In [92]:
# Смешиваем данные распределенные по разным осям
A = np.array([[1, 1], [0, 2]])  # Матрица смешивания
X = np.dot(S, A.T)  # Формируем наблюдения
XX = X / X.std()
plt.scatter(XX[:, 0], XX[:, 1], s=2, marker='o', zorder=10,
                color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
Out[92]:
(-3, 3)
In [93]:
pca = FactorAnalysis() #PCA()
S_pca_ = pca.fit(X).transform(X)

ica = FastICA(random_state=rng)
S_ica_ = ica.fit(X).transform(X)  # Estimate the sources

S_ica_ /= S_ica_.std(axis=0)
In [94]:
def plot_samples(S, axis_list=None):
    plt.scatter(S[:, 0], S[:, 1], s=2, marker='o', zorder=10,
                color='steelblue', alpha=0.5)
    if axis_list is not None:
        colors = ['orange', 'red']
        for color, axis in zip(colors, axis_list):
            axis /= axis.std()
            x_axis, y_axis = axis
            # Trick to get legend to work
            plt.plot(0.1 * x_axis, 0.1 * y_axis, linewidth=2, color=color)
            plt.quiver(0, 0, x_axis, y_axis, zorder=11, width=0.01, scale=6,
                       color=color)

    plt.hlines(0, -3, 3)
    plt.vlines(0, -3, 3)
    plt.xlim(-3, 3)
    plt.ylim(-3, 3)
    plt.xlabel('x')
    plt.ylabel('y')

plt.figure(figsize=(8, 6))
plt.subplot(2, 2, 1)
plot_samples(S / S.std())
plt.title('True Independent Sources')

axis_list = [pca.components_.T, ica.mixing_]
plt.subplot(2, 2, 2)
plot_samples(X / np.std(X), axis_list=axis_list)
legend = plt.legend(['PCA', 'ICA'], loc='upper right')
legend.set_zorder(100)

plt.title('Observations')

plt.subplot(2, 2, 3)
plot_samples(S_pca_ / np.std(S_pca_, axis=0))
plt.title('PCA recovered signals')

plt.subplot(2, 2, 4)
plot_samples(S_ica_ / np.std(S_ica_))
plt.title('ICA recovered signals')
Out[94]:
Text(0.5,1,'ICA recovered signals')

Расщепление сигнала на компоненты.

In [95]:
from scipy import signal

from sklearn.decomposition import FastICA, PCA

# #############################################################################
# Generate sample data
np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time)  # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time))  # Signal 2 : square signal
s3 = signal.sawtooth(2 * np.pi * time)  # Signal 3: saw tooth signal

S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape)  # Add noise

S /= S.std(axis=0)  # Standardize data
# Mix data
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]])  # Mixing matrix
X = np.dot(S, A.T)  # Generate observations

# Compute ICA
ica = FastICA(n_components=3)
S_ = ica.fit_transform(X)  # Reconstruct signals
A_ = ica.mixing_  # Get estimated mixing matrix

# We can `prove` that the ICA model applies by reverting the unmixing.
assert np.allclose(X, np.dot(S_, A_.T) + ica.mean_)

# For comparison, compute PCA
pca = PCA(n_components=3)
H = pca.fit_transform(X)  # Reconstruct signals based on orthogonal components

# #############################################################################
# Plot results

plt.figure(figsize=(12,12))

models = [X, S, S_, H]
names = ['Observations (mixed signal)',
         'True Sources',
         'ICA recovered signals',
         'PCA recovered signals']
colors = ['red', 'steelblue', 'orange']

for ii, (model, name) in enumerate(zip(models, names), 1):
    plt.subplot(4, 1, ii)
    plt.title(name)
    for sig, color in zip(model.T, colors):
        plt.plot(sig, color=color)
In [96]:
X.shape
Out[96]:
(2000, 3)
In [97]:
A_/(-46)
Out[97]:
array([[ 1.00046695, -0.92408907, -0.98539225],
       [ 0.95780858, -0.46065901, -1.92730583],
       [ 1.99407572, -1.3425093 , -0.96907967]])

Словари

In [98]:
from time import time

import scipy as sp

from sklearn.decomposition import MiniBatchDictionaryLearning
from sklearn.feature_extraction.image import extract_patches_2d
from sklearn.feature_extraction.image import reconstruct_from_patches_2d


try:  # SciPy >= 0.16 have face in misc
    from scipy.misc import face
    face = face(gray=True)
except ImportError:
    face = sp.face(gray=True)

# Convert from uint8 representation with values between 0 and 255 to
# a floating point representation with values between 0 and 1.
face = face / 255.

# downsample for higher speed
face = face[::2, ::2] + face[1::2, ::2] + face[::2, 1::2] + face[1::2, 1::2]
face /= 4.0
height, width = face.shape

# Distort the right half of the image
print('Distorting image...')
distorted = face.copy()
distorted[:, width // 2:] += 0.075 * np.random.randn(height, width // 2)

# Extract all reference patches from the left half of the image
print('Extracting reference patches...')
t0 = time()
patch_size = (7, 7)
data = extract_patches_2d(distorted[:, :width // 2], patch_size)
data = data.reshape(data.shape[0], -1)
data -= np.mean(data, axis=0)
data /= np.std(data, axis=0)
print('done in %.2fs.' % (time() - t0))

# #############################################################################
# Learn the dictionary from reference patches

print('Learning the dictionary...')
t0 = time()
dico = MiniBatchDictionaryLearning(n_components=100, alpha=1, n_iter=500)
V = dico.fit(data).components_
dt = time() - t0
print('done in %.2fs.' % dt)

plt.figure(figsize=(4.2, 4))
for i, comp in enumerate(V[:100]):
    plt.subplot(10, 10, i + 1)
    plt.imshow(comp.reshape(patch_size), cmap=plt.cm.gray_r,
               interpolation='nearest')
    plt.xticks(())
    plt.yticks(())
plt.suptitle('Dictionary learned from face patches\n' +
             'Train time %.1fs on %d patches' % (dt, len(data)),
             fontsize=16)
plt.subplots_adjust(0.08, 0.02, 0.92, 0.85, 0.08, 0.23)


# #############################################################################
# Display the distorted image

def show_with_diff(image, reference, title):
    """Helper function to display denoising"""
    plt.figure(figsize=(5, 3.3))
    plt.subplot(1, 2, 1)
    plt.title('Image')
    plt.imshow(image, vmin=0, vmax=1, cmap=plt.cm.gray,
               interpolation='nearest')
    plt.xticks(())
    plt.yticks(())
    plt.subplot(1, 2, 2)
    difference = image - reference

    plt.title('Difference (norm: %.2f)' % np.sqrt(np.sum(difference ** 2)))
    plt.imshow(difference, vmin=-0.5, vmax=0.5, cmap=plt.cm.PuOr,
               interpolation='nearest')
    plt.xticks(())
    plt.yticks(())
    plt.suptitle(title, size=16)
    plt.subplots_adjust(0.02, 0.02, 0.98, 0.79, 0.02, 0.2)

show_with_diff(distorted, face, 'Distorted image')

# #############################################################################
# Extract noisy patches and reconstruct them using the dictionary

print('Extracting noisy patches... ')
t0 = time()
data = extract_patches_2d(distorted[:, width // 2:], patch_size)
data = data.reshape(data.shape[0], -1)
intercept = np.mean(data, axis=0)
data -= intercept
print('done in %.2fs.' % (time() - t0))

transform_algorithms = [
    ('Orthogonal Matching Pursuit\n1 atom', 'omp',
     {'transform_n_nonzero_coefs': 1}),
    ('Orthogonal Matching Pursuit\n2 atoms', 'omp',
     {'transform_n_nonzero_coefs': 2}),
    ('Least-angle regression\n5 atoms', 'lars',
     {'transform_n_nonzero_coefs': 5}),
    ('Thresholding\n alpha=0.1', 'threshold', {'transform_alpha': .1})]

reconstructions = {}
for title, transform_algorithm, kwargs in transform_algorithms:
    print(title + '...')
    reconstructions[title] = face.copy()
    t0 = time()
    dico.set_params(transform_algorithm=transform_algorithm, **kwargs)
    code = dico.transform(data)
    patches = np.dot(code, V)

    patches += intercept
    patches = patches.reshape(len(data), *patch_size)
    if transform_algorithm == 'threshold':
        patches -= patches.min()
        patches /= patches.max()
    reconstructions[title][:, width // 2:] = reconstruct_from_patches_2d(
        patches, (height, width // 2))
    dt = time() - t0
    print('done in %.2fs.' % dt)
    show_with_diff(reconstructions[title], face,
                    title + ' (time: %.1fs)' % dt)
Distorting image...
Extracting reference patches...
done in 0.08s.
Learning the dictionary...
/data/conda/anaconda3/envs/data_analysis/lib/python3.6/site-packages/sklearn/feature_extraction/image.py:287: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  indexing_strides = arr[slices].strides
done in 8.62s.
Extracting noisy patches... 
done in 0.04s.
Orthogonal Matching Pursuit
1 atom...
done in 9.52s.
Orthogonal Matching Pursuit
2 atoms...
done in 19.84s.
Least-angle regression
5 atoms...
done in 101.48s.
Thresholding
 alpha=0.1...
done in 0.90s.
In [99]:
V.shape
Out[99]:
(100, 49)
In [ ]: