Центрирование и шкалирование

Часто при анализе требуется преобразовать исходные данные. Наиболее используемыми методами преобразования данных выступают центрирование и шкалирование каждой переменной на стандартное отклонение. В разделе 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

Содержание

5.2. SVD/PCA

Наиболее популярным способом сжатия данных в многомерном анализе является метод главных компонент (PCA). С математической точки зрения PCA — это декомпозиция исходной матрицы X, т.е. представление ее в виде произведения двух матриц T и P

X = TP^t + 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

Содержание

5.3 PCA/NIPALS

Для построения PCA счетов и нагрузок, используется рекуррентный алгоритм NIPALS, который на каждом шагу вычисляет одну компоненту. Сначала исходная матрица X преобразуется (как минимум – центрируется; см. раздел 4.3) и превращается в матрицу E₀, a=0. Далее применяют следующий алгоритм.

1. Выбрать начальный вектор t

2. p^t = t^t E_a / t^tt

3. p = p / (p^tp)^½

4. t = E_a p / p^tp

5. Проверить сходимость, если нет, то идти на 2

После вычисления очередной (a-ой) компоненты, полагаем t_a=t и p_a=p. Для получения следующей компоненты надо вычислить остатки E_a+1 = E_a – t p^t и применить к ним тот же алгоритм, заменив индекс 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=TP^t+E Y=UQ^t+F T=XW(P^tW)^–1

Проекция строится согласованно – так, чтобы максимизировать корреляцию между соответствующими векторами X-счетов t_aи Y-счетов u_a. Если блок данных Y включает несколько откликов (т.е. K>1), можно построить две проекции исходных данных – PLS1 и PLS2. В первом случае для каждого из откликов y_k строится свое проекционное подпространство. При этом и счетаT (U) и нагрузки P (W, Q) , зависят от того, какой отклик используется. Этот подход называется PLS1. Для метода PLS2 строится только одно проекционное пространство, которое является общим для всех откликов.

Детальное описание метода PLS приведено в этой книге Для построения PLS1 счетов и нагрузок, используется рекуррентный алгоритм. Сначала исходные матрицы X и Y центрируют

[E0, mX] = mc(X); [F0, mY] = mc(Y);

и они превращаются в матрицу E₀ и вектор f₀, a=0. Далее к ним применяет следующий алгоритм

1. w^t = f_a^t E_a

2. w = w / (w^tw)^½

3. t = E_a w

4. q = t^t f_a / t^tt

5. u = qf_a / q²

6. p^t = t^t E_a / t^tt

После вычисления очередной (a-ой) компоненты, полагаем t_a=t и p_a=p. Для получения следующей компоненты надо вычислить остатки E_a+1 = E_a – t p^t и применить к ним тот же алгоритм, заменив индекс 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 Add In рассказано в пособии Проекционные методы в системе Excel.

Содержание

PLS2

Для PLS2 алгоритм выглядит следующим образом. Сначала исходные матрицы X и Y преобразуют (как минимум – центрируют; см. разделе 4.3), и они превращаются в матрицы E₀ и F₀, a=0. Далее к ним применяет следующий алгоритм.

1. Выбрать начальный вектор u

2. w^t = u^t E_a

3. w = w / (w^tw)^½

4. t = E_a w

5. q^t = t^t F_a / t^tt

6. u = F_a q/ q^tq

7. Проверить сходимость, если нет, то идти на 2

8. p^t = t^t E_a / t^tt

После вычисления очередной (a-ой) PLS2 компоненты надо положить: t_a=t, p_a=p, w_a=w, u_a=u и q_a=q. Для получения следующей компоненты надо вычислить остатки E_a+1 = E_a –t p^t и F_a+1 = F_a – tq^t и применить к ним тот же алгоритм, заменив индекс 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 Add In рассказано в пособии Проекционные методы в системе Excel.

Содержание

Заключение

MatLab ­ это это очень популярный инструмент для анализа данных. По данным опроса, его используют до трети всех исследователей, тогда как программа the Unsrambler применяется только 16% ученых. Главным недостатком MatLab являются его высокая цена. Кроме того, MatLab хорош для рутинных расчетов. Отсутствие интерактивности делает его неудобным при выполнении поисковых, исследовательских расчетов для новых, неисследованных массивов данных.

Проблему цены решает альтернативное математическое обеспечение Chemometrics - специальная надстройка для системы Microsoft Excel. Подробнее о ней рассказано в пособии Проекционные методы в системе Excel.

Содержание