liuxqsmile‘s Blog: MATLAB

程序界面中包含三个部分，菜单、信息表格和底部的状态栏。软件的所有功能均通过菜单命令来操作。信息表格列出了分析数据的主要信息，如采样频率，记录长度，文件名称、锋电位数目、特征名称、各个分类包含的锋电位数目等。对于多通道数据，通过信息表格右上角的Channel下拉菜单选择当前通道。状态栏给出程序运行状态的简明信息。

2 数据文件

SPKtool没有和其它软件进行数据交换的接口，不能直接读入记录系统产生的数据。用户需要首先将数据（连续记录信号、锋电位波形矩阵等）导入MATLAB，并按照特定的变量名保存到.mat文件（MATLAB中默认的数据文件格式），才能在SPKtool中打开，进行分析。

在菜单中选择File->Load…，打开文件选择对话框。加载数据后，根据文件中的变量，自动设置各个菜单的状态。

SPKtool通过变量名区分数据，因此，对数据文件(.mat)中的变量名有严格的要求。

2.1原始记录数据

原始记录数据文件中至少应包含两个变量：

dataflow：连续记录数据，二维数组，行对应通道，列对应采样点；

fs：采样频率，数值。

默认情况下，SPKtool对读入的数据进行低通滤波，使用6阶Butterworth滤波器，截至频率4500Hz。可以打开MLoadRaw.m，找到butter函数对应的位置，修改滤波参数。

2.2锋电位波形数据

从原始信号中检测出的锋电位波形数据，至少应该包含以下三个变量：

posind：锋电位波形在连续数据中的采样位置标记。行向量。第n个波形发生的时刻为posind(n)/fs。

waveforms：锋电位波形数据，二维数组。每行对应一个锋电位的采样数据。例如，1000个波形，每个波形含40个数据点，则waveforms尺度为1000×40；

fs：采样频率。

数据文件中还可以包含如下一些变量：

featurespace：锋电位特征空间，二维数组，每行包含一个波形的特征。如1000个波形的PC1和PC2特征，则featurespace为1000×2。特征的数目不限。

featnames：特征名称，包含字符串的单元数组，每个字符串对应一个特征的名称。例如: 两个特征的名称分别为PC1和PC2，则 featnames 为{‘PC1’,’PC2’};

unitcell：包含分类信息的单元数组。每个单元对应一个类别中的锋电位序号列表。例如，unitcell{1}对应没有归类的锋电位编号。 unitcell{3} 是第二个类别中的锋电位编号数组。

对于多通道数据，上面这些变量需要放在一个名为ResultSet的结构体数组内，每个结构体对应一个通道。

刺激信息

stimmark: 刺激的开始和结束时间。行向量。例如[20 30 60 70]，表示两次刺激，第一次刺激从第20s开始，到30s结束，第二次刺激60s开始，70s结束。

上述变量可以共存于一个文件中。在程序中，如果找不到dataflow这个变量，将关闭锋电位检测功能（不能激活Detect菜单），找不到waveforms，将关闭特征提取功能（不能激活Features菜单）。

另外，程序中用到了两个全局变量，采样频率fs和定义各个类别色彩的矩阵colorary，所以在程序运行过程中，请勿执行clear all等能够清除全局变量的命令。

2.3数据保存和导出

选择菜单File->Save …，保存分析结果。文件格式同样是.mat文件，变量名称与读入数据相同。（注意，连续记录数据不包含在该数据文件内。）可以用SPKtool重新打开保存的文件，也可以直接用读入MATLAB，将特定的数据导出。

选择菜单File->Save Current Channel…，仅保存当前通道的分析数据。

菜单File-> Experimental Info …，设置试验信息，包含记录位置、时间、坐标、数据文件名称等信息。

菜单File-> Export to TXT file …，导出数据到文本文件。这将产生两个文本格式的文件，后缀名为.inf的文件包含上面的试验信息，以及数据信息（通道数、单元数等）。后缀名为.txt的文件是一个数据表格，第一列对应锋电位时间（即变量posind），第二列为类别编号，后面的各列为波形数据。

3 锋电位检测

如果读入的数据中包含原始连续记录信号（即dataflow变量），菜单Detect可用，其中包含以下两个命令：

3.1峰值密度分布曲线

菜单Detect->Plot Peak Histogram，绘出原始记录数据的峰值密度分布曲线。幅度为对数坐标。

3.2锋电位检测界面

菜单Detect->Load Detector…，加载锋电位检测，界面如下：

参数设置：

Detect On：该下拉菜单中包含了三种波形检测方式，即分别在原始信号（Raw Signal），原始信号的能量(Energy)和非线性能量(Nonlinear Energy)^¹曲线上进行检测；

