2007年12月27日星期四

Tips:2007.12工作点滴

1.MATLAB官方博客:http://blogs.mathworks.com/
其中的Inside the MATLAB Desktop 主要介绍工作环境和绘图中的新功能和技巧。推荐订阅。
最近的一篇"Focusing on Zooming" 一文中介绍了使用命令来控制figure的缩放。通过linkaxes命令,可以将一个figure中的各个axes指定相同的坐标范围,实现各个绘图同步放大、缩小。在缩放一个绘图时,其它也会响应地同步缩放。可以同时关联x、y轴,也可以只关联一个坐标轴。
Giving your code some privacy一节介绍了MATLAB2007b中引入的代码折叠功能,对编写大型的程序应该很有帮助。

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

3.判断点在离散的封闭曲线(多边形)内的方法,基本原理是,多边形内一点发出的任意射线与多边形仅有奇数个交点;
1). 求出封闭曲线的外接矩形的范围。即曲线中所有顶点的最大、最小坐标范围。
2). 判断点是否在该矩形内,否,则结束;是,则执行3)
3). 判断由该点发出的水平(或垂直)射线与多边形的交点个数。该射线与多边形的边有交点的判断(以水平正向为例):
边的两个端点的纵坐标不能同时大于或小于该点纵坐标,否则,该射线与该条边没有交点。
如果两个端点的横坐标同时小于该点的横坐标,则没有交点,如果同时大于该点的横坐标,则有交点。
如果一大一小,则需要求出该射线与边的交点的横坐标(纵坐标已知),如果交点横坐标大于点的横坐标,则有交点。(如果曲线上的点很密集,判断精度要求也不高,这一步也可以忽略)。

用上面的方法,在MATLAB中手动进行锋电位特征的聚类:


4. Word中图片使用嵌入型时,位置自动移动到页面的边缘:可能是样式引起,将图片的样式改为正文或可解决。实际上
可能是图片所在段落的行距太小造成的,可以设为单倍行距试试。

2007年12月21日星期五

sudo的简单配置

sudo的简单配置

sudo的配置文件是/etc/sudoers,它有专门的编辑工具visudo,用su切换到root用户,然后执行这个命令。
对于VI,对我这样的初学者,只要记住几个快捷键就可以了:光标移动,j上一行,k下一行,h向左,l向右。i切换到编辑状态,然后再配合End、Delete、和BackSpace键,完成这个位置的编辑,Esc退出编辑状态;移动光标到下一个位置,再切换到编辑状态,再输入,直到修改完成。然后输入冒号“:”,最下面一行会改为接受命令状态,输入w,保存,q退出。
执行visudo以后,将光标下移,可以看到一行:
root ALL=(ALL) ALL
这个命令定义了root用户的sudo权限,我们要添加的配置也具有相应的命令格式即
user host=[(Runas)] [NOPASSWD] commands
被授权用户 主机=[目标用户][不需密码] 命令
即“被授权用户”可以在“主机”上以“目标用户”的权限执行“命令”。方括号内是两个可选参数,目标用户缺省为root。命令必须用绝对路径,之间用","分开。例如
liuxq localhost=/sbin/poweroff
表示用户liuxq可以在本机上以root的权限执行"sudo /sbin/poweroff"命令,而不需要root密码(需要liuxq的用户密码)。
如果加上NOPASSWD,则表示不需要输入密码:
liuxq localhost=NOPASSWD: /sbin/poweroff

基本的定义方法就是这样,如果需要为不同的用户、不同的主机配置不同的命令权限,这样逐条配置显然太繁琐,所以这个配置文件中引入了别名Alias的功能,在配置文件的开头部分,有四个别名类型:
User_Alias; Host_Alias; Runas_Alias; Cmnd_Alias;
分别对应该命令形式中的四个部分。别名实际上定义了相同类型的对象的一个集合。例如
User_Alias LIU=liuxq,lxq
定义了一个LIU集合,其中有成员liuxq和lxq,用LIU代替前面命令中的liuxq,则用户liuxq和lxq都可以执行sudo /sbin/poweroff命令。
其它别名的定义也是类似的。这里就不多啰嗦了,推荐看看Linux 系统中的超级权限的控制,讲得比较详细。

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算法通过迭代方法,构造更好的 θ12,... ,逐步优化初始的 θ0, θ的递推优化公式通过下式得到:

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

