SAS 快速概览

以前在几个分析里用过 SAS,最近又需要用这货,然而 SAS 的语言设计实在是太陈旧了,和习惯差得有点远,官方文档又比较冗长,所以弄一个快速概览用于 rapid prototyping。

Basics

Language

SAS 的脚本是由很多个分块组成的,每个分块官方叫做步骤 (step)。 步骤以其名字开头,以运行指令 run; 结尾。

每个步骤中有若干行指令,叫做语句 (statement),以分号结尾。

下面是一个只有两个步骤的示例程序

/* Comments are enclosed like this. */

/* This is a example of a step called "data" */
data work.exampledataset;
    /* There are 3 statements in this step */
    i = 1;
/* Generally all steps require a "run;" statement 
to indicate the end of the step */
run;

/* Another step called "proc", with 2 statements */
proc freq data=work.exampledataset;
run;

步骤有很多种,不过关键的步骤只有两个,data 和 proc。 data 步用于创造数据,proc 步可以调用一些预先写好的程序用于处理数据,相当于 call。

不难发现,在 SAS 脚本的最上层,也就是步骤的这一层,从设计上是不存在真正的编程的,它实际上更类似于 shell script,只是用于表示数据处理的流程。 SAS 支持一些 macro,因此一定要说的话也确实可以在这一层编程,不过显然并不是非常友好。

主要的真正编程是在 data 步和 proc iml 中。 data 步里支持一些常见的 control flow 和 data structure,可以创建变量、数组等。 proc iml 是一个特殊的 proc 步,相当于 SAS 的线性代数库,可以完成一些简单的矩阵运算等。 这两个步中的 statement 比较接近我们一般说的编程语句。

在其它一般的 proc 步中,比如 proc freq,里面的 statement 主要就是设定参数。 官方出于某种原因故意混淆了几种非常不同的意思,容易造成一些混乱。 值得注意的是,许多 proc 步中的参数,有的是紧跟在名字后面,有的则是要放在步骤中的 statement 中。 大体上看,比较关键的,并且短一些的参数似乎更倾向于出现在名字后面,而复杂一些的参数则倾向于出现在 statement 中,不过也有很多反例,具体哪个参数要放在何处,需要参阅官方的文档。

data step

我们首先介绍 data 步。

data 步以 data example_name; 语句开头,意思是创建一个名为 example_name 的数据集。

奇妙的数据结构

data 步中可以创建变量,只要用 i=1; 这样的语句就可以。 同样也可以创造数组,用 array x{100} _temporary_; 可以创建一个 1*100 的数组。

要注意的是,SAS 中所有变量与传统意义的变量不同,它实际上是数据库表的一个列名。 换句话说,所有的变量都是数组,其中的一个元素叫做一个观测(observation)。 我们用 i=1; 创造了一个变量,实际上得到的是一个数组,其第一个元素是 1。 因此,我们创造一维数组的时候,不难想到我们创造的实际上是二维数组。

接下来当我们更改 i 的值的时候,会有一些奇异的事情。 前面说到,变量实际上是数组。 在我们初次创建这个变量之后,i 实际上指向的是数组的第 1 个元素。 因此,如果我们此时修改 i 的值,我们会修改到数组的第 1 个元素。

“output;” 语句会保存当前的所有变量,在当前所有”指针”的位置后方插入一个新元素,并且将指针指向这个新元素。

因此,如果我们运行了一次 output; 语句,则我们在数组 i 的第 1 个元素后面插入了一个新元素,并且将指针指向这个新元素即第2个元素。再修改 i 的值时,修改的是数组的第 2 个元素。

最后,当 run; 指令运行的时候,数据集会被创造,其内容就是所有的变量对应的数组。

我们来看一个例子

data work.exampledataset;
	i=1;
	output;
	i=2;
	i=3;
	output;
	j=i;
	output;
run;

它的输出是

Obsij
11.
23.
333

以上的运算规则还不算最为诡异的。 SAS 中的 data 步提供了一些读入外部数据集的办法。 “set datasetname;” 语句可以读入一个之前创造的数据集。 刚才我们说到,当新建一个变量 i 时,i 指向数组的第 1 个元素。 那么,当我们读入的数据集有一个列名是 i 时,i 指向哪个元素呢? 答案是,i 指向数组的所有元素。

一个有用的指令是 “put varname;” ,它可以将某个string或者变量的值输出到 log。 如果我们在上一个例子中 put i; 则会在 log 中输出一个数。 而如果我们将上一个例子得到的数据集读入,然后 put i; ,输出是

1
3
3

可见此时 i 选中了这个数组的所有元素。 那么,考虑一个有趣的问题。这个时候运行 output; 会发生什么? 答案是,所有被选中的元素后面都会被插入一个新元素,然后被选中的元素变成这些新插入的元素。

来看这个例子

data work.exampledataset;
	i=1;
	output;
	i=3;
	output;
run;

data work.data2;
	set work.exampledataset;
	output;
	j = i + 1;
	output;
run;

它的输出是

Obsij
11.
212
33.
434

