1.

ICA算法扩展描述

《Independent Component Analysis:Algorithms and Applications》(Aapo

Hyvärinen and Erkki Oja)文章里提到的一些内容(有些我也没看明白)。

首先里面提到了一个与“独立”相似的概念“不相关(uncorrelated)”。Uncorrelated属于部分独立，而不是完全独立，怎么刻画呢？

如果随机变量和是独立的，当且仅当。

如果随机变量和是不相关的，当且仅当

第二个不相关的条件要比第一个独立的条件“松”一些。因为独立能推出不相关，不相关推不出独立。

证明如下：

反过来不能推出。

比如，和的联合分布如下(0,1)，(0,-1)，(1,0)，(-1,0)。

因此和不相关，但是

因此和不满足上面的积分公式，和不是独立的。

上面提到过，如果是高斯分布的，A是正交的，那么也是高斯分布的，且与之间是独立的。那么无法确定A，因为任何正交变换都可以让达到同分布的效果。但是如果中只有一个分量是高斯分布的，仍然可以使用ICA。

那么ICA要解决的问题变为：如何从x中推出s，使得s最不可能满足高斯分布？

中心极限定理告诉我们：大量独立同分布随机变量之和满足高斯分布。

我们一直假设的是是由独立同分布的主元经过混合矩阵A生成。那么为了求，我们需要计算的每个分量。定义，那么，之所以这么麻烦再定义z是想说明一个关系，我们想通过整出一个来对进行线性组合，得出y。而我们不知道得出的y是否是真正的s的分量，但我们知道y是s的真正分量的线性组合。由于我们不能使s的分量成为高斯分布，因此我们的目标求是让y(也就是)最不可能是高斯分布时的w。

那么问题递归到如何度量y是否是高斯分布的了。

一种度量方法是kurtosis方法，公式如下：

如果y是高斯分布，那么该函数值为0，否则绝大多数情况下值不为0。

但这种度量方法不怎么好，有很多问题。看下一种方法：

负熵(Negentropy)度量方法。

我们在信息论里面知道对于离散的随机变量Y，其熵是

连续值时是

在信息论里有一个强有力的结论是：高斯分布的随机变量是同方差分布中熵最大的。也就是说对于一个随机变量来说，满足高斯分布时，最随机。

定义负熵的计算公式如下：

也就是随机变量y相对于高斯分布时的熵差，这个公式的问题就是直接计算时较为复杂，一般采用逼近策略。

这种逼近策略不够好，作者提出了基于最大熵的更优的公式：

之后的FastICA就基于这个公式。

另外一种度量方法是最小互信息方法：

这个公式可以这样解释，前一个H是的编码长度(以信息编码的方式理解)，第二个H是y成为随机变量时的平均编码长度。之后的内容包括FastICA就不再介绍了，我也没看懂。

8.

ICA的投影追踪解释(Projection Pursuit)

投影追踪在统计学中的意思是去寻找多维数据的“interesting”投影。这些投影可用在数据可视化、密度估计和回归中。比如在一维的投影追踪中，我们寻找一条直线，使得所有的数据点投影到直线上后，能够反映出数据的分布。然而我们最不想要的是高斯分布，最不像高斯分布的数据点最interesting。这个与我们的ICA思想是一直的，寻找独立的最不可能是高斯分布的s。

在下图中，主元是纵轴，拥有最大的方差，但最interesting的是横轴，因为它可以将两个类分开(信号分离)。

9.

ICA算法的前处理步骤

1、中心化：也就是求x均值，然后让所有x减去均值，这一步与PCA一致。

2、漂白：目的是将x乘以一个矩阵变成，使得的协方差矩阵是。解释一下吧，原始的向量是x。转换后的是。

的协方差矩阵是，即

我们只需用下面的变换，就可以从x得到想要的。

其中使用特征值分解来得到E(特征向量矩阵)和D(特征值对角矩阵)，计算公式为

下面用个图来直观描述一下：

假设信号源s1和s2是独立的，比如下图横轴是s1，纵轴是s2，根据s1得不到s2。

我们只知道他们合成后的信号x，如下

此时x1和x2不是独立的(比如看最上面的尖角，知道了x1就知道了x2)。那么直接代入我们之前的极大似然概率估计会有问题，因为我们假定x是独立的。

因此，漂白这一步为了让x独立。漂白结果如下：

可以看到数据变成了方阵，在的维度上已经达到了独立。

然而这时x分布很好的情况下能够这样转换，当有噪音时怎么办呢？可以先使用前面提到的PCA方法来对数据进行降维，滤去噪声信号，得到k维的正交向量，然后再使用ICA。

10.

小结

ICA的盲信号分析领域的一个强有力方法，也是求非高斯分布数据隐含因子的方法。从之前我们熟悉的样本-特征角度看，我们使用ICA的前提条件是，认为样本数据由独立非高斯分布的隐含因子产生，隐含因子个数等于特征数，我们要求的是隐含因子。

