оглавление ДМ экономическая информатика визуальные среды - 4GL Теория и практика обработки информации

Cистема численно-математического моделирования MatLab

Примеры программ MatLab

В этом разделе приведены наиболее употребительные алгоритмы, используемые при анализе многомерных данных. Рассмотрены как простейшие методы преобразования данных центрирование и шкалирование, так и алгоритмы для анализа данных — 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 = 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

PCA/NIPALS

Для построения PCA счетов и нагрузок, используется рекуррентный алгоритм NIPALS, который на каждом шагу вычисляет одну компоненту. Сначала исходная матрица X преобразуется (как минимум – центрируется) и превращается в матрицу 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^tf_a / t^tt

5. u = qf_a / q²

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

Приведем код этого алгоритма, взятый из книги

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 преобразуют (как минимум – центрируют), и они превращаются в матрицы 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^tF_a / t^tt

6. u = F_a q/ q^tq

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

8. p^t = t^tE_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₊₁ = 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.

оглавление ДМ экономическая информатика визуальные среды - 4GL Теория и практика обработки информации

Знаете ли Вы, что теоретическая модель - это математическая модель, описывающая структуру исследуемого объекта в общем виде, без спецификации конкретных числовых значений параметров.

НОВОСТИ ФОРУМА

Рыцари теории эфира