独立成分分析（Independent component analysis）

一 算法流程：

其中α为梯度上升速率，人为制定。

二 前置知识

主要是概率论、线性代数的知识。

设随机变量s有概率密度函数p_s(s) (连续值是概率密度函数，离散值是概率！)。我们假设s是实数，还有一个随机变量x = As，A和x都是实数。【这里的s相当于我们的原始信号，A是一个随机组合的矩阵，x即为我们观测到的信号。
令p_x是x的概率密度。现在我们要求p_x，令W=A^-1，首先将式子变换成s=Wx，然后得到p_x(x)=p_s(Ws)|W|。
推导方式：
F_x(x)=P(X≤x)=P（AS≤x）=P（S≤Wx）=F_s(Wx)
p_x(x)=F_x'(x)=F_s'(Wx)=p_s(Wx)|W|

对行列式求导，也就求行列式的梯度。设矩阵A是nxn的，我们知道行列式与代数余子式有关，!

三 算法理论推导

使用的是最大似然估计法来解释的。（突然联想EM算法）

我们假设每个s_i有概率密度p_s，那么给定时刻原信号的联合分布就是

这个公式代表一个前提：每个发声源都是各自独立的。有了p（s），我们就可以求得p（x） 左边是每个采样信号x（n维向量）的概率，右边是每个原信号概率的乘积|W|倍。

我们必须选取一个概率密度函数F（x）赋给s，才能继续求解W和s。F（x）要满足两个性质：单调递增和在[0,1]区间内，同时不能该是高斯分布的概率密度。于是我们发现了sigmoid函数很合适，定义域为负无穷到正无穷，值域为0到1，同时还是一个对称函数（即期望/均值为0）。

其导数为： 【注1】表达式中的s要取负号 限制我们就能求W了，在给定训练样本{x⁽ⁱ⁾(x⁽ⁱ⁾₁,x⁽ⁱ⁾₂,x⁽ⁱ⁾₃,...,,x⁽ⁱ⁾_n);i = 1,...,m}后，对样本的对数似然估计如下：

使用前面的概率密度函数p（x），得到