而PCA认为特征是由k个正交的特征(也可看作是隐含因子)生成的，我们要求的是数据在新特征上的投影。同是因子分析，一个用来更适合用来还原信号(因为信号比较有规律，经常不是高斯分布的)，一个更适合用来降维(用那么多特征干嘛，k个正交的即可)。有时候也需要组合两者一起使用。这段是我的个人理解，仅供参考。

% 这是一个ICA人脸识别算法实例，百度文库中的,作为第一次学习ICA算法

%---初始化--------------------------------

clc;

clear all;

close all;

% -----读入原始图像，混合，并输出混合图像---------------------------

% 读入混合前的袁术图像图片并显示

t = 0:1/100:9

T1 = sin(t);

T2 = randn(1,901);

T3 = square(4*t);

subplot(4,3,1),plot(T1);title('输入信号1');

subplot(4,3,2),plot(T2);title('输入信号2');

subplot(4,3,3),plot(T3);title('输入信号3');

%将其组成矩阵

S = [T1;T2;T3];% 图片个数即为变量数，图片的像素数即为采样数

% 因此 S_all 是一个变量个数 * 采样个数的矩阵

a = size(S,1); % 得到S 的行数

Sweight = randn(size(S,1)); % 只有一个数的时候，默认产生 一个方阵

% 上边是 去一个随机矩阵，作为信号混合的权矩阵

MixedS = Sweight * S; % 得到三个混合信号矩阵

% 将 混合矩阵重新排列并输出

subplot(4,3,4),plot(MixedS(1,:)),title('混合信号1');

subplot(4,3,5),plot(MixedS(2,:)),title('混合信号2');

subplot(4,3,6),plot(MixedS(3,:)),title('混合信号3');

MixedS_bak = MixedS; % 将混合后的数据备份，以便在恢复时直接调用

% -------标准化----------------------------------------

MixedS_mean = zeros(3,1);

for i = 1:3

MixedS_mean(i) = mean(MixedS(i,:));

end

% 上边是计算MixedS 的均值

for i = 1:3

for j = 1:

size(MixedS,2)

MixedS(i,j) = MixedS(i,j) - MixedS_mean(i);

end

end

%

---------白化----------------------------------------------

MixedS_cov = cov(MixedS'); % cov为求协方差的函数

[E,D] = eig(MixedS_cov); % 对图片矩阵的协方差函数进行特征值分解

Q = inv(sqrt(D))*(E)'; % Q为白化矩阵

MixedS_white = Q*MixedS; % MixedS_white为白化后的图片矩阵

IsI = cov(MixedS_white'); % IsI 应为单位阵

% ---------FASTICA算法---------------------------------

X = MixedS_white; % 这也是编程的小技巧，下边是对X进行操作， 不会改变MixedS_whiteD的值

[VariableNum,SampleNum] = size(X);

numofIC = VariableNum; %在此 应用中，独立元个数 等于变量的个数

B = zeros(numofIC,VariableNum);

for r = 1:numofIC % 迭代求取每一个独立元

i =

1; maxIterationsNum = 100; % 设置最大迭代次数(即对于每个独立分量而言迭代均不超过此次数)

IterationsNum = 0;

b =

rand(numofIC,1) - 0.5; % 随机设置b 初值

b =

b/norm(b); % 对b标准化

while i

<= maxIterationsNum + 1

if i == maxIterationsNum % 循环结束处理

fprintf('\n第%d分量在%d次迭代内并不收敛。',r,maxIterationsNum);

break;

end

bOld = b;

a2 = 1;

u = 1;

t = X'*b;

g = (exp(2.*t) - 1)./(exp(2.*t) +1);

dg = 4*exp(2.*t)./(exp(2.*t) + 1).^2;

b = ((1 -u)*t'*g*b + u*X*g)/ SampleNum - mean(dg)*b;

b = b - B*B'*b; % 对b进行正交化

b = b/norm(b);

if abs(abs(b'*bOld)-1) < 1e-9 % 如果收敛，则保存b

B(:,r) = b;

break;

end

i = i +1;

end

end

% --------------数据复原并构图---------------------------

ICAedS = B'*Q*MixedS_bak;

% 将混合矩阵重新排列并输出

subplot(437),plot(ICAedS(1,:)),title('ICA输出信号1');

subplot(438),plot(ICAedS(2,:)),title('ICA输出信号2');

subplot(439),plot(ICAedS(3,:)),title('ICA输出信号3');

%

----------------PCA计算并构图----------------------------------------------

[V,D] = eig(MixedS_cov);

Vtmp = zeros(size(V,1),1);

for j = 1:2

for i =

1:2

if D(i,i)