ICA - Independencia estadística y valores propios de la matriz de covarianza

14

Actualmente estoy creando diferentes señales usando Matlab, mezclándolas multiplicándolas por una matriz de mezcla A, y luego tratando de recuperar las señales originales usando FastICA .

Hasta ahora, las señales recuperadas son realmente malas en comparación con las originales, que no era lo que esperaba.

Estoy tratando de ver si estoy haciendo algo mal. Las señales que estoy generando son las siguientes:

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

Señales originales

Una condición para que ICA tenga éxito es que, como máximo, una señal es gaussiana, y he observado esto en mi generación de señal.

Sin embargo, otra condición es que todas las señales son estadísticamente independientes.

Todo lo que sé es que esto significa que, dadas dos señales A y B, conocer una señal no proporciona ninguna información con respecto a la otra, es decir: P (A | B) = P (A) donde P es la probabilidad .

Ahora mi pregunta es esta: ¿son mis señales estadísticamente independientes? ¿Hay alguna forma de determinar esto? Quizás alguna propiedad que debe ser observada?

Otra cosa que he notado es que cuando calculo los valores propios de la matriz de covarianza (calculada para la matriz que contiene las señales mixtas), el espectro propio parece mostrar que solo hay un componente principal (principal) . ¿Qué significa esto realmente? ¿No debería haber 5, ya que tengo 5 (supuestamente) señales independientes?

Por ejemplo, cuando se usa la siguiente matriz de mezcla:

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

Los valores propios son: 0.0000 0.0005 0.0022 0.0042 0.0345(¡solo 4!)

Cuando se utiliza la matriz de identidad como la matriz de mezcla (es decir, las señales mixtas son los mismos que los originales), el eigenspectrum es: 0.0103 0.0199 0.0330 0.0811 0.1762. Todavía hay un valor mucho más grande que el resto.

Gracias por tu ayuda.

Pido disculpas si las respuestas a mis preguntas son dolorosamente obvias, pero realmente soy nuevo en estadísticas, ICA y Matlab. Gracias de nuevo.

EDITAR

Tengo 500 muestras de cada señal, en el rango [0.2, 100], en pasos de 0.2, es decir, x = 0: 0.1: 100.

