В этом разделе приведены наиболее употребительные алгоритмы,
используемые при анализе многомерных
данных. Рассмотрены как простейшие
методы преобразования данных центрирование и шкалирование, так и
алгоритмы для анализа данных — PCA, PLS.
Центрирование и шкалирование
Часто при
анализе требуется преобразовать
исходные данные. Наиболее
используемыми методами преобразования
данных выступают центрирование и
шкалирование каждой переменной на
стандартное отклонение. В разделе 4.3 приводился код функции для
центрирования матрицы. Поэтому ниже показан только код функции, которая
шкалирует данные. Обратите внимание, что исходная матрица должна быть
центрирована
function Xs = scaling(X)
% scaling: the output matrix is Xs
% matrix X must be centered
Xs = X * inv(diag(std(X)));
%end of scaling
SVD/PCA
Наиболее популярным способом
сжатия данных в многомерном анализе
является метод главных компонент (PCA). С
математической точки зрения PCA — это
декомпозиция исходной матрицы X, т.е.
представление ее в виде произведения
двух матриц T и P
X = TPt + E
Матрица T называется
матрицей счетов (scores) , матрица P —
матрицей нагрузок (loadings), а E —
матрицей остатков.
Простейший
способ найти матрицы T и P —
использовать SVD разложение через стандартную
функцию MatLab, называемую svd.
function [T, P] = pcasvd(X)
% pcasvd: calculates PCA components.
% The output matrices are T and P.
% T contains scores
% P contains loadings
[U,D,V] = svd(X);
T = U * D;
P = V;
%end of pcasvd
PCA/NIPALS
Для построения PCA счетов и нагрузок,
используется рекуррентный алгоритм NIPALS,
который на каждом шагу вычисляет одну
компоненту. Сначала исходная матрица X
преобразуется (как минимум – центрируется)
и превращается в матрицу E0, a=0.
Далее применяют следующий алгоритм.
1. Выбрать начальный вектор t
2. pt = ttEa
/ ttt
3. p = p / (ptp)½
4. t = Eap / ptp
5. Проверить сходимость, если нет,
то идти на 2
После вычисления очередной (a-ой)
компоненты, полагаем ta=t
и pa=p. Для получения
следующей компоненты надо вычислить
остатки Ea+1 = Ea
– tpt и применить к ним тот
же алгоритм, заменив индекс a на a+1.
Код
алгоритма NIPALS может быть написан и
самими читателями, в данном же пособии
авторы приводят свой вариант. При
расчете PCA, можно вводить число главных
компонент (переменная numberPC).
Если же не известно, сколько необходимо
компонент, следует написать в командной
строке [P,T] =
pcanipals (X) и тогда программа задаст
число компонент равным наименьшему из
показателей размерности исходной
матрицы X.
function [T, P] = pcanipals(X, numberPC)
% pcanipals: calculates PCA components.
% The output matrices are T and P.
% T contains scores
% P contains loadings
% calculation of number of
components
[X_r, X_c] = size(X); P=[]; T=[];
if lenfth(numberPC) > 0 pc = numberPC{1};
elseif (length(numberPC) == 0) & X_r < X_c pc = X_r;
else pc = X_c;
end;
% calculation of scores and loadings for each component
for k = 1:pc P1 = rand(X_c,
1); T1 = X * P1; d0 = T1'*T1; P1 = (T1' * X/(T1'
* T1))'; P1 = P1/norm(P1); T1 = X * P1; d = T1' * T1;
while d - d0
> 0.0001; P1 = (T1' * X/(T1' * T1)); P1 = P1/norm(P1); T1 = X * P1; d0
= T1'*T1; P1 = (T1' * X/(T1' * T1)); P1 = P1/norm(P1); T1 = X * P1; d
= T1'*T1; end
X = X - T1 * P1;
P = cat(1, P, P1'); T = [T,T1];
end
Самым популярным способом для
многомерной калибровки является метод
проекции на латентные структуры (PLS). В
этом методе проводится одновременная
декомпозиция матрицы предикторов X и
матрицы откликов Y:
X=TPt+E Y=UQt+FT=XW(PtW)–1
Проекция строится согласованно – так,
чтобы максимизировать корреляцию
между соответствующими векторами X-счетов
ta и Y-счетов ua.
Если блок данных Y включает
несколько откликов (т.е. K>1), можно
построить две проекции исходных данных
– PLS1 и PLS2. В первом случае для каждого из
откликов yk строится свое
проекционное подпространство. При этом
и счета T (U) и нагрузки P (W, Q)
, зависят от того, какой отклик
используется. Этот подход называется PLS1.
Для метода PLS2 строится только одно
проекционное пространство, которое
является общим для всех откликов.
Детальное
описание метода PLS приведено в этой
книге Для построения PLS1 счетов и
нагрузок, используется рекуррентный
алгоритм. Сначала исходные матрицы X
и Y центрируют
[E0, mX] = mc(X);
[F0, mY] = mc(Y);
и они
превращаются в матрицу E0 и
вектор f0, a=0. Далее к ним
применяет следующий алгоритм
1. wt = fatEa
2. w = w / (wtw)½
3. t = Eaw
4. q = tt fa /
ttt
5. u = qfa / q2
6. pt = tt Ea
/ ttt
После вычисления очередной (a-ой)
компоненты, полагаем ta=t
и pa=p. Для получения
следующей компоненты надо вычислить
остатки Ea+1 = Ea
– tpt и применить к ним тот
же алгоритм, заменив индекс a на a+1.
function [w, t, u, q, p] = pls(x, y)
%PLS: calculates a PLS component.
%The output vectors are w, t, u, q and p.
%
% Choose a vector from y as starting vector u.
u = y(:, 1);
% The convergence criterion is set very high.
kri = 100;
% The commands from here to end are repeated
until convergence.
while (kri > 1e - 10)
% Each starting vector u is saved as uold.
uold = u; w = (u' * x)'; w = w/norm(w);
t = x * w; q = (t' * y)'/(t' * t);
u = y * q/(q' * q);
% The convergence criterion is the norm of
u-uold divided by the norm of u.
kri = norm(uold - u)/norm(u);
end;
% After convergence, calculate p.
p = (t' * x)'/(t' * t);
Для PLS2 алгоритм выглядит следующим
образом. Сначала исходные матрицы X
и Y преобразуют (как минимум – центрируют), и
они превращаются в матрицы E0
и F0, a=0. Далее к ним
применяет следующий алгоритм.
1. Выбрать начальный вектор u
2. wt = utEa
3. w = w / (wtw)½
4. t = Eaw
5. qt = tt Fa
/ ttt
6. u = Faq/ qtq
7. Проверить сходимость, если нет,
то идти на 2
8. pt = tt Ea
/ ttt
После вычисления очередной (a-ой)
PLS2 компоненты надо положить: ta=t,
pa=p, wa=w, ua=u
и qa=q. Для получения
следующей компоненты надо вычислить
остатки Ea+1 = Ea
– t pt и Fa+1 =
Fa – tqt и
применить к ним тот же алгоритм,
заменив индекс a на a+1.
Приведем код, которой также
заимствован из из книги.
function [W, T, U, Q, P, B, SS] = plsr(x, y,
a)
% PLS: calculates a PLS component.
% The output matrices are W, T, U, Q and P.
% B contains the regression coefficients and SS the sums of
% squares for the residuals.
% a is the numbers of components.
%
% For a components: use all commands to end.
for i=1:a
% Calculate the sum of squares. Use the
function ss.
sx = [sx; ss(x)];
sy = [sy; ss(y)];
% Use the function pls to calculate one
component.
[w, t, u, q, p] = pls(x, y);
% Calculate the residuals.
x = x - t * p';
y = y - t * q';
% Save the vectors in matrices.
W = [W w];
T = [T t];
U = [U u];
Q = [Q q];
P = [P p];
end;
% Calculate the regression coefficients
after the loop.
B=W*inv(P'*W)*Q';
% Add the final residual SS to the sum of
squares vectors.
sx=[sx; ss(x)];
sy=[sy; ss(y)];
% Make a matrix of the ss vectors for X and
Y.
SS = [sx sy];
%Calculate the fraction of SS used.
[a, b] = size(SS);
tt = (SS * diag(SS(1,:).^(-1)) - ones(a, b)) * (-1)
%End of plsr
function [ss] = ss(x)
%SS: calculates the sum of squares of a
matrix X.
%
ss=sum(sum(x. * x));
%End of ss
Знаете ли Вы, что методы Рунге - Кутты - это важное семейство численных алгоритмов решения обыкновенных дифференциальных уравнений и их систем. Данные итеративные методы явного и неявного приближённого вычисления были разработаны около 1900 года немецкими математиками К. Рунге и М. В. Куттой. Формально, методом Рунге - Кутты является модифицированный и исправленный метод Эйлера, они представляют собой схемы второго порядка точности. Существуют стандартные схемы третьего порядка, не получившие широкого распространения. Наиболее часто используется и реализована в различных математических пакетах (Maple, MathCAD, Maxima) стандартная схема четвёртого порядка. Иногда при выполнении расчётов с повышенной точностью применяются схемы пятого и шестого порядков. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями. Методы седьмого порядка должны иметь по меньшей мере девять стадий, в схему восьмого порядка входит 11 стадий. Хотя схемы девятого порядка не имеют большой практической значимости, неизвестно, сколько стадий необходимо для достижения этого порядка точности. Аналогичная задача существует для схем десятого и более высоких порядков.