ICA

         本总结是是个人为防止遗忘而作，不得转载和商用。

      ICA的著名应用是盲源分离，于是这里就以盲源分离为例子进行说明。

 

题目

      假设n个人面前有n个话筒，然后这n个人说话时这n个话筒进行录音，这n个人说了m句话，最后从这n个话筒中收集一些录音，目标：从这些录音中分离出每个人的声音。

      如下图所示：

             

      下面开始解题。

 

题目整理

      首先将信息整理如下：

      源信号，如下图所示：

      PS：源信号是未知的（知道还求啥），但为了分析还是需要把源信号表示出来。

      第一个人是s₁，他第一个时刻发送的信号是s₁₁，第二个时刻发送的信号是s₁₂，以此类推第m个时刻发送的时刻是s_1m。第二个人，第三个人，….，第n个人同理。

      因此列向量中的每个元素就代表每个人，行向量的每个元素代表每个时刻。

 

      观测信号，如下图所示：

      解读方式和源信号一样。

      这里我们需要的是源信号和观测信号行向量，因为我们只知道每个时刻大伙说了什么话。于是s_m和x_m都是n维的向量。

      然后，我们假设源信号经过线性加权得到了观测信号， 即：s_m乘上一个权值A后得到了x_m，即：x_m = As_m。接下来为了表示所有的s_m和x_m，就用s代表s_m组成的矩阵，x同理，于是式子更新为x = A·s。

       到此，模型出来了，即：

              x =A·s

      后面要做的是如何求解该模型，但因为除了观测信号x之外，A未知(不知道怎么混合的)，s也未知（源信号不知道）。

      别忘了我们的目标：分离源信号。所以要把上面的式子转换成源信号=什么什么，即，转换成：

           s = A^-1·x，

      为了方便书写，用W表示A的逆矩阵(或广义逆矩阵)，于是上式变成：

                   s = W·x

          之后，我们的核心目标就是如何求W（x已知，s是目标）。

          到此题目整理完毕。

      PS：这里可以看出ICA的一个不足：顺序和振幅不稳定，即：对于ICA模型x = A·s，甲求得A=A₁，S=S₁；乙求得A=-A₁/2，S=-2S₁，那甲乙的结果都能得出源信号，但甲和乙求得的结果的顺序就是反的（有个负号），振幅甲是乙的1/2。

 

进行假设

      首先，先做两个假设：

              假设1：源信号彼此间独立。

              假设2：源信号是非高斯分布。

      上面的假设1不用解释，假设2是的根据是中心极限定理，即：一组有确定方差的独立随机变量的和趋近于高斯分布。

不过假设2可以多解释些，嗯….如果打个比方的话：假设源信号是一个个颜色时，高斯分布就是颜色混合后的那一坨脏色，任何颜色（源信号）一旦和其他颜色混合（成为高斯分布），那该颜色就再也无法被识别。换句话说：给定随机变量A和B时，A+B比A或B更接近高斯分布，再换句话说：如果能够找到一组独立信号，或者说找到了一组“最不像高斯分布”的信号，则它们极有可能是源信号。

 

开始推导ICA

推导的工具就是最大似然估计。

假定第i个源信号的概率密度函数为p_i(s)，第j时刻的n个源信号记做向量s_j，则在j时刻向量s_j 的联合密度为：

            

      根据x_j = A·s_j和http://blog.csdn.net/xueyingxue001/article/details/51915390中“复合函数的概率密度”的知识，得：

          

     

      从而得到似然函数：

          

      进一步得到对数似然：

            

      从而，目标函数为：

          

      但是，对于目标函数，除了我们相求的W外，我们还不知道每个人说话的概率密度，即：p_i(W_i^T·x_j)，这点比较郁闷，不过我们可以做一个有很强根据的假定：

           假定源信号的累积概率密度函数是Logistics/Sigmoid函数。

      为什么说这个假定有很强的依据？

      如下图所示：

          

      这是Sigmoid的函数图像，从图像中可以很清楚的看到，在负无穷时Sigmoid函数接近0，正无穷时接近1，所以，我们就可以将其看成是某个随机变量的累积概率。所以在没有任何先验信息时，我们不妨认为源信号的累积概率密度大体符合Sigmoid函数。（当然也可以假设成其他函数，最后在对比效果。）

      接下来给出Sigmoid函数的公式和求导（即：概率密度）：

           公式：F(x) = 1/(1+e^-x)

           概率密度/求导：f(x) = F’(x) =F(x)(1-F(x))

      为了之后的步骤方便，再求次导数：

           概率密度在求导：f’(x) = F’’(x) = f(x)(1-2F(x))

      那按照极大似然估计的知识，接下来我们就直接对J(W)对w_ij求偏导：

          

      PS：上面第二行利用了http://blog.csdn.net/xueyingxue001/article/details/51915390中“标量对矩阵求导”的知识

      既然对w_ij求偏导是上面的结果，那对整个W求偏导就是把上面的结果写成向量的形式，即：

                

      于是梯度下降的公式就有了，即：

                

           PS1：如果W不是方阵，则只要把W^-1替换成W⁺就可以了。

           PS2：α是学习率。

 

