ICA - статистическая независимость и собственные значения ковариационной матрицы

14

В настоящее время я создаю различные сигналы с использованием Matlab, смешиваю их, умножая их на матрицу микширования A, а затем пытаюсь вернуть исходные сигналы с помощью FastICA .

Пока что восстановленные сигналы действительно плохи по сравнению с оригинальными, что оказалось не таким, как я ожидал.

Я пытаюсь понять, делаю ли я что-то не так. Сигналы, которые я генерирую, следующие:

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

Оригинальные сигналы

Одним из условий успеха ICA является то, что не более одного сигнала является гауссовским, и я наблюдал это при генерации моего сигнала.

Однако другое условие заключается в том, что все сигналы статистически независимы.

Все, что я знаю, это то, что это означает, что, учитывая два сигнала A & B, знание одного сигнала не дает никакой информации относительно другого, то есть: P (A | B) = P (A), где P - вероятность .

Теперь мой вопрос заключается в следующем: являются ли мои сигналы статистически независимыми? Есть ли способ, которым я могу определить это? Возможно, какое-то свойство, которое необходимо соблюдать?

Еще одна вещь, которую я заметил, это то, что когда я вычисляю собственные значения ковариационной матрицы (рассчитанной для матрицы, содержащей смешанные сигналы), собственный спектр, похоже, показывает, что существует только один (основной) главный компонент . Что это на самом деле значит? Разве не должно быть 5, так как у меня есть 5 (предположительно) независимых сигналов?

Например, при использовании следующей матрицы смешения:

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

Собственные значения: 0.0000 0.0005 0.0022 0.0042 0.0345(только 4!)

При использовании матрицы тождественности в качестве матрицы смешения (т.е. смешанные сигналы , такие же , как и оригинальные), то есть спектр собственных: 0.0103 0.0199 0.0330 0.0811 0.1762. Там все еще есть одно значение, намного большее, чем остальные.

Спасибо за помощь.

Я извиняюсь, если ответы на мои вопросы до боли очевидны, но я действительно плохо знаком со статистикой, ICA и Matlab. Еще раз спасибо.

РЕДАКТИРОВАТЬ

У меня есть 500 выборок каждого сигнала в диапазоне [0,2, 100] с шагом 0,2, то есть х = 0: 0,1: 100.