Además, dado el modelo ICA: X = As + n (no estoy agregando ningún ruido en este momento), me estoy refiriendo al espectro propio de la transposición de X, es decir, eig (cov (X ')).

ACTUALIZAR

Como se sugiere (consulte los comentarios), probé FastICA con solo 2 señales. Los resultados fueron bastante buenos (ver foto abajo). La matriz de mezcla utilizada fue A = [0.75 0.25; 0.25 0.75]. Sin embargo, el espectro propio 0.1657 0.7732todavía mostraba solo un componente principal principal.

Por lo tanto, mi pregunta se reduce a lo siguiente: ¿Qué función / ecuación / propiedad puedo usar para verificar si varios vectores de señal son estadísticamente independientes?

Seno y gaussiano - FastICA

Rachel
fuente
1
Excelente pregunta Pregunté cómo podemos saber cuándo dos señales son independientes aquí ( dsp.stackexchange.com/questions/1242/… ) pero no llegué demasiado lejos con eso. :-) También soy nuevo en ICA pero podría arrojar algo de luz.
Spacey
@Mohammad ¿Todavía estás interesado en una respuesta a esa pregunta? Con mucho gusto pondré una recompensa para atraer el interés.
Phonon
@Mohammad I votó por tu pregunta. Espero que obtengas una buena respuesta, está realmente relacionada con la mía. He estado leyendo los comentarios hasta el momento, y hay muchas estadísticas que no entiendo. ¿Has logrado encontrar una forma definitiva de concluir si dos señales son independientes o no?
Rachel
@Rachel Todavía no por el momento, pero lo investigaré un poco más y te lo haré saber. Es un concepto muy importante, uno que siento que, por desgracia, suele estar satinado.
Spacey
Gracias @ Mohammad. Estoy de acuerdo. Las señales independientes observan la propiedad de que E (s1, s2) = E (s1) x E (s2), pero no sé cómo calcularlo para señales reales.
Rachel

Respuestas:

8

Las señales 3 y 5 parecen estar bastante correlacionadas: comparten su primer armónico. Si me dieran dos mezclas de esas, no podría separarlas, estaría tentado a poner el armónico común como una señal y los armónicos más altos como una segunda señal. ¡Y estaría equivocado! Esto podría explicar el valor propio que falta.

Las señales 1 y 2 tampoco se ven independientes.

Un "control de cordura" rápido y sucio para la independencia de dos series es hacer un diagrama (x, y) de una señal contra la otra:

plot (sig3, sig5)

y luego hacer el mismo diagrama (x, y) con una señal mezclada:

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

Si las dos parcelas tienen un aspecto diferente, sus señales no son independientes. En términos más generales, si la gráfica (x, y) de los datos muestra "características", disimetrías, etc., es un mal presagio.

Las pruebas adecuadas de independencia (y esas son las funciones objetivas utilizadas en el ciclo de optimización ICA) incluyen, por ejemplo, información mutua.

ICA se está recuperando las señales más independientes, una mezcla lineal de las cuales produce sus datos de entrada . Funcionará como un método de separación de señal y recuperará las señales originales solo si esas eran máximamente independientes de acuerdo con el criterio de optimización utilizado en su implementación de ICA.

pichenettes
fuente
1
Pregunta: Si las 5 señales en su caso fueran de hecho todas independientes, entonces esperaríamos que NO hubiera componentes principales correctos. (En otras palabras, todos los valores propios serían más o menos lo mismo). Geométricamente, tendríamos una 'nube' guasiana en 5 dimensiones, ¿de acuerdo?
Spacey
También me puse en contacto con un autor en ICA sobre la eliminación de dos sinusoides de una mezcla, y dijo que, de hecho, podría hacerse con ICA. Esto me confunde un poco en función de lo que está diciendo con respecto a las señales 3 y 5 porque (estoy de acuerdo con usted) parecen correlacionadas.
Spacey
@pichenettes Tracé esas gráficas como usted sugirió, y las parcelas tienen un aspecto diferente. Lamentablemente, todavía estoy atascado en cuanto a cómo probar la independencia. Realmente necesito una forma de generar señales que sean estadísticamente independientes, para poder evaluar el rendimiento de FastICA.
Rachel
x1[n]x2[n]
@Mohammad No grabé mi propia voz, pero intenté usar FastICA en una mezcla de señales sinusodiales y gaussianas. Creo que son independientes. FastICA funcionó bastante bien, pero el espectro propio aún era extraño. Actualizaré mi pregunta para mostrar los resultados.
Rachel
7

No soy un experto en ICA, pero puedo contarles un poco sobre la independencia.

Como algunos de los comentarios han mencionado, la independencia estadística entre dos variables aleatorias puede interpretarse aproximadamente como "la cantidad de información que la observación de una variable proporciona sobre otra".

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

Aquí hay un código matlab que generará dos señales independientes de una distribución conjunta construida, y dos de una distribución conjunta no independiente, y luego calculará la información mutua de las juntas.

La función "computeMIplugin.m" es una función simple que escribí que calcula la información mutua usando la fórmula de suma anterior.

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

Nuevamente, esto supone que tiene una buena estimación de la distribución conjunta (junto con otros supuestos desenfadados), pero debería ser útil como regla general.


fuente
Esa es una buena respuesta sydeulissie gracias, Voy a tener que investigar un poco más profundo.
Spacey
En primer lugar, gracias por la larga respuesta, fue muy informativo. Sólo tengo un par de preguntas. Usted mencionó el uso de una prueba de chi-cuadrado. Lo he investigado y parece realmente interesante, pero ¿cómo puedo usarlo en las señales? ¿No se puede aplicar solo a datos categóricos?
Rachel
Además, está utilizando Pj1 = P1 '* P2 para calcular la distribución conjunta, ¿correcto? Pero, técnicamente, creo que esto no se puede hacer. ¿Tal vez lo estás haciendo porque estás asumiendo que las señales originales son independientes y el resultado se mantiene? Pero, ¿cómo puede calcular la información mutua, ya que su resultado depende de la distribución conjunta ...? Puede ser que haya entendido mal algo, pero me gustaría una aclaración, por favor.
Rachel
Estaré encantado de hacerlo, aunque será un poco antes de que tenga tiempo :).
Gracias @sydeulissie. Quisiera una respuesta que no suponga que tengo conocimiento de la distribución conjunta, por favor.
Rachel
3

Como se mencionó anteriormente, ambas señales 3 y 5 parecen estar bastante correlacionadas y tienen un período similar.

Podemos pensar en dos señales correlacionadas si podemos desplazar una de las fuentes hacia la izquierda o hacia la derecha y aumentar o disminuir su amplitud para que se ajuste a la otra fuente. Tenga en cuenta que no estamos cambiando la frecuencia de la fuente, solo estamos realizando un cambio de fase y amplitud.

En el caso anterior, podemos cambiar la fuente 3 para que sus picos coincidan con la fuente 5. Este es el tipo de cosa que afectará la extracción de la fuente cuando se usa ICA debido al supuesto de independencia.

Nota : Una buena ilustración del concepto anterior es pensar en dos ondas sinusoidales. Ambos son completamente deterministas. Si ambos tienen la misma frecuencia (incluso con una fase diferente), entonces están perfectamente correlacionados e ICA no podrá separarlos. Si, en cambio, tienen diferentes frecuencias (que no son múltiplos enteros entre sí), entonces son independientes y pueden separarse.

A continuación se muestra un código de Matlab para que pueda verlo usted mismo.

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

Tenga en cuenta que para las ondas de la misma frecuencia, ICA solo devuelve las señales de entrada, pero para diferentes frecuencias devuelve las fuentes originales.

rwolst
fuente
2

Rachel

De mi investigación, hasta ahora he podido encontrar algo llamado ' Prueba de independencia de Chi-cuadrado ', pero no estoy seguro de cómo funciona en este momento, pero podría valer la pena echarle un vistazo.

Spacey
fuente
Encontré estos dos tutoriales que explican cómo realizar la prueba de chi-cuadrado: ling.upenn.edu/~clight/chisquared.htm y math.hws.edu/javamath/ryan/ChiSquare.html . Sin embargo, la prueba solo se puede realizar con datos categóricos. No sé si esto se puede aplicar a nuestras observaciones de señal ...
Rachel