我们从上面的例子还可以看到,任何算数符号都是 element-wise 的,比如说等号就会将当前选中的所有值依次拷贝到另一个变量中,并且保持位置不变。

最后,有没有办法把指针向前移动呢? 也就是说我们怎样返回去修改某个变量数组之前的值? 答案是没有办法,因为从逻辑上看,每次 output; 之后,所有当前变量的值都已经被保存为了一个 observation,既然已经输出出去了,那就不应该再重新拿回来修改了。

奇妙的数组

我们再来看 SAS 的数组,我们刚刚说定义数组的语句是

array x{100} _temporary_ ;

那么其中的 _temporary_ 是什么意思呢? 实际上,按官方的文档,定义数组的方法是

array x{3} name1 name2 name3;

也就是说,其实数组是要给每个元素起名字的。那显然这个根本不是什么数组,一般我们好像把这种结构叫做字典。

更骚的是,我们定义了数组之后,如何访问其中的元素呢? 实际上我们根本不用管数组名 x ,直接用元素的名字比如说 name2 就可以访问到数组 x 的第二个元素。 所以我们这个数组实际上可能都不能算字典,数组指令的用处实际上是用来批量创建很多个变量。 这种奇妙的数组实际上是因为 SAS 保存数据集的时候不允许某个列没有列名造成的。

可以看这个例子

data work.asdf;
	array x(3) x1 x2 x3 (1,2,3) ;
	x2 = 4;
    output;
run;

它的输出是

Obsx1x2x3
1143

那我们怎么才能定义一个比较像是真的数组的东西呢? 官方提供了一个选项,也就是我们上面用的 _temporary_ ,加上这个选项之后,就不用为每个元素指定名字了,要访问其中的元素时也只能用下标,并且,在输出时这个数组也不会被输出出去。

作为例子,筛法求质数

data whatever;
	ARRAY pn{10000} _TEMPORARY_;

	DO i=2 TO dim(pn);
		pn(i)=i;
	END;

	DO i=2 TO dim(pn);

		IF pn(i)=. THEN
			CONTINUE;

		DO j=i*2 TO dim(pn) BY i;
			pn(j)=.;
		END;
		output;
	END;

	DROP j;
run;

打包windows下软件的一些坑

前一阵子答应帮了个忙,开发一套仪器软件。写了一天之后写好了,交付源码,但是对方几乎没有懂计算机的人,所以要求我把所有东西再打包成一个 windows 上的可执行文件。这几天刚好回来,稍微有一点空闲时间,随手打了一下包。不过,因为很多年没有用windows,实在是很不熟悉,打包windows软件的时候遇到了很多奇怪的小坑,稍微总结一些,方便自己以后参考。

Continue reading “打包windows下软件的一些坑”

解决 VSCode format tex 文件时报错 Can’t locate YAML/Tiny.pm in @INC

这个错误的原因很简单,就是 YAML::Tiny 没安装,VSCode 安装 LaTeX Workshop 插件的时候不会自动安装,因此手动安装一下就可以解决。

cpan Unicode::GCString
cpan App::cpanminus
cpan YAML::Tiny
perl -MCPAN -e 'install "File::HomeDir"'

注意如果 perl 是 brew 安装的,每次 perl 升级后,会需要重新来一次。

By default non-brewed cpan modules are installed to the Cellar. If you wish
for your modules to persist across updates we recommend using `local::lib`.

为了避免麻烦,可以

PERL_MM_OPT="INSTALL_BASE=$HOME/perl5" cpan local::lib
echo 'eval "$(perl -I$HOME/perl5/lib/perl5 -Mlocal::lib=$HOME/perl5)"' >> ~/.bash_profile

如何DIY一台信噪比和分辨率不输于商业仪器的低成本激光Raman光谱仪

Raman光谱是什么我就不多解释了,大家都懂,相信大家也都已经体验过用Raman光谱仪做各种分析的方便之处了,比如说自己新搞来或者合成了什么物质想简单地看一下纯度,再比如说用来快速地鉴定各种塑料材料是什么,或者鉴定一些样品里面有什么东西等等。总之,拉曼光谱最爽的几个点就在于它不需要对样品做过多的处理,晶体、固体粉末、液体、溶液、表面等等都可以直接打谱,可以说是极其的方便。

然而普通的人想要搞一台Raman放在家里用并不是一件容易的事,因为目前来看商业的Raman光谱仪价格都还比较高,便宜的在几十万,贵一点的一套房就没有了。然而,实际上商业的拉曼里面有很多功能我们都是用不上的,我们日常使用的话对拉曼的要求并没有那么高,最重要的参数就是分辨率和信噪比(SNR)。实际来看,分辨率做到10个波数左右,SNR上百就已经非常够用了(作为参考,Wasatch Photonics 的 WP532 系列的分辨率有3种,分别是7、9、15个波数,常温探测器的信噪比在500:1,低温探测器的信噪比在2400:1)。

