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

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

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

  1. Центрирование и шкалирование
  2. SVD/PCA
  3. PCA/NIPALS
  4. PLS1
  5. PLS2

В этом разделе приведены наиболее употребительные алгоритмы, используемые при анализе многомерных данных. Рассмотрены как простейшие методы преобразования данных центрирование и шкалирование, так и алгоритмы для анализа данных — 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 = tt Ea / ttt 
3. p = p / (ptp)½
4. t = Ea p / ptp
5. Проверить сходимость, если нет, то идти на 2

После вычисления очередной (a-ой) компоненты, полагаем ta=t и pa=p. Для получения следующей компоненты надо вычислить остатки Ea+1 = Eat pt и применить к ним тот же алгоритм, заменив индекс 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+E            Y=UQt+F                T=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 = fat Ea 
2. w = w / (wtw)½ 
3. t = Ea w
4. q = tt fa / ttt 
5. u = qfa / q2 
6. pt = tt Ea / ttt 

После вычисления очередной (a-ой) компоненты, полагаем ta=t и pa=p. Для получения следующей компоненты надо вычислить остатки Ea+1 = Eat pt и применить к ним тот же алгоритм, заменив индекс 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 преобразуют (как минимум – центрируют), и они превращаются в матрицы E0 и F0, a=0. Далее к ним применяет следующий алгоритм.

1. Выбрать начальный вектор u 
2. wt = ut Ea 
3. w = w / (wtw)½ 
4. t = Ea w
5. qt = tt Fa / ttt 
6. u = Fa q/ 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 = Fatqt  и применить к ним тот же алгоритм, заменив индекс 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   Теория и практика обработки информации

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

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

Форум Рыцари теории эфира


Рыцари теории эфира
 10.11.2021 - 12:37: ПЕРСОНАЛИИ - Personalias -> WHO IS WHO - КТО ЕСТЬ КТО - Карим_Хайдаров.
10.11.2021 - 12:36: СОВЕСТЬ - Conscience -> РАСЧЕЛОВЕЧИВАНИЕ ЧЕЛОВЕКА. КОМУ ЭТО НАДО? - Карим_Хайдаров.
10.11.2021 - 12:36: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от д.м.н. Александра Алексеевича Редько - Карим_Хайдаров.
10.11.2021 - 12:35: ЭКОЛОГИЯ - Ecology -> Биологическая безопасность населения - Карим_Хайдаров.
10.11.2021 - 12:34: ВОЙНА, ПОЛИТИКА И НАУКА - War, Politics and Science -> Проблема государственного терроризма - Карим_Хайдаров.
10.11.2021 - 12:34: ВОЙНА, ПОЛИТИКА И НАУКА - War, Politics and Science -> ПРАВОСУДИЯ.НЕТ - Карим_Хайдаров.
10.11.2021 - 12:34: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Вадима Глогера, США - Карим_Хайдаров.
10.11.2021 - 09:18: НОВЫЕ ТЕХНОЛОГИИ - New Technologies -> Волновая генетика Петра Гаряева, 5G-контроль и управление - Карим_Хайдаров.
10.11.2021 - 09:18: ЭКОЛОГИЯ - Ecology -> ЭКОЛОГИЯ ДЛЯ ВСЕХ - Карим_Хайдаров.
10.11.2021 - 09:16: ЭКОЛОГИЯ - Ecology -> ПРОБЛЕМЫ МЕДИЦИНЫ - Карим_Хайдаров.
10.11.2021 - 09:15: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Екатерины Коваленко - Карим_Хайдаров.
10.11.2021 - 09:13: ВОСПИТАНИЕ, ПРОСВЕЩЕНИЕ, ОБРАЗОВАНИЕ - Upbringing, Inlightening, Education -> Просвещение от Вильгельма Варкентина - Карим_Хайдаров.
Bourabai Research - Технологии XXI века Bourabai Research Institution