Кроме того, учитывая модель ICA: X = As + n (в данный момент я не добавляю шума), я имею в виду собственный спектр транспонирования X, то есть eig (cov (X ')).

ОБНОВИТЬ

Как и предлагалось (см. Комментарии), я попробовал FastICA только на 2 сигналах. Результаты были довольно хорошими (см. Рис ниже). Используемая матрица смешивания была A = [0.75 0.25; 0.25 0.75]. Тем не менее, eigenspectrum по- 0.1657 0.7732прежнему показал только один основной основной компонент.

Поэтому мой вопрос сводится к следующему: какую функцию / уравнение / свойство я могу использовать, чтобы проверить, является ли число векторов сигнала статистически независимым?

Синус & Гаусс - FastICA

Рейчел
источник
1
Отличный вопрос. Я спросил о том, как мы можем узнать, когда два сигнала здесь независимы ( dsp.stackexchange.com/questions/1242/… ), но не зашел слишком далеко. :-) Я также новичок в ICA, но я мог бы пролить свет.
Спейси
@ Мохаммад Вы все еще заинтересованы в ответе на этот вопрос? Я с радостью положу награду за это, чтобы привлечь интерес.
Фонон
@ Мохаммед Я проголосовал за твой вопрос. Надеюсь, вы получите хороший ответ, это действительно связано с моим. До сих пор я читал комментарии к нему, и сейчас идет много статистики, которую я не понимаю. Удалось ли вам придумать определенный способ сделать вывод, являются ли два сигнала независимыми или нет?
Рэйчел
@Rachel Пока что нет, но я еще немного разберусь и сообщу. Это очень важная концепция, которая, к сожалению, обычно застекляется.
Спейси
Спасибо, Мохаммед. Я согласен. Независимые сигналы соблюдают свойство E (s1, s2) = E (s1) x E (s2), но я не знаю, как на самом деле рассчитать его для реальных сигналов.
Рэйчел

Ответы:

8

Сигналы 3 и 5 выглядят достаточно коррелированными - они имеют свою первую гармонику. Если бы мне дали две их смеси, я бы не смог их разделить, я бы соблазнил поставить общую гармонику в качестве одного сигнала и высшую гармонику в качестве второго сигнала. И я был бы неправ! Это может объяснить отсутствующее собственное значение.

Сигналы 1 и 2 тоже не выглядят независимыми.

Быстрая и грязная «проверка работоспособности» для независимости двух рядов состоит в том, чтобы сделать (x, y) график одного сигнала против другого:

plot (sig3, sig5)

и затем сделать тот же (x, y) график с одним перетасованным сигналом:

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

Если два графика имеют разный вид, ваши сигналы не являются независимыми. В более общем смысле, если график (x, y) данных показывает «особенности», несимметрии и т. Д., Это плохое предзнаменование.

Надлежащие тесты на независимость (и это целевые функции, используемые в цикле оптимизации ICA) включают, например, взаимную информацию.

ICA восстанавливает большинство независимых сигналов, линейное микширование которых дает ваши входные данные . Он будет работать как методы разделения сигналов и восстанавливать исходные сигналы, только если они были максимально независимыми в соответствии с критерием оптимизации, используемым в вашей реализации ICA.

pichenettes
источник
1
Вопрос: Если бы 5 сигналов в ее случае были на самом деле все независимы, то мы ожидали бы, что НЕТ главных компонентов правильно? (Другими словами, все собственные значения будут одинаковыми). Геометрически, у нас было бы гусианское «облако» в 5 измерениях, согласны?
Спейси
Я также связался с автором по ICA по поводу удаления двух синусоид из смеси, и он сказал, что на самом деле это можно сделать с помощью ICA. Это немного смущает меня, основываясь на том, что вы говорите в отношении сигналов 3 и 5, потому что (я согласен с вами) они выглядят коррелированно.
Спейси
@pichenettes Я построил эти графики, как вы и предлагали, и они действительно выглядят по-разному. К сожалению, я все еще застрял в том, как проверить независимость. Мне действительно нужен способ генерировать сигналы, которые являются статистически независимыми, чтобы я мог оценить производительность FastICA.
Рэйчел
Икс1[N]Икс2[N]
@ Мохаммед Я не записывал свой собственный голос, но я пытался использовать FastICA на смеси синусодиального и гауссовского сигналов. Я бы подумал, что они независимы. FastICA выступил довольно хорошо, но собственный спектакль все еще был странным. Я обновлю свой вопрос, чтобы показать результаты.
Рэйчел
7

Я не эксперт по ICA, но могу немного рассказать о независимости.

Как упоминалось в некоторых комментариях, статистическую независимость между двумя случайными переменными можно приблизительно интерпретировать как «количество информации, которую дает одна переменная о другой».

ИксYИксYп(Икс,Y)ИксYп(Икс,Y)знак равноп(Икс)п(Y)

п(Икс,Y)

ИксYИксYп(Иксзнак равноя,Yзнак равноJ)знак равнопяJп(Иксзнак равноя)знак равнопяп(Yзнак равноJ)знак равнопJ

я(Икс,Y)знак равноΣяΣJпяJжурналпяJпяпJ

Вот код Matlab, который будет генерировать два независимых сигнала из построенного совместного распределения и два из независимого совместного распределения, а затем вычислять взаимную информацию о соединениях.

Функция «computeMIplugin.m» - это простая функция, которую я написал, которая вычисляет взаимную информацию, используя формулу суммирования выше.

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

Опять же, это предполагает, что у вас есть хорошая оценка совместного распределения (вместе с другими легкими предположениями), но это должно быть полезно, как правило.

Ага
источник
Это хороший ответ sydeulissie спасибо, я должен изучить это немного глубже.
Спейси
Прежде всего, спасибо за длинный ответ, он был очень информативным. У меня просто пара вопросов. Вы упомянули об использовании теста хи-квадрат. Я посмотрел на это и выглядит действительно интересно, но как я могу использовать его на сигналах? Разве это не может быть применено только к категориальным данным?
Рэйчел
Кроме того, вы используете Pj1 = P1 '* P2 для вычисления совместного распределения, правильно? Но технически я считаю, что это невозможно сделать. Может быть, вы делаете это, потому что вы предполагаете, что исходные сигналы независимы, и результат, таким образом, сохраняется? Но как вы можете рассчитать взаимную информацию - так как ее результат зависит от совместного распространения ..? Может быть, я что-то неправильно понял, но я хотел бы получить разъяснения, пожалуйста.
Рэйчел
Я буду счастлив - хотя это будет немного прежде, чем я получу время :).
Да,
Спасибо @sydeulissie. Я хотел бы получить ответ, который не предполагает, что я знаю о совместном распространении, пожалуйста.
Рейчел
3

Как упомянуто выше, оба сигнала 3 и 5 выглядят достаточно коррелированными и имеют одинаковый период.

Мы можем думать о том, что два сигнала коррелируют, если мы можем сместить один из источников влево или вправо и увеличить или уменьшить его амплитуду так, чтобы он подходил поверх другого источника. Обратите внимание, что мы не меняем частоту источника, мы просто выполняем сдвиг фазы и амплитуды.

В приведенном выше случае мы можем сместить источник 3 так, чтобы его пики совпадали с источником 5. Это такая вещь, которая будет мешать извлечению источника при использовании ICA из-за предположения о независимости.

Примечание : хорошая иллюстрация вышеупомянутой концепции - это думать о двух синусоидальных волнах. Они оба полностью детерминированы. Если они оба имеют одинаковую частоту (даже с разной фазой), то они идеально коррелируют и ICA не сможет их разделить. Если вместо этого они имеют разные частоты (которые не являются целыми числами, кратными друг другу), то они независимы и могут быть разделены.

Ниже приведен код Matlab, чтобы вы могли убедиться в этом сами.

%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')

Обратите внимание, что для волн одинаковой частоты ICA просто возвращает входные сигналы, но для разных частот возвращает исходные источники.

rwolst
источник
2

Рэйчел,

Из моих исследований я до сих пор смог найти что-то под названием « Тест хи-квадрат на независимость », но я не уверен, как он работает в данный момент, но, возможно, стоит посмотреть.

ошалевший
источник
Я нашел эти два учебника, объясняющих, как выполнить тест хи-квадрат: ling.upenn.edu/~clight/chisquared.htm & math.hws.edu/javamath/ryan/ChiSquare.html . Тем не менее, тест может быть выполнен только на категориальных данных. Я не знаю, может ли это быть применено к нашим наблюдениям сигналов ...
Рэйчел