Ответ 1
Вот краткое пошаговое руководство. Сначала мы создаем матрицу ваших скрытых переменных (или "факторов" ). Он имеет 100 наблюдений и есть два независимых фактора.
>> factors = randn(100, 2);
Теперь создайте матрицу нагрузок. Это будет отображать скрытые переменные на ваши наблюдаемые переменные. Скажите, что ваши наблюдаемые переменные имеют четыре функции. Тогда ваша матрица нагрузок должна быть 4 x 2
>> loadings = [
1 0
0 1
1 1
1 -1 ];
Это говорит о том, что первая наблюдаемая переменная нагрузка на первый коэффициент, вторая нагрузка на второй коэффициент, третья переменная нагрузка на сумму факторов и четвертая переменная нагрузка на разность факторов.
Теперь создайте свои наблюдения:
>> observations = factors * loadings' + 0.1 * randn(100,4);
Я добавил небольшое количество случайных шумов для имитации экспериментальной ошибки. Теперь мы выполняем PCA с помощью функции pca из панели инструментов статистики:
>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);
Переменная score - это массив показателей основного компонента. Они будут ортогональными по построению, которые вы можете проверить -
>> corr(score)
ans =
1.0000 0.0000 0.0000 0.0000
0.0000 1.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000
0.0000 0.0000 0.0000 1.0000
Комбинация score * coeff' будет воспроизводить центральную версию ваших наблюдений. Среднее значение mu вычитается до выполнения PCA. Чтобы воспроизвести ваши первоначальные наблюдения, вам нужно добавить их обратно,
>> reconstructed = score * coeff' + repmat(mu, 100, 1);
>> sum((observations - reconstructed).^2)
ans =
1.0e-27 *
0.0311 0.0104 0.0440 0.3378
Чтобы получить приблизительное количество исходных данных, вы можете начать отбрасывать столбцы из вычисленных основных компонентов. Чтобы понять, какие столбцы нужно отбросить, мы рассмотрим переменную explained
>> explained
explained =
58.0639
41.6302
0.1693
0.1366
Записи указывают, какой процент от дисперсии объясняется каждым из основных компонентов. Мы можем ясно видеть, что первые два компонента более значительны, чем два вторых (они объясняют более 99% дисперсии между ними). Использование первых двух компонентов для восстановления наблюдений дает приближение ранга-2,
>> approximationRank2 = score(:,1:2) * coeff(:,1:2)' + repmat(mu, 100, 1);
Теперь мы можем попытаться выполнить печать:
>> for k = 1:4
subplot(2, 2, k);
hold on;
grid on
plot(approximationRank2(:, k), observations(:, k), 'x');
plot([-4 4], [-4 4]);
xlim([-4 4]);
ylim([-4 4]);
title(sprintf('Variable %d', k));
end
Мы получаем почти идеальное воспроизведение оригинальных наблюдений. Если бы мы хотели получить более грубое приближение, мы могли бы просто использовать первый главный компонент:
>> approximationRank1 = score(:,1) * coeff(:,1)' + repmat(mu, 100, 1);
и нарисуйте его,
>> for k = 1:4
subplot(2, 2, k);
hold on;
grid on
plot(approximationRank1(:, k), observations(:, k), 'x');
plot([-4 4], [-4 4]);
xlim([-4 4]);
ylim([-4 4]);
title(sprintf('Variable %d', k));
end
На этот раз реконструкция не так хороша. Это потому, что мы сознательно построили наши данные, чтобы иметь два фактора, и мы только восстанавливаем их из одного из них.
Обратите внимание, что, несмотря на наводящее сходство между тем, как мы построили исходные данные и их воспроизведение,
>> observations = factors * loadings' + 0.1 * randn(100,4);
>> reconstructed = score * coeff' + repmat(mu, 100, 1);
нет никакого соответствия между factors и score, или между loadings и coeff. Алгоритм PCA не знает ничего о том, как строятся ваши данные - он просто пытается объяснить как можно большую общую дисперсию с каждым последующим компонентом.
Пользователь @Mari спросил в комментариях, как она могла бы построить ошибку восстановления в зависимости от количества основных компонентов. Используя переменную explained выше, это довольно просто. Я буду генерировать некоторые данные с более интересной структурой факторов, чтобы проиллюстрировать эффект -
>> factors = randn(100, 20);
>> loadings = chol(corr(factors * triu(ones(20))))';
>> observations = factors * loadings' + 0.1 * randn(100, 20);
Теперь все наблюдения накладываются на значительный общий фактор, с другими факторами уменьшения значения. Мы можем получить разложение PCA по-прежнему
>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);
и укажите процент объясненной дисперсии следующим образом:
>> cumexplained = cumsum(explained);
cumunexplained = 100 - cumexplained;
plot(1:20, cumunexplained, 'x-');
grid on;
xlabel('Number of factors');
ylabel('Unexplained variance')