Threshold：阈值。输入阈值后左图中的红色水平线也相应调整。信号中与阈值线有交点的波形将根据后面的参数进行截取；

Spike length：每个锋电位波形的长度，单位为ms。如果采样频率为40KHz，1ms长度的锋电位波形包含40个数据点；

Prethreshold ratio：对齐点以前数据长度比例。对齐点前保留的数据长度与波形数据总长度的比值。例如包含40个点的锋电位波形，对齐方式为valley时，在波谷位置前保留10个采样点，则该参数需要设置为10/40=0.25；

Dead time：不应期长度，单位为ms。例如，如果采样频率为40KHz，该参数设置为2时，程序在检测到一个波形后，将跳过该波形后紧邻的80个点，然后才开始寻找下一个波形。

Alignment：波形对齐方式，阈值位置（gate）、波谷 (valley)、和波峰 (peak) 对齐三种形式。阈值和波谷的对比如下：

a) gate b) valley

使用方法：

中间的坐标轴根据Detect On的设置，显示原始信号或其（非线性）能量曲线。底部坐标轴在（非线性）能量检测方式下显示原始信号波形。左侧坐标轴显示中间坐标轴中数据的峰值分布曲线（横坐标为对数坐标）。

可以利用顶部的工具栏，对中间坐标轴进行缩放和平移，左侧坐标轴的纵坐标和底部坐标轴的横坐标也将进行相应的调整。

参数设置完毕后，选择Detect进行波形检测，检测完成后，将打开一个单独的绘图窗口显示锋电位波形。如果不符合要求，可以修改参数，重新检测。最后，选择Done完成，并关闭该界面。

需要说明的是，由于原始记录信号中的数据量通常很大，为了加快运行速度，坐标轴中只绘出前0.4e6个采样点。当然，锋电位波形检测还是针对所有的采样数据。

右侧的参数中除Detect On外，具有自动保存功能，再次运行Detector时将自动加载前一次设置的参数，即使你前一次选择Cancel退出。

在切换不同的检测方式时，阈值将根据中坐标轴中的数据自动调整为：均值+/-3倍标准差。

4 锋电位特征提取

通过阈值检测出锋电位波形后，或者读入的数据中保护锋电位数据，这时，Feature菜单将可用。在该菜单中，可以进行锋电位特征的选择、计算、修改、绘图等功能。通常需要定义特征后，才能打开Sort菜单进行聚类操作。

4.1特征定义

菜单 Features->Feature Define…，打开如下的特征选择对话框：

左侧All Features 列表中列出了所有可用的特征。中间有两个特征列表：Features Vs Features 中的特征用于绘制特征比较图（菜单 Features-> Features Vs Features）；Features For Cluster 才是用于聚类分析的特征。本说明中，除非特别指明，特征均指用于聚类的特征。在右侧可以修改一些特征相关参数。

可用特征：

peak：锋电位波峰处的幅值；

peak.pos：峰值对应的位置；

valley：锋电位波谷处的幅值；

valley.pos：波谷对应的位置；

peak-valley：峰峰值，波峰到波谷的振幅。

width：波形宽度，即波峰位置到波谷位置的距离：peak.pos-valley.pos；

energy：波形的能量，即采样点的平方和。

nonlinear.energy：所有点的非线性能量之和，幅度平方与频率平方的乘积，其计算方法为：NEO(t) = v(t)²- v(t-1)*v(t+1)，（v(t)为记录信号）；

Slice：锋电位波形中某个时刻的幅值，对应参数Slice Parameters。例如每个锋电位含32个采样点，Slice Parameters为0.5，则Slice特征为锋电位中第32×0.5=16个采样点的数值。Slice Parameters可以包含多个值，提取多个采样点，例如 0.5，0.75 则第16和第24个采样点的数值都作为特征。

pca：主成份分析方法，可以在右侧PCA Parameters中设置要提取的主成份序号；

ridge.pca：基于小波脊的PCA特征提取方法，首先对波形逐个进行小波变换（使用mexh小波），然后提取出小波脊，再对小波脊上的小波系数进行PCA，提取特征。需要三个参数：小波分析的频率范围Frequency limits，小波脊的数目Number of Ridges，以及需要提取的主成份序号PCA Parameters。

完成特征选取后，选择Done，开始提取用于聚类分析的特征，并给出特征分布散点图。

4.2特征修改

特征比例缩放

Features->Scale Features，打开如下图的对话框，在其中，可以调整特征的相对范围，可能有利于改善K-Means的聚类结果。

删除特征

Features->Remove Features，删除部分使用中的特征。一般特征的计算不需要花很长时间，重新选择并计算即可。在同时使用Ridge.pca和PCA作为特征时，由于共用PCA参数，可以通过该工具重新选择。

4.3特征分布散点图