其中大的小括号即为p（x）取对数 其中的导数涉及到了对行列式|W|的求导，属于矩阵微积分，在上述已经解释过了。 ，最后求导之后的公式如下，【注2】其中我计算出来log(g'(s))的导数是-1-2g(s)，参考博客是1-2g(s). 其中α是梯度上升速率，人为指定。
当迭代求出W后，便可得到s⁽ⁱ⁾=Wx⁽ⁱ⁾求出原始信号。

括号中前一部分为直接求导后的形式，后一部分|W|约分后，求转置与求逆交换后的结果。

【没有代码，其实这里我对迭代结束判定条件还存有疑惑。

四 算法的前处理步骤

我本来以为前置处理在算法流程之中，在此本来准备写理解的，不过看来得写在一起了。

1.中心化：理同PCA，使数据处于坐标中心。
2.白化/球化：
白化是一种重要的预处理过程，其目的是为了降低输入数据的冗余性，使得经过白化处理的输入数据（i）特征之间相关性较低 （ii）所有特征具有相同的方差。
实际上，如果我们不用PCA降维
①仅仅用PCA求出特征向量，然后把数据映射到新的特征空间，就能够满足我们的的第一个性质
②对新的坐标做方差归一化操作，就满足第二个性质，这就叫PCA白化。
③如果我们PCA白化后的结果，变换回原来的坐标系，那么就叫做ZCA白化。

def zca_whitening(inputs):  
    sigma = np.dot(inputs, inputs.T)/inputs.shape[1] #inputs是经过中心化处理的，所以这边就相当于计算协方差矩阵  
    U,S,V = np.linalg.svd(sigma) #奇异分解  
    epsilon = 0.1                #白化的时候，防止除数为0  
    ZCAMatrix = np.dot(np.dot(U, np.diag(1.0/np.sqrt(np.diag(S) + epsilon))), U.T) #计算zca白化矩阵  
    return np.dot(ZCAMatrix, inputs)   #白化变换  

参考：
独立成分分析（Independent Component Analysis）
PCA 白化 ZCA白化
独立成分分析ICA系列2：概念、应用和估计原理．

五 FastICA/固定点（Fixed-Point）算法

FastICA有基于峭度、基于似然最大、基于负熵最大等形式，这里介绍的是基于负熵最大的FastICA算法。

由信息论理论：在所有等方差的随机变量中，高斯变量的熵最大，因而我们可以利用熵来度量非高斯性，常用熵的修正形式，即负熵。根据中心极限定理，若一随机变量X由许多相互独立的随机变量S_i（i=1,2,...,N）之和组成，只要S_i具有有限的均值和方差，则不论其为何种分布，随机变量X较S_i更接近高斯分布。因此，在分离过程汇总，可通过对分离结果的非高斯性度量来表示分离结果间的互相独立性，当非高斯性度量达到最大时，则表明已完成对各个独立分量的分离。

负熵的定义：

上式中其中Y_Gauss是与Y具有相同方差的高斯随机变量。H（-）为随机变量的微分熵，定义如下： 根据信息理论，在具有相同方差的随机变量中，高斯分布的随机变量具有最大的微分熵。当Y具有高斯分布时，N_g=0；Y的非高斯性越强，其微分熵越小，N_g（Y）的值越大，所以N_g（Y）可以作为随机变量Y的非高斯性的测度。采用负熵定义求解需要知道Y的概率分布函数，但是实际不可能，于是采用近似公式： 其中E[-]为均值运算，g(.)为非线性函数，可取g1(y)=tanh(a1y)或g2(y)=y exp(-y2/2)或g3(y)=y3等非线性函数，这里1≤a1≤2，通常取a1=1

快速ICA的规则就是找到一个方向以便W^TX（Y=W^TX）具有最大的非高斯性，非高斯性用Ng(Y)={E[g(Y)]-E[g(Y_Gauss)]}²给出的负熵的近似值来度量。WTX的方差约束为1，对于白化数据，等于约束W的范数为1.

FastICA的推导：

①我们知道W^TX是我们还原的函数，所以目标是时W^TX的负熵最大，又知最大近似值能通过对E{G(W^TX)}进行优化取得。根据Kuhn-Tucker条件，在E{( W^TX)²}=||W||²=1的约束下，E{G(W^TX)}的最优值能在满足下式的点上获得。

E{Xg(W^TX)}+βW=0

其中β=E{W₀^TXg(W^TX)} 是一个恒定值，W₀是优化后的W值。

②利用牛顿迭代法解①的方程。

用F表示左边的函数，得到F的雅克比矩阵JF(W)如下：

为了简化可以近似为第一项，即忽略βI。

由于数据被球化，所以E{XX^T}=I,所以E{XX^Tg’(W^TX)}≈E{XX^T}*E{g’(W^TX)}=E{g’(W^TX)}I。

从而雅克比矩阵变成了对角阵，并且比较容易求逆。因而得到下面的近似牛顿迭代公式：

这里的W*是W的新值，β= E{W^TXg(W^TX)}，规格化能提高稳定性。

实践中，FastICA算法中用的期望必须用他们的估计值代替。最好的估计是相应的样本平均。理想情况下，所有的有效数据都应该参与计算，但是会降低运算速度，所以通常选取一部分样本的平均来估计，样本数目的多少对最后估计的精确度有很大影响。迭代中的样本点应该分别选取，加入收敛不理想，可以增加样本数量。

简化后得到FastICA的迭代公式：

%FastICA的matlab代码
function Z = ICA( X )
 
%去均值
[M,T]=size(X);   %获取输入矩阵的行列数，行数为观测数据的数目，列数为采样点数
average=mean(X')';  %均值
for i=1:M
    X(i,:)=X(i,:)-average(i)*ones(1,T);
end
 
%白化/球化
Cx=cov(X',1);  %计算协方差矩阵Cx
[eigvector,eigvalue]=eig(Cx);   %计算Cx的特征值和特征向量
W=eigvalue^(-1/2)*eigvector';   %白化矩阵
Z=W*X;    %正交矩阵
 
%迭代
Maxcount=10000;  %最大迭代次数
Critical=0.00001;  %判断是否收敛
m=M;
W=rand(m);
for n=1:m
    WP=W(:,n);  %初始权矢量（任意）
    %Y=WP'*Z;
    %G=Y.^3;%G为非线性函数，可取y^3等
    %GG=3*Y.^2；   %G的导数
    count=0;
    LastWP=zeros(m,1);
    W(:,n)=W(:,n)/norm(W(:,n));
    while abs(WP-LastWP)&abs(WP+LastWP)>Critical
        count=count+1;  %迭代次数
        LastWP=WP;      %上次迭代的值
        %WP=1/T*Z*((LastWP'*Z).^3)'-3*LastWP;
        for i=1:m
            WP(i)=mean(Z(i,:).*(tanh((LastWP)'*Z)))-(mean(1-(tanh((LastWP))'*Z).^2)).*LastWP(i);
        end
        WPP=zeros(m,1);
        for j=1:n-1
            WPP=WPP+(WP'*W(:,j))*W(:,j);
        end
        WP=WP-WPP;
        WP=WP/(norm(WP));
        if count==Maxcount
            fprintf('未找到相应的信号');
            return;
        end
    end
    W(:,n)=WP;
end
Z=W'*Z;
end

补充：在做胎心电项目的时候，发现ICA对信号源非常敏感，所以在实际操作中应该尽量消除原始信号的各种噪音。
例如：在MIT-BIH 2013年比赛【Noninvasive Fetal ECG】中，我们取a04数据，直接作图是这样的：

原始信号 其实可以看得出来，这四路信号还是非常规律的，但是其中有不少基线漂移，环境噪音（我仅通过毛刺来看的。。）所以直接运用ICA的结果是： 对原始信号直接做ICA 做完之后反而更看不到规律了，这是因为我们用默认是只有四个源信号，导致噪音无法归类/排除。

在 A Multi-step Approach for Non-invasive Fetal ECG Analysis 中，对原始信号做了各种除噪音处理后，得到的图更加好看：

除噪音后原始数据
具体使用的方法是：
1.For each ECG channel a median filtering with a 60ms window was applied.
2.The baseline signal was computed applying alow pass first order Butterwort filter in forwad and backward direction to avoid phase distortion(cut off frequence at 3.17Hz)
3.a notch filter(forward-backward,zero phase,1Hz bandwidth) was applied to remove its characteristic frequency and its next three harmonics 对这个数据进行ICA后： 对除噪音数据进行ICA

一些高斯噪音很明显分出来了，最重要的是母亲和胎儿的心电信号提取出来了！

我自己做了一个简单的信号模拟来测试ICA的具体性能。
结果表明ICA会将噪音提取出来，并且将噪音作为一个信号源，这样你的输入路数最好是等于信号源数+噪音数，对于比较好的信号源分离效果非常好，反之，如果噪音过多、输入路数过多等，分离效果尤其地差！

Fs1=8;
Fs2=5;
T1=1/Fs1;
T2=1/Fs2;
x=(1:1000);
y1=sin(x*T1*2/pi);
y2=2*cos(x*T2*2/pi);
x=rand(1,1000);
% 两个信号源，四路信号，结果有问题
% data=[y1+2*y2+x;0.5*y1+1.3*y2+x;0.7*y1+1.2*y2+2*x;1.7*y1+5.9*y2+x];
% data=[y1+2*y2+x;0.5*y1+1.3*y2+x;0.7*y1+1.2*y2+2*x;1.4*y1+2.1*y2+x];
% 两个信号源，有噪音，三路信号,能正确得到结果
% data=[y1+2*y2+x;0.5*y1+1.3*y2+x;0.7*y1+1.2*y2+2*x];
% 两个信号源，无噪音，三路信号,不能正确得到结果
data=[y1+2*y2;0.5*y1+1.3*y2;0.7*y1+1.2*y2];
% 两个信号源，有噪音，两路信号,不能得到正确结果
% data=[y1+2*y2+x;0.5*y1+1.3*y2+x];
% 两个信号源，无噪音，两路信号,能得到正确结果
% data=[y1+2*y2;0.5*y1+5.4*y2];
Z=ICA(data);
subplot(3,1,1)
plot(Z(1,:));
subplot(3,1,2)
plot(Z(2,:))
subplot(3,1,3)
plot(Z(3,:))

% 4
subplot(4,1,1)
plot(Z(1,:));
subplot(4,1,2)
plot(Z(2,:))
subplot(4,1,3)
plot(Z(3,:))
subplot(4,1,4)
plot(Z(4,:))

补1：中心极限定理

中心极限定理是说：不论原来的总体是否服从正态分布，样本均值的抽样分布都将趋于正态分布。

更多的细节我希望在重学概率论与数理统计之后再开一个基础知识方面的专题。

补2：投影追踪（Projection Pursuit）线性变换

本方法是摘自《独立分量分析原理与应用》-杨福生，因为网上无法免费阅读，所以在我借到学长的书之后看看再填坑。