Matlab-analiza PCA i rekonstrukcja danych wielowymiarowych

Oto krótki przegląd. Najpierw tworzymy macierz ukrytych zmiennych (lub"czynników"). Ma 100 obserwacji i istnieją dwa niezależne czynniki.

>> factors = randn(100, 2);

Teraz Utwórz macierz obciążenia. To będzie mapować Ukryte zmienne NA obserwowanych zmiennych. Powiedzmy, że zaobserwowane zmienne mają cztery cechy. Następnie macierz obciążenia musi być 4 x 2

>> loadings = [
      1   0
      0   1
      1   1
      1  -1   ];

To mówi, że pierwsze zaobserwowane obciążenia zmienne na pierwszym współczynniku, drugie obciążenia na drugi czynnik, trzeci ładunek zmienny na sumę czynników i czwarty ładunek zmienny na różnicę czynników.

Teraz stwórz swoje obserwacje:

>> observations = factors * loadings' + 0.1 * randn(100,4);

Dodałem niewielką ilość losowego szumu, aby symulować eksperymentalny błąd. Teraz wykonujemy PCA używając funkcji pca z narzędzia stats toolbox:

>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);

Zmienna score jest tablicą głównych wyników składowych. Będą prostopadłe do konstrukcji, co można sprawdzić -

>> 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

Kombinacja score * coeff' odtworzy wyśrodkowaną wersję twoich obserwacji. Średnia mu jest odejmowana przed wykonaniem PCA. Aby odtworzyć swoje oryginalne obserwacje, musisz dodać je z powrotem,

>> reconstructed = score * coeff' + repmat(mu, 100, 1);
>> sum((observations - reconstructed).^2)
ans =
   1.0e-27 *
    0.0311    0.0104    0.0440    0.3378

Aby uzyskać przybliżenie do oryginalnych danych, możesz zacząć opuszczać kolumny z obliczonych głównych składników. Aby dowiedzieć się, które kolumny należy upuścić, badamy zmienną explained

>> explained
explained =
   58.0639
   41.6302
    0.1693
    0.1366

Wpisy mówią jaki procent wariancja jest wyjaśniona przez każdy z głównych składników. Możemy wyraźnie zobaczyć, że dwa pierwsze składniki są bardziej znaczące niż dwa drugie (wyjaśniają ponad 99% wariancji między nimi). Użycie dwóch pierwszych składników do rekonstrukcji obserwacji daje przybliżenie rank-2,

>> approximationRank2 = score(:,1:2) * coeff(:,1:2)' + repmat(mu, 100, 1);

Możemy teraz spróbować spiskować:

>> 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

Otrzymujemy niemal doskonałą reprodukcję oryginalnych obserwacji. Gdybyśmy chcieli grubszego przybliżenie, możemy po prostu użyć pierwszego głównego komponentu:

>> approximationRank1 = score(:,1) * coeff(:,1)' + repmat(mu, 100, 1);

I wykreślić,

>> 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

Tym razem rekonstrukcja nie jest zbyt dobra. To dlatego, że celowo skonstruowaliśmy nasze dane tak, aby miały dwa czynniki, a rekonstruujemy je tylko z jednego z nich.

>> observations  = factors * loadings'  +  0.1 * randn(100,4);
>> reconstructed = score   * coeff'     +  repmat(mu, 100, 1);

Nie musi być żadnych korespondencja pomiędzy factors a score, lub między loadings a coeff. Algorytm PCA nie wie nic o sposobie konstruowania danych - stara się jedynie wyjaśnić jak najwięcej całkowitej wariancji z każdym kolejnym składnikiem.

Użytkownik @ Mari zapytał w komentarzach, jak może wykreślić błąd rekonstrukcji jako funkcję liczby głównych komponentów. Użycie powyższej zmiennej explained jest dość proste. Wygeneruję trochę danych z bardziej interesująca struktura czynnika w celu zilustrowania efektu -

>> factors = randn(100, 20);
>> loadings = chol(corr(factors * triu(ones(20))))';
>> observations = factors * loadings' + 0.1 * randn(100, 20);

Teraz wszystkie obserwacje obciążają istotny wspólny czynnik, z innymi czynnikami o malejącym znaczeniu. Możemy uzyskać rozkład PCA jak wcześniej

>> [coeff, score, latent, tsquared, explained, mu] = pca(observations);

I wykreślić procent wyjaśnionej wariancji w następujący sposób,

>> cumexplained = cumsum(explained);
   cumunexplained = 100 - cumexplained;
   plot(1:20, cumunexplained, 'x-');
   grid on;
   xlabel('Number of factors');
   ylabel('Unexplained variance')