В этом примере показано, как использовать Анализ основных компонентов (PCA), чтобы соответствовать линейной регрессии. PCA минимизирует перпендикулярные расстояния от данных до подобранной модели. Это - линейный случай того, что является известным как Ортогональную Регрессию или Общие Наименьшие квадраты, и является соответствующим, когда нет никакого естественного различия между переменными прогноза и переменными отклика, или когда все переменные измеряются с ошибкой. Это в отличие от обычного предположения регрессии, что переменные предикторы измеряются точно, и только переменная отклика имеет ошибочный компонент.
Например, учитывая двух X и Y векторов данных, можно соответствовать линии, которая минимизирует перпендикулярные расстояния от каждой из точек (x (i), y (i)) к линии. В более общем плане, с p наблюдал переменные, можно приспособить r-dimensional гиперплоскость в p-мерном-пространстве (r <p). Выбор r эквивалентен выбору количества компонентов, чтобы сохранить в PCA. Это может быть основано на ошибке прогноза, или это может просто быть прагматический выбор уменьшать данные до управляемого количества размерностей.
В этом примере мы соответствуем плоскости и линии через некоторые данные по трем наблюдаемым переменным. Легко сделать то же самое для любого количества переменных, и для любой размерности модели, несмотря на то, что визуализация подгонки в более высоких размерностях, очевидно, не была бы прямой.
Во-первых, мы генерируем некоторые 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 = 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
Это - все, которое существует к подгонке. Но давайте посмотрим ближе на результаты и построим подгонку наряду с данными.
Поскольку первые два компонента объясняют такое большое отклонение в данных, как возможно с двумя размерностями, плоскость является лучшей 2D линейной аппроксимацией к данным. Эквивалентно, третий компонент объясняет наименьшее количество объема изменения данных, и это - остаточный член в регрессии. Скрытые корни (или собственные значения) от 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);
Зеленые точки выше плоскости, красные точки ниже.
Подбор кривой прямой линии к данным еще более прост, и из-за вложенного свойства PCA, мы можем использовать компоненты, которые были уже вычислены. Вектор направления, который задает линию, дан коэффициентами для первого основного компонента. Вторые и третьи PC являются ортогональными к первому, и их коэффициенты задают направления, которые перпендикулярны линии. Самым простым уравнением, чтобы описать линию является meanX + t*dirVect
, где t
параметризовал положение вдоль линии.
dirVect = coeff(:,1)
dirVect = 3×1
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®
окно рисунка, вы могли в интерактивном режиме вращать график к другим точкам зрения, чтобы проверить, что проекции действительно перпендикулярны, и получить лучшее ощущение того, как линия соответствует данным.