这篇文章主要是从原理出发,带你走完一个设计到搭建一台基本的拉曼光谱仪的完全流程,包括光学部分,电路部分和软件部分等等。在这台基本的拉曼光谱仪上,可以做各种改动来做更高级的拉曼,比如探测器冷却技术,表面增强技术等等。

Continue reading “如何DIY一台信噪比和分辨率不输于商业仪器的低成本激光Raman光谱仪”

简单介绍一下生物体系中的相分离和凝胶化,它们之间的相互影响,以及Intrinsically Disordered Linkers和蛋白结构域在其中的作用

生物体系相分离是最近生物物理研究中的新兴方向,目前理论研究还较少。比较成体系的是 Michael Rosen 和 Rohit Pappu 等人的工作,我校的来鲁华老师组最近也对这方面有一些关注。这篇文章主要介绍Pappu等人在去年所做的计算模拟工作,稍微结合一些他之前的文章。

Harmon et al., eLife, 2017, 6:e30294

S. Banjade et al., Proc. Natl. Acad. Sci., 2015, 112 (47), E6426-E6435

Rahul K. Das,  Rohit V. Pappu, Proc. Natl. Acad. Sci.2013,  110(33), 13392–13397

文中所引图片、资料均为 CC License 或不超过合理限度的评述性引用。

Continue reading “简单介绍一下生物体系中的相分离和凝胶化,它们之间的相互影响,以及Intrinsically Disordered Linkers和蛋白结构域在其中的作用”

约化密度矩阵理论简介


\(
\def \bra#1{{\langle #1 \rvert}}
\def \ket#1{{\lvert #1 \rangle}}
\def \braket#1#2{{\langle #1 | #2 \rangle}}
\def \expect#1{{\langle #1 \rangle}}
\)如何用约化密度矩阵方法计算量子体系能量等问题,最早在约六七十年前就已经提出。但由于N-representability问题的悬而未决,在过去的几十年中一直没有解决。自上世纪90年代来,随着控制理论中semidefinite programming优化方法的发展,虽然N-representability问题依然没有完全解决,但利用p-positivity条件已经可以由2-RDM计算得到较好的结果。另一种不同但与之相关的方法,即采用Contracted Schrödinger equation的SCF解来进行计算,在3-RDM和4-RDM近似展开在2-RDM上的条件下,也能得到不错的结果。

虽然这种方法在近年发展缓慢,甚至有可能是一种不太可行的方法(非常感谢刘文剑老师课上提醒,N-representable 条件可能是无穷多的!),但这种方法的formulation其实非常漂亮,其背后思想十分简洁。作为一种 alternative viewpoint,也许能对我们有一些启发。

这篇文章主要是结合一下这方面工作的代表性人物 Mazziotti 教授 (JFI, UChicago) 写的一些综述,以及他在07年出的一本书《Reduced Density Matrix Mechanics: with Application to Many-Electron Atoms and Molecules》(见Advances in Chemical Physics第134卷),简单介绍一些这方面的基础。

Continue reading “约化密度矩阵理论简介”

一个入门水平的简单统计力学教程[未完工]


这篇文章大概是从高中的时候就打算写,不过考虑到长度过长,来北大之后就一直没有时间动笔。 最近刚好考完了一波期末,两学期的物理化学课程也结课,准备趁着摸鱼的这小段时间简单地介绍一下统计力学当中比较浅薄的概念。 (顺便:强力推荐化院高毅勤和刘剑两位老师合授的物理化学英文版课程, 010-01035200-1006172506-1 ,化院两年以来上的最有价值的课,虽然我是个菜逼,依然学到很多)

对于我个人来说,中学时期接触统计力学思想之前,学习经典热力学漫长又缺乏清晰内核的经历是相当痛苦的。 因此我认为从微观的角度出发,尤其是从分子动力和量子力学的角度出发,结合微观与宏观的关系来理解热力学量的成因是一条更为清楚、更加富有内涵的道路。 并且,理解微观与宏观的联系方法 (eg. Fokker-Planck equation),对于社会统计、微观/宏观经济学、生态学等其它学科的理解也是至关重要的 (as an interesting example, see this and this

作为一个初级的统计力学教程兼一篇安利性质的文章,我认为不应该有过多的具体实现细节,并且希望尽量写得轻松一点。 作最简陋的打算,文章主要会把握住微观与宏观的联系这条主线,考虑相应简单地介绍数学背景和物理图像。 在此之外,如果我还有时间完善文章,作最理想的打算,则物理化学内容方面,希望尽量包括化院本科的化数、物化、中级物化课的初级内容和其它学校一些研究生课程的内容; 算法和代码方面,打算尽量讲清思路,并结合我自己的编码经验,偶尔会给出python, julia, MATLAB 或者类 C 代码的伪实现; 数学方面,希望能摆脱初级微积分的技巧性操作,尽量从抽代、实变、泛函、几何之类的角度,建立高观点来看待复变、概率、ODE、PDE等等需要用到的工具中的技巧,试图寻找比较优雅的理解问题的方法。

Continue reading “一个入门水平的简单统计力学教程[未完工]”