Q(θ)是对数似然函数lnp( 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,每一个 yj 所从属的分布用 zj 表示,zj 的取值范围为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


2007年12月10日星期一

MAXScript:利用3DSMAX脚本产生原子运动动画

利用3DSMAX脚本产生原子运动动画

很久,很久以前,曾经做过一些关于分子动力学的工作,最后无疾而终,几乎什么都没有留下。分子动力学(Molecular Dynamic Simulation)计算微观原子间的相互作用,产生其运动轨迹,所有原子的运动形态最终表现为物体宏观的运动和变形等。

如果将这些原子的运动以动画的形式渲染出来,无疑是非常形象的。我曾经用3DSMAX的脚本MAXScript编写了一段程序来完成这个功能,在这里贴出来,供有需要的参考。
以下图中的模型为例,上面是Au探针,压入下面的Ni平板。原子的初始坐标值(x,y,z)分别放在文本文件posau0.txt和posni0.txt;
计算过程中产生一系列的原子坐标文件,这些文件以编号区别,如posau1.txt~posau100.txt;
在脚本文件中,首先读入原子的坐标信息,建立模型,每个原子对应一个球模型;
然后,在时间轴的不同时间坐标上,读入并指定每个原子的新坐标。
脚本的执行我也记不清了。可以在网上搜一下。
完整的程序如下(“--”的注释是把程序粘贴过来后补上去的,不知道有没有什么问题。我是在3DSMAX 5上运行的。):

--指定数组变量;
posarray = #()
postemp = #()
spAu = #()
spNi = #()
lay = #()
--字符串变量
strti = " "
--Au模型中原子个数,坐标数据量
atomau = 768
dataau = atomau*3
--Ni模型中原子个数,坐标数据量
atomni = 616
datani = atomni*3
--位置坐标文件的个数,不含初始坐标文件
framenum = 14
--原子半径
Rau = 1.44
Rni = 1.3
--指定材质,Ni平板每一层具有不同的颜色
--textMatAu = standardMaterial diffuse:[125, 125, 255]
textMatNi = standardMaterial diffuse:[92, 216, 92]
textMat1 = standardMaterial diffuse:[126, 173, 126]
textMat2 = standardMaterial diffuse:[233, 114, 67]
textMat3 = standardMaterial diffuse:[236, 224, 92]
textMat4 = standardMaterial diffuse:[92, 92, 216]
textMat5 = standardMaterial diffuse:[60, 151, 159]
textMat6 = standardMaterial diffuse:[179, 121, 216]
--读入Au的初始坐标
filename = "F:\data\posau0.txt"
in_name = (filename)
in_file = openFile in_name

if in_file != undefined then
(
for v=1 to dataau do
(
vert = readValue in_file
posarray[v]=vert
)
close in_file
--根据原子坐标建立每个原子的球模型
for v=1 to atomau do
(
ti=(v-1)*3+1
px = posarray[ti]
ti=ti+1
py = posarray[ti]
pz = posarray[ti+1]

spAu[v] = sphere radius:Rau position:[px, py, pz] segs:12 smooth:true
--判断原子所在的层,指定不同的材质
lay = (v-1)/288 as integer
lay = lay+1
if lay == 1 then
(spAu[v].mat = textMat1 )
if lay ==2 then
(spAu[v].mat = textMat2 )
if lay ==3 then
(spAu[v].mat = textMat3 )
if lay ==4 then
(spAu[v].mat = textMat4 )
if lay ==5 then
(spAu[v].mat = textMat5 )
if lay ==6 then
(spAu[v].mat = textMat6 )

)

)

--读入Ni的原子位置,并建立模型
filename = "F:data\pos\ni0.txt"
in_name = (filename)
in_file = openFile in_name
if in_file != undefined then
(
for v=1 to datani do
(
vert = readValue in_file
posarray[v]=vert
)
close in_file

for v=1 to atomni do
(
ti=(v-1)*3+1
px = posarray[ti]
ti=ti+1
py = posarray[ti]
pz = posarray[ti+1]

spNi[v] = sphere radius:Rni position:[px, py, pz] segs:12 smooth:true
spNi[v].mat = textmatNi
)

)

--开始设置坐标-时间动画
animate on (
at time 0
for tj in 1 to framenum do
( ti=tj*10
--指定新坐标对应的时间
at time ti (
strti = tj as string
filename = "F:\data\posau"+strti+".txt"
in_name = (filename)
in_file = openFile in_name
if in_file != undefined then
(
for v=1 to dataau do
(
vert = readValue in_file
posarray[v]=vert
)
close in_file

for v=1 to atomau do
(
ti=(v-1)*3+1
px = posarray[ti]
ti=ti+1
py = posarray[ti]
pz = posarray[ti+1]
--设置Au原子模型的新位置
spAu[v].position = [px, py, pz]
)
)

filename = "F:\data\posni"+strti+".txt"
in_name = (filename)
in_file = openFile in_name
if in_file != undefined then
(
for v=1 to datani do
(
vert = readValue in_file
posarray[v]=vert
)
close in_file

for v=1 to atomni do
(
ti=(v-1)*3+1
px = posarray[ti]
ti=ti+1
py = posarray[ti]
pz = posarray[ti+1]

spNi[v].position = [px, py, pz]

)

)
)

)
)

--max time play


再贴一个效果图(模拟计算结果本身可能不正确):