ICA

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

题目

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

如下图所示:

下面开始解题。

题目整理

首先将信息整理如下

源信号,如下图所示:

PS:源信号是未知的(知道还求啥),但为了分析还是需要把源信号表示出来。

第一个人是s1,他第一个时刻发送的信号是s11,第二个时刻发送的信号是s12,以此类推第m个时刻发送的时刻是s1m。第二个人,第三个人,….,第n个人同理。

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

观测信号,如下图所示:

解读方式和源信号一样。

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

然后,我们假设源信号经过线性加权得到了观测信号, 即:sm乘上一个权值A后得到了xm,即:xm = Asm。接下来为了表示所有的sm和xm,就用s代表sm组成的矩阵,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=A1,S=S1;乙求得A=-A1/2,S=-2S1,那甲乙的结果都能得出源信号,但甲和乙求得的结果的顺序就是反的(有个负号),振幅甲是乙的1/2。

进行假设

首先,先做两个假设:

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

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

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

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

开始推导ICA

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

假定第i个源信号的概率密度函数为pi(s),第j时刻的n个源信号记做向量sj,则在j时刻向量sj 的联合密度为:

根据xj = A·sjhttp://blog.csdn.net/xueyingxue001/article/details/51915390中“复合函数的概率密度”的知识,得:

从而得到似然函数:

进一步得到对数似然:

从而,目标函数为:

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

假定源信号的累积概率密度函数是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)对wij求偏导:

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

既然对wij求偏导是上面的结果,那对整个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 sys
import time
import math
import datetime
import threading
import numpy asnp
importmatplotlib.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)
 
时间: 2024-10-05 12:57:06

ICA的相关文章

ICA (独立成分分析)

介绍 独立成分分析(ICA,Independent Component Correlation Algorithm)简介 X=AS X为n维观测信号矢量,S为独立的m(m<=n)维未知源信号矢量,矩阵A被称为混合矩阵. ICA的目的就是寻找解混矩阵W(A的逆矩阵),然后对X进行线性变换,得到输出向量U. U=WX=WAS 过程 编辑 (1)对输入数据进行中心化和白化预处理 X*=X-u 经过白化变换后的样本数据为 Z=Wz X* (2)从白化样本中求解出解混矩阵W 通过优化目标函数的方法得到W

ICA Proxy数据访问流程

在Citrix的解决方案当中,都采用Citrix的NetScaler设备进行ICAProxy,将虚拟桌面和虚拟应用发布到外网.接下来详细说明ICA Proxy的访问数据流程. 如上图所示,其表示的步骤详细说明如下: 客户端通过HTTP GET请求发送访问的FQDN 地址到NetScaler的VPN虚拟服务器. NetScaler的VPN虚拟服务器经过检测确定传入的连接类型是未认证的客户端HTTP请求,然后将该客户端的HTTP请求重定向到认证页面(状态HTTP 302:其认证页面为:/vpn/in

UFLDL上的ICA为什么一定要做PCA whiten

UFLDL上的ICA为什么一定要做PCA whiten Andrew Ng先生的UFLDL教程真可谓deep learning入门的首选课程.在两年前我看教程里讲ICA部分的(链接)时候,里面提到使用教程所述的ICA模型时,输入数据必须经过PCA白化操作,页面上有个TODO问为什么要这样做.以当年的我对机器学习的理解并不能解答这个问题,就只是按照教程上讲的写完了代码,后来就一直没有看过了.  今天在与人讨论无监督学习的几种损失函数的时候,提到了PCA的损失函数: max∥Wx∥2s.t.WWT=

Citrix ICA协议简要介绍

关于Citrix的ICA协议,他的英文全称,网上都称呼其为Independent ComputingArchitecture,翻译为中文就是独立计算体系结构.但是根据Citrix的内部材料显示,也可以称为是ICA = Intelligent Console Architecture!翻译为中文就是智能控制台架构! 1.       历史 ICA 1.0 – 1992 在ICA1.0的版本中,最初是基于串行连接开放的,后来添加了IPX和NetBIOS的支持.所以在ICA1.0版本中,支持串行.IP

Windows 2008 想自动下载, default.ica 文件

Windows 2008  想自动下载, default.ica 文件,Windows 2008 ,IE8 默认不允许下载将Internet,下载启用即可,如图:

【转载】独立成分分析(Independent Component Analysis)ICA

独立成分分析(Independent Component Analysis) 1. 问题: 1.上节提到的PCA是一种数据降维的方法,但是只对符合高斯分布的样本点比较有效,那么对于其他分布的样本,有没有主元分解的方法呢? 2.经典的鸡尾酒宴会问题(cocktail party problem).假设在party中有n个人,他们可以同时说话,我们也在房间中一些角落里共放置了n个声音接收器(Microphone)用来记录声音.宴会过后,我们从n个麦克风中得到了一组数据,i表示采样的时间顺序,也就是说

ICA学习笔记

1. 独立成分分析(ICA)的起源: 独立成分分析(Independent Component Analysis),最早应用于盲源信号分离(Blind Source Separation,BBS).起源于“鸡尾酒会问题”,描述如下:在嘈杂的鸡尾酒会上,许多人在同时交谈,可能还有背景音乐,但人耳却能准确而清晰的听到对方的话语.这种可以从混合声音中选择自己感兴趣的声音而忽略其他声音的现象称为“鸡尾酒会效应”. 2. ICA的定义:(ICA的不确定性——在没有任何先验知识的前提下,是无法同时确定w和s

机器学习 —— 基础整理(四):特征提取之线性方法——主成分分析PCA、独立成分分析ICA、线性判别分析LDA

本文简单整理了以下内容: (一)维数灾难 (二)特征提取--线性方法 1. 主成分分析PCA 2. 独立成分分析ICA 3. 线性判别分析LDA (一)维数灾难(Curse of dimensionality) 维数灾难就是说当样本的维数增加时,若要保持与低维情形下相同的样本密度,所需要的样本数指数型增长.从下面的图可以直观体会一下.当维度很大样本数量少时,无法通过它们学习到有价值的知识:所以需要降维,一方面在损失的信息量可以接受的情况下获得数据的低维表示,增加样本的密度:另一方面也可以达到去噪

使用ICA进行EEG滤噪

数据预处理 假设6通道EEG(4通道EEG+2通道EOG),生成ndarray或矩阵S,row为time point,col为通道.对S进行标准化: S /= S.std(axis=0) ICA分析 模块导入: from sklearn.decomposition import FastICA 构建ICA对象,选择成分数: ica = FastICA(n_components=6) (成分排列随机,必要时使用random_state进行可重复性分析) S对应的成分为S_: S_ = ica.fi

网页客户端无法打开citrix ica后缀文件

网页浏览器无法打开citrix ica后缀文件,一直会转到下载页面,要求下载 解决方法 该情况是由于安装了迅雷等下载软件,默认点击是下载而不是直接打开. 可选择迅雷禁止浏览器监测,或者在迅雷下载界面有个直接用ie打开. 假设ica 文件没有设置默认打开程序,则也会出先这种情况.出现这种情况,请将文件保存至桌面,右键选择打开方式,手工选择默认打开程序默认程序为\Citrix\ICA Client\wfcrun32.exe