Подбор плоскости для 3-D данных

Во-первых, мы генерируем некоторые трехмерные нормальные данные для примера. Две переменные довольно сильно коррелированы.

rng(5,'twister');
X = mvnrnd([0 0 0], [1 .2 .7; .2 1 0; .7 0 1],50);
plot3(X(:,1),X(:,2),X(:,3),'bo');
grid on;
maxlim = max(abs(X(:)))*1.1;
axis([-maxlim maxlim -maxlim maxlim -maxlim maxlim]);
axis square
view(-9,12);

Затем мы подгоняем плоскость к данным с помощью PCA. Коэффициенты для первых двух главных компонентов определяют векторы, которые образуют базис для плоскости. Третий ПК ортогональен первым двум, и его коэффициенты определяют нормальный вектор плоскости.

[coeff,score,roots] = pca(X);
basis = coeff(:,1:2)

basis = 3×2

    0.6774   -0.0790
    0.2193    0.9707
    0.7022   -0.2269

normal = coeff(:,3)

normal = 3×1

    0.7314
   -0.0982
   -0.6749

Это все, что есть до конца. Но давайте посмотрим поближе на результаты и построим график соответствия вместе с данными.

Поскольку первые два компонента объясняют такое большое различие в данных, как возможно с двумя размерами, самолет - лучшее 2-е линейное приближение к данным. Эквивалентно, третий компонент объясняет наименьшее количество вариаций в данных и является термином ошибки в регрессии. Скрытые корни (или собственные значения) из PCA определяют величину объясненной дисперсии для каждого компонента.

pctExplained = roots' ./ sum(roots)

pctExplained = 1×3

    0.6226    0.2976    0.0798

Первые две координаты баллов главного компонента дают проекцию каждой точки на плоскость в системе координат плоскости. Чтобы получить координаты аппроксимированных точек в терминах исходной системы координат, мы умножаем каждый вектор коэффициента PC на соответствующую оценку и добавляем обратно в среднее значение данных. Остатки являются просто исходными данными минус совпадающие точки.

[n,p] = size(X);
meanX = mean(X,1);
Xfit = repmat(meanX,n,1) + score(:,1:2)*coeff(:,1:2)';
residuals = X - Xfit;

Уравнение аппроксимируемой плоскости, удовлетворяемое каждой из аппроксимируемых точек в Xfitявляется ([x1 x2 x3] - meanX)*normal = 0. Плоскость проходит через точку meanX, и его перпендикулярное расстояние до начала координат равно meanX*normal. Перпендикулярное расстояние от каждой точки в X к плоскости, т.е. норме остатков, - скалярное произведение каждой центрированной точки с нормалью к плоскости. Подогнанная плоскость минимизирует сумму квадратичных ошибок.

error = abs((X - repmat(meanX,n,1))*normal);
sse = sum(error.^2)

sse = 15.5142

Для визуализации аппроксимации можно построить график плоскости, исходных данных и их проекции на плоскость.

[xgrid,ygrid] = meshgrid(linspace(min(X(:,1)),max(X(:,1)),5), ...
                         linspace(min(X(:,2)),max(X(:,2)),5));
zgrid = (1/normal(3)) .* (meanX*normal - (xgrid.*normal(1) + ygrid.*normal(2)));
h = mesh(xgrid,ygrid,zgrid,'EdgeColor',[0 0 0],'FaceAlpha',0);

hold on
above = (X-repmat(meanX,n,1))*normal < 0;
below = ~above;
nabove = sum(above);
X1 = [X(above,1) Xfit(above,1) nan*ones(nabove,1)];
X2 = [X(above,2) Xfit(above,2) nan*ones(nabove,1)];
X3 = [X(above,3) Xfit(above,3) nan*ones(nabove,1)];
plot3(X1',X2',X3','-', X(above,1),X(above,2),X(above,3),'o', 'Color',[0 .7 0]);
nbelow = sum(below);
X1 = [X(below,1) Xfit(below,1) nan*ones(nbelow,1)];
X2 = [X(below,2) Xfit(below,2) nan*ones(nbelow,1)];
X3 = [X(below,3) Xfit(below,3) nan*ones(nbelow,1)];
plot3(X1',X2',X3','-', X(below,1),X(below,2),X(below,3),'o', 'Color',[1 0 0]);

hold off
maxlim = max(abs(X(:)))*1.1;
axis([-maxlim maxlim -maxlim maxlim -maxlim maxlim]);
axis square
view(-9,12);

Зеленые точки находятся над плоскостью, красные - ниже.

Подбор линии для 3-D данных

Подгонка прямой линии к данным еще проще, и из-за свойства вложенности PCA можно использовать уже вычисленные компоненты. Вектор направления, который определяет линию, задается коэффициентами для первого главного компонента. Второй и третий ПК ортогональны первому, и их коэффициенты определяют направления, перпендикулярные прямой. Самое простое уравнение для описания линии meanX + t*dirVect, где t параметризует положение вдоль линии.

dirVect = coeff(:,1)

dirVect = 3×1

    0.6774
    0.2193
    0.7022

Первая координата баллов главного компонента дает проекцию каждой точки на прямую. Как и в случае 2-D аппроксимации, векторы коэффициентов ПК, умноженные на оценки, дают подогнанные точки в исходной системе координат.

Xfit1 = repmat(meanX,n,1) + score(:,1)*coeff(:,1)';

Постройте график линии, исходных данных и их проекции на линию.

t = [min(score(:,1))-.2, max(score(:,1))+.2];
endpts = [meanX + t(1)*dirVect'; meanX + t(2)*dirVect'];
plot3(endpts(:,1),endpts(:,2),endpts(:,3),'k-');

X1 = [X(:,1) Xfit1(:,1) nan*ones(n,1)];
X2 = [X(:,2) Xfit1(:,2) nan*ones(n,1)];
X3 = [X(:,3) Xfit1(:,3) nan*ones(n,1)];
hold on
plot3(X1',X2',X3','b-', X(:,1),X(:,2),X(:,3),'bo');
hold off
maxlim = max(abs(X(:)))*1.1;
axis([-maxlim maxlim -maxlim maxlim -maxlim maxlim]);
axis square
view(-9,12);
grid on

Хотя многие проекции на этом графике не перпендикулярны линии, это просто потому, что мы строим график 3-D данных в двух измерениях. В прямом эфире MATLAB® в окне рисунка можно в интерактивном режиме повернуть график в различные точки зрения, чтобы убедиться, что проекции действительно перпендикулярны, и получить лучшее представление о том, как линия подходит к данным.