Features->Plot Features 绘制特征分布散点图，如果定义了两个特征，散点图是二维的，如果有三个或更多的特征，则只使用前三个，以三维形式显示，并同时画出特征向三个平面上的投影分布，如下图。

4.4特征密度图

Features-> Features Density：特征分布密度图。如下图，该图只能绘出前两个特征的分布密度。

4.5特征对比图

Features-> Features Vs Features，绘制特征分布的对比图，注意该功能中使用的特征与聚类分析的特征是独立的。并且，每次绘图都需要重新计算特征。

如果已经划分类别，不同类别的特征将以不同颜色区别显示。

4.6主成份

Features-> Principal Components：给出特征提取中主成份曲线。例如特征定义中使用第一个和第二个主成份，则给出PC1和PC2的形状曲线。

5 聚类分析

聚类功能集中于菜单Sort中。在选择具体的聚类方法之前，应该先设置聚类的参数。

5.1参数设置

菜单Sort->Options…, 打开参数设置对话框：

其中共有四个参数：

Unit Number：划分类别的数目，K-Means，期望最大化(EM)方法，以及模板匹配方法使用该参数；

K-Means Radius：K-Means方法的聚类半径，值越小，聚类的半径越小。模板匹配方法也受该参数影响；

Parzen Mult：寻谷聚类算法的初始半径大小，初始半径越小，获得的类别数目越多。

Sort On：设置聚类的锋电位范围，选择All Spikes时对所有锋电位的特征进行聚类；选择Unsorted Spikes时仅对未分类的锋电位进行聚类。

5.2手动聚类

菜单Sort->Contour…，对应手动聚类，选择该菜单后，弹出如下图的手动聚类对话框。

对话框中间给出了未分类的特征的分布图。底部有四个控制按钮。首先单击New Contour按钮，光标将变成图中所示的十字形，在特征分布图中单击左键，将出现一个红色的十字形标记点，这样依次点击，就可以确定一个类别的大致轮廓，之后，鼠标右击，这些点将自动连接成一条封闭的曲线。然后选择底部的Update按钮，曲线内的点就被划分到一个类别。再选择New Contour绘制下一个类别的边界，依此操作，直到聚类完成。

也可以连续绘制多条边界线，然后再选择Update，显示聚类结果。最后选择Done将显示的聚类结果传回主程序，选择Cancel放弃。

手工聚类的类别数目不受前面聚类参数对话框中的类别数目的影响。

5.3自动聚类方法

SPKtool中包含了几种常见的聚类算法，简单说明如下：

K-Means

根据特征之间的距离来进行聚类。用户指定聚类的数目（聚类参数对话框中），该算法首先从样本（特征）中随机抽取确定聚类中心，然后计算每个样本到各个聚类中心的距离，并将其归到距离最近的类中，第一次聚类完成后，计算各个样本的均值作为新的聚类中心，重新计算每个样本到新的聚类中心的距离，并将其划分到聚类最近的类别中，如此循环，直到聚类中心的位置稳定^[^²^]。

Template 模板匹配法

在SPKtool中，模板方法的与K-Means基本上是一样的，不同之处仅在于直接用锋电位的波形数据作为特征，而K-Means使用锋电位降维后的特征数据。

Gaussian EM 基于高斯混合模型的期望最大化算法

该算法首先用K-Means方法确定初始聚类，然后将锋电位特征的分布看成符合混合高斯分布，用期望最大化算法确定高斯混合模型的均值、协方差等参数，然后计算每个特征在各个高斯分布中的概率，并划分到概率最大的一个类别（高斯分布）中^[2]。

Valley Seeking 寻谷聚类算法

在特征空间中寻找特征分布密度的极小位置，作为聚类的边界。算法的流程可以简述如下^³:

1. 确定初始半径，可以取为样本标准差的某个倍数；初始半径的设置影响聚类的数目；

2. 在该半径范围内，确定每个特征的近邻数目，以及近邻的集合。如果特征的数目为N，需要初始化一个N×N的数组来存放每个特征的近邻编号，计算完成后，可以根据最大近邻集合中元素的数目，删去多余的元素。

3. 如果一个特征的近邻数目，超过它每一个近邻各自的近邻数目，则该特征确定了一个特征密度峰。如此，设置一个最小近邻标准，可以找出全部峰。

4. 五个峰确定5个类别，给每个峰的近邻特征设置类别标记；

5. 将所有特征按照近邻数目从大到小排序；

6. 从近邻数目最大的特征开始，看其近邻中属于哪一类别的最多，就将该特征也归入这一类；如此，完成全部特征的聚类。

5.4聚类结果调整

移除未分类的锋电位

菜单Sort->Edit units …中，选择复选框Remove unsorted spikes。未分类的锋电位将从数据中删去。

