Подбор кривой ортогональной регрессии Используя анализ основных компонентов

Этот пример показывает, как использовать Анализ основных компонентов (PCA), чтобы соответствовать линейной регрессии. PCA минимизирует перпендикулярные расстояния от данных до подобранной модели. Это - линейный случай того, что является известным как Ортогональную Регрессию или Общие Наименьшие квадраты, и является соответствующим, когда нет никакого естественного различия между переменными прогноза и переменными отклика, или когда все переменные измеряются с ошибкой. Это в отличие от обычного предположения регрессии, что переменные прогноза измеряются точно, и только переменная отклика имеет ошибочный компонент.

Например, учитывая двух X и Y векторов данных, можно соответствовать строке, которая минимизирует перпендикулярные расстояния от каждой из точек (x (i), y (i)) к строке. В более общем плане, с p наблюдал переменные, можно приспособить r-dimensional гиперплоскость в p-мерном-пространстве (r <p). Выбор r эквивалентен выбору количества компонентов, чтобы сохранить в PCA. Это может быть основано на ошибке прогноза, или это может просто быть прагматический выбор уменьшать данные до управляемого количества размерностей.

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

Подбор кривой плоскости к 3-D данным

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

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 =

    0.6774   -0.0790
    0.2193    0.9707
    0.7022   -0.2269

normal = coeff(:,3)
normal =

    0.7314
   -0.0982
   -0.6749

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

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

pctExplained = roots' ./ sum(roots)
pctExplained =

    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, мы можем использовать компоненты, которые были уже вычислены. Вектор направления, который задает строку, дан коэффициентами для первого основного компонента. Вторые и третьи PC являются ортогональными к первому, и их коэффициенты задают направления, которые перпендикулярны строке. Самым простым уравнением, чтобы описать строку является meanX + t*dirVect, где t параметризовал положение вдоль строки.

dirVect = coeff(:,1)
dirVect =

    0.6774
    0.2193
    0.7022

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

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® вы могли в интерактивном режиме вращать график к другим точкам зрения, чтобы проверить, что проекции действительно перпендикулярны, и получить лучшее ощущение того, как строка соответствует данным.