ICA - statystyczna niezależność i wartości własne macierzy kowariancji

14

Obecnie tworzę różne sygnały za pomocą Matlaba, miksuję je, mnożąc je przez macierz miksowania A, a następnie próbuję odzyskać oryginalne sygnały za pomocą FastICA .

Jak dotąd odzyskane sygnały są naprawdę złe w porównaniu z oryginalnymi, czego nie oczekiwałem.

Próbuję sprawdzić, czy robię coś źle. Generowane przeze mnie sygnały to:

s1 = (-x.^2 + 100*x + 500) / 3000; % quadratic
s2 = exp(-x / 10); % -ve exponential
s3 = (sin(x)+ 1) * 0.5; % sine
s4 = 0.5 + 0.1 * randn(size(x, 2), 1); % gaussian
s5 = (sawtooth(x, 0.75)+ 1) * 0.5; % sawtooth

Oryginalne sygnały

Jednym z warunków powodzenia ICA jest to, że co najwyżej jeden sygnał jest gaussowski i zaobserwowałem to podczas generowania sygnału.

Jednak innym warunkiem jest to, że wszystkie sygnały są statystycznie niezależne.

Wiem tylko, że oznacza to, że biorąc pod uwagę dwa sygnały A i B, znajomość jednego sygnału nie daje żadnych informacji w odniesieniu do drugiego, tj .: P (A | B) = P (A), gdzie P jest prawdopodobieństwem .

Teraz moje pytanie brzmi: czy moje sygnały są statystycznie niezależne? Czy jest jakiś sposób, aby to ustalić? Może jakaś własność, której należy przestrzegać?

Inną rzeczą, którą zauważyłem, jest to, że kiedy obliczam wartości własne macierzy kowariancji (obliczonej dla macierzy zawierającej sygnały mieszane), widmo własne wydaje się pokazywać, że istnieje tylko jeden (główny) główny składnik . Co to tak naprawdę znaczy? Czy nie powinno być 5, skoro mam 5 (podobno) niezależnych sygnałów?

Na przykład przy użyciu następującej matrycy mieszania:

A =

0.2000    0.4267    0.2133    0.1067    0.0533
0.2909    0.2000    0.2909    0.1455    0.0727
0.1333    0.2667    0.2000    0.2667    0.1333
0.0727    0.1455    0.2909    0.2000    0.2909
0.0533    0.1067    0.2133    0.4267    0.2000

Wartościami własnymi są: 0.0000 0.0005 0.0022 0.0042 0.0345(tylko 4!)

Przy użyciu macierzy tożsamości jako matrycy do mieszania (tj sygnały mieszane są takie same jak oryginalne), przy czym eigenspectrum jest: 0.0103 0.0199 0.0330 0.0811 0.1762. Jest jeszcze jedna wartość znacznie większa niż reszta.

Dziękuję za pomoc

Przepraszam, jeśli odpowiedzi na moje pytania są boleśnie oczywiste, ale naprawdę jestem nowy w statystyce, ICA i Matlabie. Dzięki jeszcze raz.

EDYTOWAĆ

Mam 500 próbek każdego sygnału, w zakresie [0,2, 100], w krokach co 0,2, tj. X = 0: 0,1: 100.