合并分类

同样在菜单Sort->Edit units …中，在Merge Units框中输入欲合并的单元编号；例如，输入“2 3”将单元2和单元3合并为单元2，输入“0 3”删除单元3（其中的波形转为未归类状态）。

移除特征离散度很大的锋电位

可以首先通过轮廓线聚类方法，将要保留的特征划分到单元1中，然后按照上面的方法，先移除没有归类的边缘特征，再将单元1与0（未分类）合并。

6 绘图功能

6.1锋电位波形

菜单：Plot->Spikes：绘出所有锋电位波形；

菜单：Plot-> Sorted spikes：分布绘出各个类别的波形。每个波形图左下的数字标明了该分类中波形的数目。

6.2特征分布图

菜单：Plot-> Features：绘出特征分布散点图；

菜单：Plot->Sorted Features：画出分类后的特征散点图；如下图。

6.3时间间隔(ISI)分布

菜单：Plot->ISI …，打开如下的ISI参数设置窗口：

Min. Interval: 最小间隔时间；

Max. Interval: 最大间隔时间；

Number of Bins: 阶梯图的分段数目，数值越大，ISI曲线越平滑；

Units: 选择单元，只画出指定单元的ISI分布图。

6.4 间隔-时间分布图

菜单：Plot->ISI versus time…，绘制相邻放电间隔的时间分布图。参数设置对话框：

X Min. :横坐标（时间）最小值；

X Max.：横坐标最大值；

Units：选择需要进行绘图的单元编号。

6.5庞加莱图

菜单：Plot-> Poincare maps：该图以与前一个放电的时间间隔为纵坐标，与后一个放电的时间间隔为横坐标，标出每一个放电的位置。坐标轴的范围（即最大时间间隔）由ISI中的Min. Interval和Max Interval 指定。如下图。

6.6放电频次阶梯图

菜单：Plot->Fire Rate…：绘制放电的频次直方图。参数设置对话框：

X Min. :横坐标（时间）最小值；

X Max.：横坐标最大值；

Time per Bin：每个阶梯表示的时间长度；

Units：选择需要进行绘图的单元编号。

Add stimulation markers：选中该项后，可以在下面的Markers中输入刺激标记（对应数据文件中的simmark变量）。例如，刺激开始的时刻为49.89885s，结束的时刻为262.4599s，在Markers中填写如下：

放电频次分布直方图如下，刺激时间将以一个黑框标示出来。

6.7放电光栅图

菜单：Plot->Rasters …：绘制放电的光栅分布图。参数设置对话框：

X Min. :横坐标（时间）最小值；

X Max.：横坐标最大值；

Units：选择需要进行绘图的单元编号。

Add stimulation markers：与频次直方图中相同。

6.8相关函数图

菜单：Plot->CrossCorrelograms…：绘制单元间的互相关和自相关函数。参数设置对话框：

Pre-ref. Time：横坐标（时间）最小值，计算相关函数时信号的偏移范围；

Post-ref. Time：横坐标最大值；

Number of Bins：相关函数分段数目；

Units：参与计算的单元编号；

Ref Units：作为参考的单元编号。

6.9刺激前后放电直方图PSTH

菜单：Plot->Peristimulus Histogram … ：打开如下的参数设置对话框：

T. pre-Event：刺激事件前保留的时间长度；

T. post-Event：刺激事件后保留的时间长度；

Time/Bin：直方图的每段包含的时间长度；

Units：选择需要进行绘图的单元编号。

Stimulus Phase：选择刺激事件的相位，“on”表示刺激开始；“off”表示刺激结束；

Stimulus No.：当有多个刺激周期时，可以选择其中的一个刺激周期。设置为0时，将所有刺激周期的的PSTH叠加。

6.10刺激前后放电光栅图

菜单：Plot->Peristimulus Raster … ：打开如下的参数设置对话框：

T. pre-Event：刺激事件前时间长度；

T. post-Event：刺激事件后时间长度；

Units：选择需要进行绘图的单元编号。

Stimulus Phase：选择刺激事件的相位，“on”表示刺激开始；“off”表示刺激结束；

Stimulus No.：当有多个刺激周期时，可以选择其中的一个刺激周期。设置为0时，将所有刺激前后放电都画出。

在该组信号中，共有三个刺激事件，Stimulus No.设置为0时，三个刺激前后的放电都绘出，如下图。纵坐标名称中，E*表示事件的编号，U*表示单元编号：

部分参考文献：

^{1Mukhopadhyay S, Ray GC. A new interpretation of nonlinear energy operator and its efficacy in spike detection. IEEE Transactions on Bio-Medical Engineering. 1998 Feb ;45(2):180-7.}

2 （美）迪达等著，李宏东等译. 模式分类（原书第2版）[M]. 北京: 机械工业出版社，2003.