使用时：

      虽然教科书中说在使用ica时，对于一个数据要经过如下步骤：

           原始数据->中心化去均值->白化->PCA降维->ICA

      但实时上：

           若噪声不太强，PCA可以忽视

           对于白化，有些数据不做此步骤会更好。

 

最后，从矩阵的角度总结ICA：

      对于ICA的模型：x = A·s

       假设x是个m*n的矩阵，A是个m*k的矩阵，那s就需要是个k*n的矩阵，于是ICA就是把x分解成2个矩阵的乘积，这两个矩阵就是求出来的源信号。

     

代码：

#-*-coding:utf-8-*-import sysimport timeimport mathimport datetimeimport threadingimport numpy asnpimportmatplotlib.pyplot as plt def show_data(x,y):    num_list_x = np.arange(len(x))    plt.figure(figsize=(20, 10))    plt.xlim(0, len(num_list_x)+50)    plt.plot(num_list_x, x, color='b',linestyle='-', label='x')    plt.legend()    plt.show()    if y != None:        plt.figure(figsize=(20, 10))        num_list_y = np.arange(len(y))        plt.xlim(0, len(num_list_y)+50)        plt.plot(num_list_y, y, color='r',linestyle='-', label='y')        plt.legend()        plt.show() def logistic(t):    return 1.0/(1+np.exp(-t)) defn_multiply(t, x):    return t * x deftrans_inverse(x):    if not isinstance(x, np.ndarray):        x = np.array(x)    return np.linalg.inv(x.transpose()) def ica(x):    m = len(x)      # 样本数目，这里是 2 个    n = len(x[0])   # mic 数目，这里是 1000 个    # 建立个 1000*1000 的列表    w = [[0.0]*n for t in range(n)]    # 建立个 1000*1000 的列表    iw = [[0.0]*n for t in range(n)]    # 建立个 1000*1000 的列表    w1 = [[0.0]*n for t in range(n)]    # 将对角线初始化为1    for i in range(n):        w[i][i] = 1    # 初始化学习率为    alpha = 0.001    # shuffle(x)    # 迭代次数最多不超过 200 次    for time in range(200):        # 分离出"样本数目"个样本        for i in range(m):            for j in range(n):                t = np.dot(w[j], x[i])                t = 1 - 2*logistic(t)                w1[j] = n_multiply(t,x[i])  # w1[j] = t*x[i]            iw = trans_inverse(w)    # iw = w^T^(-1)            iw = w1 + iw            w1 = np.add(w1, iw) #w1 += iw            w1 = np.dot(w1, alpha)            w = np.add(w, w1)        print time, ":\t", w    return w def mix(A, x,y):    mix = np.dot(A, np.array([x, y]))     for i in xrange(len(mix[0])):        mix[0][i] = mix[0][i] + mix[1][i]     #show_data(mix[0], None)     return mix def decode(w,x):    ps = np.dot(x, w)    return ps[0], ps[1]  if __name__ =="__main__":    s1 =[math.sin(float(x)/20) for x in range(0, 1000, 1)]    s2 = [float(x)/50 for x in range(0, 50, 1)]* 20    #s1 = [math.sin(float(x)/20) for x inrange(0, 100, 1)]    #s2 = [float(x)/50 for x in range(0, 50,1)] * 2    #show_data(s1, s2)     # 假定有N个人在说话，则在任何一个时刻，原始的源信号是N维的，混合信号也是N维的。从而，混合矩阵A和解混矩阵W都是N*N维的。——这个维度，和观测时间M无关。    # 即：混合矩阵A必须是N*N维，才能使得n维的原始信号乘上A后可以得到个n维的混合信号，解混矩阵同理。    A = [[0.6, 0.4], [0.45, 0.55]]  # 混合矩阵    x = mix(A, s1, s2)  #s1/s2线性加权得到输入数据x    # 去均值，这个老师的效果不错，但我这里反而变差了....    # 实际运用中还是要根据实际场合看看是否选取    x_mean = np.mean(x)    for i in xrange(x.shape[0]):        for j in xrange(x.shape[1]):            x[i][j] = x[i][j] - x_mean    # ica    w = ica(x)    # 解混    [ps1, ps2] = decode(w, x)    show_data(ps1, ps2)