Ponadto, biorąc pod uwagę model ICA: X = As + n (w tej chwili nie dodam żadnego hałasu), mam na myśli spektrum własne transpozycji X, tj. Eig (cov (X ')).

AKTUALIZACJA

Zgodnie z sugestią (patrz komentarze), wypróbowałem FastICA tylko na 2 sygnałach. Wyniki były dość dobre (patrz zdjęcie poniżej). Zastosowano matrycę mieszającą A = [0.75 0.25; 0.25 0.75]. Jednak spektrum własne 0.1657 0.7732 nadal wykazywało tylko jeden główny główny składnik.

Moje pytanie sprowadza się zatem do następującego: jakiej funkcji / równania / właściwości mogę użyć do sprawdzenia, czy pewna liczba wektorów sygnałowych jest statystycznie niezależna?

Sinus i Gaussian - FastICA

Rachel
źródło
1
Doskonałe pytanie. Zapytałem o to, skąd możemy wiedzieć, kiedy dwa sygnały są tutaj niezależne ( dsp.stackexchange.com/questions/1242/… ), ale nie posunąłem się za daleko. :-) Jestem także nowy w ICA, ale mogę rzucić trochę światła.
Spacey
@Mohammad Czy nadal interesuje Cię odpowiedź na to pytanie? Z przyjemnością naliczę nagrodę, aby przyciągnąć zainteresowanie.
Phonon
@Mohammad Głosowałem za twoim pytaniem. Mam nadzieję, że dostaniesz dobrą odpowiedź, to jest naprawdę związane z moją. Do tej pory czytałem komentarze na ten temat i dzieje się wiele statystyk, których nie rozumiem. Czy udało ci się znaleźć konkretny sposób na stwierdzenie, czy dwa sygnały są niezależne, czy nie?
Rachel
@Rachel W tej chwili jeszcze nie, ale przyjrzę się temu i dam znać. Jest to bardzo ważna koncepcja, która wydaje mi się niestety przeszklona.
Spacey
Dziękuję @Mohammad. Zgadzam się. Niezależne sygnały obserwują właściwość, że E (s1, s2) = E (s1) x E (s2), ale nie wiem, jak to właściwie obliczyć dla rzeczywistych sygnałów.
Rachel

Odpowiedzi:

8

Sygnały 3 i 5 wydają się być dość skorelowane - mają wspólną pierwszą harmoniczną. Gdybym otrzymał dwie mieszanki tych, nie byłbym w stanie ich rozdzielić, pokusiłbym się, aby umieścić wspólną harmoniczną jako jeden sygnał, a wyższe harmoniczne jako drugi sygnał. I myliłbym się! To może wyjaśniać brakującą wartość własną.

Sygnały 1 i 2 również nie wyglądają na niezależne.

Szybkim i brudnym „sprawdzeniem rozsądku” dla niezależności dwóch serii jest wykonanie (x, y) wykresu jednego sygnału względem drugiego:

plot (sig3, sig5)

a następnie wykonać ten sam wykres (x, y) z jednym tasowanym sygnałem:

indices = randperm(length(sig3))
plot(sig3(indices), sig5)

Jeśli dwie wykresy mają odmienny wygląd, twoje sygnały nie są niezależne. Mówiąc bardziej ogólnie, jeśli wykres danych (x, y) pokazuje „cechy”, nierówności itp., To zły znak.

Właściwe testy niezależności (i są to funkcje celu stosowane w pętli optymalizacyjnej ICA) obejmują na przykład wzajemną informację.

ICA odzyskuje najbardziej niezależne sygnały, których liniowe miksowanie daje dane wejściowe . Będzie działać jako metoda separacji sygnałów i odzyskać oryginalne sygnały tylko wtedy, gdy były one maksymalnie niezależne zgodnie z kryterium optymalizacji zastosowanym w implementacji ICA.

fenenety
źródło
1
Pytanie: Gdyby 5 sygnałów w jej przypadku było w rzeczywistości wszystkich niezależnych, to spodziewalibyśmy się, że NIE będzie żadnych poprawnych głównych elementów? (Innymi słowy, wszystkie wartości własne byłyby trochę takie same). Geometrycznie, mielibyśmy „chmurę” Guassiana w 5 wymiarach, zgadzasz się?
Spacey
Skontaktowałem się również z autorem ICA w sprawie usunięcia dwóch sinusoid z mieszanki, a on powiedział, że można to faktycznie zrobić za pomocą ICA. To trochę mnie dezorientuje w oparciu o to, co mówisz w odniesieniu do sygnałów 3 i 5, ponieważ (zgadzam się z tobą) wyglądają na skorelowane.
Spacey
@ Pikhenetki Wykreśliłem te wykresy, jak zasugerowałeś - a wykresy rzeczywiście wyglądają inaczej. Niestety nadal utknąłem, jak sprawdzić niepodległość. Naprawdę potrzebuję sposobu na generowanie sygnałów, które są statystycznie niezależne, aby móc ocenić wydajność FastICA.
Rachel
x1[n]x2[n]
@Mohammad Nie nagrałem własnego głosu, ale próbowałem użyć FastICA na połączeniu sygnałów sinusodialnych i gaussowskich. MYŚLĘ, że są niezależni. FastICA wypadła całkiem dobrze, ale spektrum własne było nadal dziwne. Zaktualizuję moje pytanie, aby wyświetlić wyniki.
Rachel
7

Nie jestem ekspertem od ICA, ale mogę powiedzieć trochę o niezależności.

Jak wspomniano w niektórych komentarzach, statystyczną niezależność między dwiema zmiennymi losowymi można z grubsza interpretować jako „ilość informacji, jaką zaobserwowanie jednej zmiennej daje o drugiej”.

XYXYp(x,y)XYp(x,y)=p(x)p(y)

p(x,y)

XYXYp(X=i,Y=j)=pijP(X=i)=piP(Y=j)=pj

I(X,Y)=ijpijlogpijpipj

Oto kod Matlaba, który wygeneruje dwa niezależne sygnały ze skonstruowanego rozkładu połączeń, a dwa z niezależnego rozkładu połączeń, a następnie obliczy wzajemną informację o połączeniach.

Funkcja „computeMIplugin.m” to prosta funkcja, którą napisałem, która oblicza wzajemne informacje przy użyciu powyższej formuły sumowania.

Ndist = 25;
xx = linspace(-pi, pi, Ndist);

P1 = abs(sin(xx)); P2 = abs(cos(xx)); 
P1 = P1/sum(P1); P2 = P2/sum(P2); % generate marginal distributions

%% Draw independent samples.
Nsamp = 1e4;
X = randsample(xx, Nsamp, 'true', P1);
Y = randsample(xx, Nsamp, 'true', P2);

Pj1 = P1'*P2;
computeMIplugin(Pj1)

% I get approx 8e-15 ... independent!

% Now Sample the joint distribution 
cnt = {}; cnt{1} = xx; cnt{2} = xx; % bin centers
Pj1_samp= hist3([X' Y'],cnt); Pj1_samp = Pj1_samp/sum(Pj1_samp(:));
computeMIplugin(Pj1_samp)
% I get approx .02; since we've estimated the distribution from
% samples, we don't know the true value of the MI. This is where
% a confidence interval would come in handy. We'd like to know 
% whether value of MI is significantly different from 0. 

% mean square difference between true and sampled?
% (this is small for these parameter settings... 
% depends on the sample size and # bins in the distribution).
mean( (Pj1_samp(:) - Pj1(:)).^2)

%% Draw samples that aren't independent. 

tx = linspace(0,30,Nsamp);
X = pi*sin(tx);
Y = pi*cos(tx);

% estimate the joint distribution
cnt = {}; cnt{1} = xx; cnt{2} = xx; % bin centers
Pj2= hist3([X' Y'],cnt); Pj2 = Pj2/sum(Pj2(:));
computeMIplugin(Pj2)

% I get 1.9281  - not independent!

%% make figure
figure(1); 
colormap gray
subplot(221)
imagesc(xx,xx,Pj1_samp)
title('sampled joint distribution 1')
subplot(222)
imagesc(xx,xx,Pj2)
title('sampled joint distribution 2')
subplot(223)
imagesc(xx,xx,Pj1)
title('true joint distribution 1')

Ponownie, zakłada to, że masz dobre oszacowanie wspólnego rozkładu (wraz z innymi prostymi założeniami), ale powinno to być przydatne jako ogólna zasada.

tak
źródło
To dobra odpowiedź sydeulissie dzięki, Będę musiał zajrzeć w to trochę głębiej.
Spacey,
Przede wszystkim dziękuję za długą odpowiedź, która była bardzo pouczająca. Mam tylko kilka pytań. Wspomniałeś o użyciu testu chi-kwadrat. Zajrzałem do tego i wygląda naprawdę interesująco, ale jak mogę używać tego na sygnałach? Czy nie można go zastosować tylko do danych kategorycznych?
Rachel
Ponadto do obliczenia rozkładu połączeń używasz Pj1 = P1 '* P2, prawda? Ale technicznie uważam, że nie da się tego zrobić. Być może robisz to, ponieważ zakładasz, że oryginalne sygnały są niezależne, a wynik w ten sposób się utrzymuje? Ale jak możesz obliczyć wzajemną informację - ponieważ jej wynik zależy od wspólnego podziału ..? Możliwe, że coś źle zrozumiałem, ale proszę o wyjaśnienie.
Rachel
Będę szczęśliwy - choć minie trochę czasu :).
tak,
Dziękuję @sydeulissie. Chciałbym uzyskać odpowiedź, która nie zakłada, że ​​mam wiedzę na temat wspólnej dystrybucji, proszę.
Rachel
3

Jak wspomniano powyżej, oba sygnały 3 i 5 wydają się być dość skorelowane i mają podobny okres.

Możemy myśleć o korelacji dwóch sygnałów, jeśli możemy przesunąć jedno ze źródeł w lewo lub w prawo i zwiększyć lub zmniejszyć jego amplitudę, aby pasowała do drugiego źródła. Uwaga: nie zmieniamy częstotliwości źródła, po prostu wykonujemy przesunięcie fazy i amplitudy.

W powyższym przypadku możemy przesunąć źródło 3, aby jego szczyty pokrywały się ze źródłem 5. Jest to rodzaj rzeczy, która zepsuje wydobycie źródła podczas korzystania z ICA z powodu założenia niezależności.

Uwaga : Dobrym przykładem powyższej koncepcji jest pomyślenie o dwóch falach sinusoidalnych. Oba są całkowicie deterministyczne. Jeśli oba mają tę samą częstotliwość (nawet z inną fazą), wówczas są doskonale skorelowane i ICA nie będzie w stanie ich rozdzielić. Jeśli zamiast tego mają różne częstotliwości (które nie są całkowitymi wielokrotnościami), to są niezależne i można je rozdzielić.

Poniżej znajduje się kod Matlaba, abyś sam to mógł zobaczyć

%Sine waves of equal frequency
X = 1:1000;
Y(1,:) = sin(2*pi*X*10/1000);
Y(2,:) = sin(1+2*pi*X*10/1000);

figure
subplot(3,2,1)
plot(Y(1,:))
title('Initial Source 1')
subplot(3,2,2)
plot(Y(2,:))
title('Initial Source 2')
A = [1, 2; 4, -1];
Y = A*Y;
subplot(3,2,3)
plot(Y(1,:))
title('Signal 1')
subplot(3,2,4)
plot(Y(2,:))
title('Signal 2')

Z = fastica(Y);

subplot(3,2,5)
plot(Z(1,:))
title('Source 1')
subplot(3,2,6)
plot(Z(2,:))
title('Source 2')

%Sine waves of different frequency
X = 1:1000;
Y(1,:) = sin(2*pi*X*10/1000);
Y(2,:) = sin(1+2*pi*X*8/1000);

figure
subplot(3,2,1)
plot(Y(1,:))
title('Initial Source 1')
subplot(3,2,2)
plot(Y(2,:))
title('Initial Source 2')
A = [1, 2; 4, -1];
Y = A*Y;
subplot(3,2,3)
plot(Y(1,:))
title('Signal 1')
subplot(3,2,4)
plot(Y(2,:))
title('Signal 2')

Z = fastica(Y);

subplot(3,2,5)
plot(Z(1,:))
title('Source 1')
subplot(3,2,6)
plot(Z(2,:))
title('Source 2')

Zauważ, że dla fal o tej samej częstotliwości ICA po prostu zwraca sygnały wejściowe, ale dla różnych częstotliwości zwraca oryginalne źródła.

rwolst
źródło
2

Rachel,

Z moich badań do tej pory udało mi się znaleźć coś, co nazywa się „ testem chi-kwadrat dla niezależności ”, ale nie jestem pewien, jak to działa w tej chwili, ale może warto to sprawdzić.

Spacey
źródło
Znalazłem dwa samouczki wyjaśniające, jak wykonać test chi-kwadrat: ling.upenn.edu/~clight/chisquared.htm & math.hws.edu/javamath/ryan/ChiSquare.html . Test można jednak wykonać tylko na danych kategorycznych. Nie wiem, czy można to zastosować do naszych obserwacji sygnałów.
Rachel