3 Zhang C, Zhang X, Zhang MQ, Li Y. Neighbor number, valley seeking and clustering. Pattern Recogn. Lett. 2007 ;28(2):173-180.

其它文献：

梁培基，陈爱华. 神经元活动的多电极同步记录及神经信息处理 [M]. 北京：北京工业出版社, 2003.

Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potential [J]. Net Comput Neural Syst, 1998, 9(4): 53-78.

Answer.com: Expectation-maximization algorithm

2008年6月18日星期三

如何将曲线以不连续的形式画出来

最近有些“忙”， blog也没有什么内容好更新。弄点简单的东西，凑个篇幅。
MATLAB中的数据，用plot画出来一般是一条连续的曲线，怎么用将它画成离散的形式呢？
例如下面的连续正弦曲线，怎么画成后面的离散形式呢？

这个方法很简单，这个上图中的连续曲线：
x = 1:100;
y = sin(x/20);
figure, plot(x,y);
不连续曲线：
y1 = y(1:2:100);
y2 = y(2:2:100);
x1 = x(1:2:100);
x2 = x(2:2:100);
figure, plot([x1;x2], [y1;y2]);
简单点来说，就是将每段曲线的起点和终点的横纵坐标分别放在一个2×N的数组中。
在这个不连续曲线中，每个线段都是一个绘图对象，所以具有不同的颜色。这种方法的限制是每条连续线段必须包含相同的数据点数。否则的话，只能通过循环来逐段绘制。
实际上，如果数据中有NaN（即：Not a Number），该点是不会绘出来的，这样也可以形成不连续效果。
仍然是上面的曲线：
yn = y;
yn(3:3:100) = NaN;
figure, plot(x, yn)
得到下图。这种方法就没有了各段尺度相同的限制。简单点来说，就是在需要断开的数据点（包括横纵坐标）前面插入NaN。这样的不连续线仍然是作为一个整体的。

2008年4月13日星期日

MATLAB中的正态分布检验

要对一组样本进行正态性检验，在MATLAB中，一种方法是用normplot画出样本，如果都分布在一条直线上，则表明样本来自正态分布，否则是非正态分布。
MATLAB中也提供了几种更正式的检验方法：
kstest Kolmogorov-Smirnov 正态性检验，将样本与标准正态分布（均值为0，方差为1）进行对比，不符合正态分布返回1，否则返回0；该函数也可以用于其它分布类型的检验；
lillietest Lilliefors test。与kstest不同，检验目标不是标准正态，而是具有与样本相同均值和方差的正态分布。
jbtest Jarque-Bera test。与 Lilliefors test 类似，但不适用于小样本的情况。

2008年3月31日星期一

MATLAB中的一些小技巧(2)

MATLAB中的一些小技巧(2)

1. Ctrl+D打开子程序
在MATLAB的Editor中，将输入光标放到一个子程序名称中间，然后按Ctrl+D可以打开该子函数的m文件。当然这个子程序要在路径列表中（或在当前工作路径中）。实际上该快捷键对应右键菜单中“Open Selection”命令。

2. 列出依赖函数
在Editor菜单中，选择Tools->Show Dependency Report，可以列出当前m文件调用的子程序。如果需要将程序移动到别的电脑上，使用这个工具可以确保不会漏掉一些东西。这个工具也适用于文件夹。

3. 显示鼠标数据内容
调试程序过程中，程序运行到断点后，将鼠标移到一个变量名上，将会弹出一个方框，显示该变量的内容。不在调试状态下，怎么使用这个功能呢？在File菜单中打开设置选项（Preferences），选择 Editor/Debugger -> Display，在 "General Display Options" 中选择 "Enable datatips in Edit Mode "。

4. 单元数组的基本操作
清除cell中的元素

cellname(n1:n2) = [];

合并如下的两个cell数组

a = {'str1', 'str2', 'str3'};
b = [1:3];

c = [a(1:2), b, a(3)];

得到：c = { 'str1' 'str2' [1x3 double] 'str3'}

而 c = {a(1:2), b, a(3)}

则得到 c = {1x2 cell} [1x3 double] {1x1 cell}

cell array本质上是数组，前一种操作方式与普通数组的操作方式相同，而后一种方式则是将逗号分隔的各个元素分别一个cell，整体作为一个cell数组。

cell的两种引用方式
a(1) 得到 {'str1'}, 因为a是一个cell数组，所以其中的每一个元素都是一个cell，
而a{1} 得到 'str1'。

5. MATLAB 中的字体平滑
在File菜单中打开设置选项（Preferences），在 Fonts 中，选择底部的 “Using antialiasing to smooth desktop fonts” ，MATLAB 中默认字号比较小，可以调大一些（可以在下拉菜单中输入数字，更精细地控制文字大小）。有衬线字体打开反锯齿后可能会看起来“发虚”。

