Для построения PCA счетов и нагрузок, используется рекуррентный алгоритм NIPALS, который на каждом шагу вычисляет одну компоненту. Сначала исходная матрица X преобразуется (как минимум – центрируется; см. раздел 4.3) и превращается в матрицу 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
О вычислении PCA с помощью надстройки Chemometrics рассказано в пособии Проекционные методы в системе Excel.
PLS1
Самым популярным способом для многомерной калибровки является метод проекции на латентные структуры (PLS). В этом методе проводится одновременная декомпозиция матрицы предикторов Xи матрицы откликов Y:
X=TPt+EY=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 = ttfa / ttt
5. u = qfa / q2
6. pt = ttEa / 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);
% End of pls
О вычислении PLS1 с помощью надстройки Chemometrics рассказано в пособии Проекционные методы в системе Excel.
PLS2
Для PLS2 алгоритм выглядит следующим образом. Сначала исходные матрицы X и Y преобразуют (как минимум – центрируют; см. разделе 4.3), и они превращаются в матрицы E0 и F0, a=0. Далее к ним применяет следующий алгоритм.
1. Выбрать начальный вектор u
2. wt = utEa
3. w = w / (wtw)½
4. t = Eaw
5. qt = ttFa / ttt
6. u = Faq/ qtq
7. Проверить сходимость, если нет, то идти на 2
8. pt = ttEa / 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
О вычислении PLS2 с помощью надстройки Chemometrics рассказано в пособии Проекционные методы в системе Excel.
