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

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

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?

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.

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

$X$$Y$$X$$Y$$p(x,y)$$X$$Y$$p(x,y) = p(x)p(y)$

$p(x,y)$

$X$$Y$$X$$Y$$p(X=i, Y=j) = p_{ij}$$P(X=i) = p_i$$P(Y=j) = p_j$

$I(X,Y) = \sum_i \sum_j p_{ij} \log \frac{p_{ij}}{p_i p_j}$

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.

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.

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