6. 同步缩放坐标轴linkaxes
linkaxes命令，可以将一个figure中的各个axes指定相同的坐标范围，实现各个绘图同步放大、缩小。在缩放一个绘图时，其它也会同步缩放。可以同时关联x、y轴，也可以只关联一个坐标轴。

7. 鼠标绘图
MATLAB Help中的Interactive Plotting一文提供了一种在figure中用鼠标绘图的方法，主要使用ginput函数。

8. GUI程序关闭按钮的callback
为GUI窗口右上角的关闭按钮添加回调函数：在GUIDE中打开GUI，在没有控件的地方单击，选中整个figure，右键菜单中选择 View Callbacks -> CloseRequestFcn，m文件中将自动添加该回调函数，输入代码就可以了。

2008年3月25日星期二

MATLAB 中GUI子程序的参数传递

输入参数传递：
比如子GUI的名称为subGUI, 设想的参数输入输出为：[out1, out2] = subGUI(in1, in2)
在subGUI的m文件中（由GUIDE自动产生）：
1.第一行的形式为：function varargout = subGUI(varargin)
该行不用做任何修改；varargin 和 varargout 分别是一个可变长度的cell数组(MATLAB帮助文件中有说明)。输入参数in1和in2保存在varargin中，输出参数out1，out2包含在varargout中；
2.在subGUI的OpeningFcn中，读入参数，并用guidata保存，即：
handles.in1 = varargin{1};
handles.in2 = varargin{2};
guidata(hObject, handles);

返回参数的设置：
1. 在GUI子程序的OpeningFcn函数的结尾加上uiwait(handles.figure1); figure1是subGUI的Tag；
2. subGUI中控制程序结束（如"OK”和"Cancel"按钮）的callback末尾加上uiresume(handles.figure1)，不要将delete命令放在这些callback中；
3. 在子GUI的OutputFcn中设置要传递出去的参数，如 varargout{1} = handles.out1；varargout{2} = handles.out2;末尾添加 delete(handles.figure1); 结束程序。

在GUI的OpenFcn中，如果不加uiwait，程序会直接运行到下面，执行OutputFcn。也就是说程序一运行，返回值就确定了，再在其它部分对handles.output作更改也没有效果了。

加上uiwait后，只有执行了uiresume后，才会继续执行到OutputFcn，在此之前用户有充分的时间设置返回值。

通过以上设置以后，就可以通过 [out1, out2] = subGUI(in1, in2) 的形式调用该子程序。

在一个GUI中调用另一个GUI时，主GUI不需要特别的设置，同调用普通的函数一样。在打开子GUI界面的同时，主程序还可以响应其它的控件。不需要担心子GUI的返回值被传错了地方。

2008年2月19日星期二

MATLAB中的一些小技巧

MATLAB中的一些小技巧

1. Ctrl+C 中断正在执行的操作
如果程序不小心进入死循环，或者计算时间太长，可以在命令窗口中使用Ctrl+c来中断。MATLAB这时可能正疲于应付，响应会有些滞后。

2. figure命令新建一个绘图窗口
figure 可以打开一个空的绘图窗口，接下的绘图命令可以将图画在它里面，而不会覆盖以前的绘图窗口。当有多个figure窗口时，在命令窗口中执行如Plot等命令将覆盖当前figure窗口中的对象。所谓的当前figure窗口，也就是最后一次查看的窗口(可以用命令gcf得到)。
figure(N)，将编号为N的窗口置为当前figure，出现在所有窗口的最前面，如果该窗口不存在，则新建一个编号为N的空白figure。

3. 注释掉一段程序：%{、%}。
经典方法是用 if 0，但缺点是不够直观，注释掉的内容仍然保持代码的颜色。现在可以用 %和{ 的组合。输入%{后，后面的代码都会变绿，在注释结束的地方再加上 %} 。
也可以选中要注释内容，在右键菜单中选择Comment (Uncomment去掉注释标记)，或使用快捷键Ctrl+R。
将光标放在需要注释的行中，按Ctrl+R，将该行变为注释。取消注释也是一样的，快捷键为Ctrl+T。

4. doc 命令名，打开命令的帮助文档
与help命令不同，帮助文档MATLAB Help中对命令的描述更详细，往往还有一些例子，更便于理解。

5. clc 清屏
清除命令窗口中的所有输入和输出信息，不影响命令的历史记录。

6. clear 和clear all
clear 变量名：可以清除workspace中的无用的变量，尤其是一些特别大的矩阵，不用时及时清理，可以减少内存占用。
clear all 清除所有的变量，使workspace一无所有，当重新开始一次算法验证时，最好执行一次，让workspace中的变量一目了然。

7. close all
如果开了多个绘图窗口，用这个命令将它们一起关掉。

8. 上下光标键↑↓
在命令窗口中，上下光标键可以将历史记录中的命令复制到输入位置，便于快速重新执行。如果输入命令的前几个字母如 [row, col] =，再使用光标键，则只会选择以这些字母开始的命令。

9. Tab补全
对名称记得不太清楚的命令，可以输入开头的几个字母，然后按Tab键，当只有一个以这些字母开头的命令时，将自动补全命令名，否则显示一个命令名列表，方便从中选择。当然，只在命令窗口中有效。

10. cell模式
在一个长长的脚本m文件中，可能需要对其中的一段反复修改，查看执行效果，这时，cell模式就非常有用了。cell模式相当于将其中的代码拷贝到命令窗口中运行。两个%后接一个空格(%% )开始一个cell。将输入光标放到一个cell中时，背景将变为浅黄色，Ctrl+Enter执行cell中的代码。
执行cell中代码时不需要保存m文件，该m文件可以不在路径列表中。
cell模式中，断点不起作用，当然，调用的子程序中的断点还是正常的。

11. 获取文件列表，批处理
MATLAB内置了一些文件操作命令，例如cd(切换工作目录)，dir(同ls，显示目录内文件列表)等。dir命令可以返回目录中的文件和文件夹列表，存在一个结构体数组中。如果需要对一些数据文件进行批处理，而文件名又没有一定的规律，你可能需要借助于这个命令。

2007年12月13日星期四

Expectation-maximization：高斯混合模型的期望最大化算法

高斯混合模型的期望最大化算法

本文大部分内容翻译自Answers.com的Expectation-maximization algorithm 词条。

对于一个不完整的样本集，假设y为观测（已知）数据，z为丢失数据，z和y构成了完整的数据。z可以为丢失的观测量，也可以是隐含的变量，如果这些变量已知的话将会简化问题的求解。例如，在混合模型中，如果产生样本点的各个成分的分布为已知，则可以大大简化似然函数的计算。

假设p为完整数据的联合概率密度函数，该函数具有参数θ： $p( mathbf y, mathbf z | theta)$ ，它也可以看成是完整数据的似然函数（即样本集(y,z)具有θ分布函数的相似性），一个关于θ的函数。进一步，根据贝叶斯公式和全概率公式，丢失数据对于观测数据的条件分布为：

$p(mathbf z |mathbf y, theta) = frac{p(mathbf y, mathbf z | theta)}{p(mathbf y | theta)} = frac{p(mathbf y|mathbf z, theta) p(mathbf z |theta) }{int p(mathbf y|mathbf hat{z}, theta) p(mathbf hat{z} |theta) dmathbf hat{z}}$

这个式子仅要求计算观测数据对于丢失数据的条件概率（似然函数） $p(mathbf y|mathbf z, theta)$ 和丢失数据的分布概率 $p(mathbf z |theta)$ 。（分母仅为归一化需要）

EM算法通过迭代方法，构造更好的 θ₁,θ₂,... ，逐步优化初始的 θ₀， θ的递推优化公式通过下式得到：

$theta_{n+1} = argmax_{theta}Q(theta)$

Q(θ)是对数似然函数ln $p( mathbf y, mathbf z | theta)$ 的期望。换句话说，我们不知道所有数据的值，所以我们不知道似然函数的准确值，但是对于给定的观测数据y，我们可以得到未知的z的概率分布，对于每一个丢失的观测（z中的一个元素），都有一个关于θ的似然值，这样我们就可以计算似然函数的数学期望。这个期望显然受到θ的值的影响。

Q(θ)通过下式得到：

$Q(theta) = sum_z p left(z ,|, y, theta_n right) log p left(y, z ,|, theta right)$

或者，更普遍地，可以写成条件期望的形式

$Q(theta) = E_{mathbf z} ! ! left[ log p left(mathbf y, mathbf z ,|, theta right) Big| mathbf y right],$

应该清楚，这里关于 $log p left( mathbf y, mathbf z ,|, theta right)$ 的条件期望中用到的z的条件分布（前一个式子中的第一项）是在固定θ_n的情况下取得的。用对数似然代替直接的似然函数，可以在保证极值位置不变的情况下简化计算。也就是说，θ_{n + 1}是通过在θ_n 的基础上，最大化所有数据关于已知数据的条件期望E而取得的。

当样本集合的最大似然估计比较容易获得时，EM算法更为有用。如果似然估计具有给定的形式，则M步就非常简单。一个经典的例子是用于有限的高斯混合模型的最大似然估计，如果混合分布已知，则可以很容易得到各个成分的分布。

高斯混合模型

假设n个高斯分布，从中产生m个样本 $mathbf y_1, dots, textbf{y}_m$ ,每一个 y_j 所从属的分布用 z_j 表示，z_j 的取值范围为1到n。对于任意的y，其来自第i个高斯分布的概率为

$P(mathbf y | z=i,theta) = mathcal{N}(mu_i,sigma_i) = (2pi)^{-l/2} {left| sigma_i right|}^{-1/2} expleft(-frac{1}{2}(mathbf y - mathbf mu_i)^T sigma_i^{-1} (mathbf y - mathbf mu_i)right)$

我们的任务是估计未知的参数 $theta = left{ mu_1, dots, mu_n, sigma_1, dots, sigma_n, P(z=1), dots, P(z=n) right}$ ，即每一个高斯分布的均值、方差和先验概率。

首先我们必须设定初始值。

对于无监督的聚类，我们首先必须为每一个样本指定一个类别，这一步可以通过其他的聚类方法实现，如K-means方法，求出各个类别的中心和每一个样本的类别。然后求出各个类别中样本的协方差阵，可以用每个类别中样本的个数来表示该类别的权重（先验概率）。

E-step

使用上一次M步中得到的参数θ，计算未知的z（即每一个分类）对于观测数据的条件分布。即，在当前参数θ_n 下，对每一个样本，计算各个类别的后验概率。下式中第三个等号后面的分子的第一项，为每一个样本在类别i中的分布概率，分子中的第二项为每一个类别的先验概率。
分母为归一化，即对每一个样本，其在各个类别中的分布概率的和为1.

$p(z_j=i|mathbf y_j,theta_t) = frac{p(z_j=i, mathbf y_j | theta_t)}{p(mathbf y_j|theta_t)} = frac{p(mathbf y_j|z_j=i,theta_t) p(z_j=i|theta_t)}{sum_{k=1}^n p(mathbf y_j | z_j=k, theta_t) p(z_j=k|theta_t)}$

for ti = 1:Ncluster
PyatCi = mvnpdf(featurespace,Mu(ti,:),Sigma(:,:,ti));
E(:,ti) = PyatCi*Weight(ti);
end
Esum = sum(E,2);
for ti = 1:Ncluster
E(:,ti) = E(:,ti)./Esum;
end

M-step

为了最大化联合分布的对数似然函数的期望

上式中最后一项的联合分布可以写成条件分布和边缘分布的乘积的形式，得到

$Q(theta) = sum_{j=1}^m sum_{i=1}^n p(z_j=i | mathbf y_j, theta_t) ln left( p(mathbf y_j | z_j=i, theta) p(z_j=i | theta) right)$

上式需要满足先验概率归一化的要求：

$sum_{i=1}^{n} p(z_j=i|theta) = 1$

为了求得Q(θ)的极大值的位置，引入拉格朗日算子，然后带入概率密度函数，

$mathcal{L}(theta) = left( sum_{j=1}^m sum_{i=1}^n p(z_j=i | mathbf y_j, theta_t) left( - frac{l}{2} ln (2pi) - frac{1}{2} ln left| sigma_i right| - frac{1}{2}(mathbf y_j - mathbf mu_i)^T sigma_i^{-1} (mathbf y_j - mathbf mu_i) + ln p(z_j=i | theta) right) right) - lambda left( sum_{i=1}^{n} p(z_j=i | theta) - 1 right)$

为了求得新的 θ_{t + 1},我们需要找出上式取得极值时对应的θ，即其一阶导数 $frac{partial mathcal{L}(theta)}{partial theta} = 0$ 。

下面推导满足上式的均值，

for ti = 1:Ncluster
muti = E(:,ti)'*featurespace;
muti = muti/sum(E(:,ti));
Mu(ti,:) = muti;
end

产生新的协方差矩阵

新的各个分类的权重

进行归一化处理

Inserting λ into our estimate:

将新的参数，均值、方差、分类的先验概率，带入E步，重新计算，直到样本集合对各个分类的似然函数不再有明显的变化为止。

下面的程序提供了一种快速的似然函数计算方法，摘自Patrick Tsui的EM_GM样例程序：
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8636

function L = Likelihood(X,k,W,M,V)
% Compute L based on K. V. Mardia, "Multivariate Analysis", Academic Press, 1979, PP. 96-97
% to enchance computational speed
% a subfunction from Patrick P. C. Tsui's EM_GM code.
[n,d] = size(X);
U = mean(X)';
S = cov(X);
L = 0;
for i=1:k,
iV = inv(V(:,:,i));
L = L + W(i)*(-0.5*n*log(det(2*pi*V(:,:,i))) ...
-0.5*(n-1)*(trace(iV*S)+(U-M(:,i))'*iV*(U-M(:,i))));
end

订阅：博文 (Atom)