Mar 052015
 

设置代理服务器 Setting up proxy 出于对患者信息的保护,学校内网接入Internet时必须使用代理(Proxy)。 这令很多软件的安装使用过程变得复杂。 这个帖子会列出给常用软件设置代理的方法。 R 最直接的设置代理的方法是设置环境变量: http_proxy, https_proxy, 例如: 如果上述方法不起作用,也有可能需要: 这是在curl不能正常使用代理,但wget可以用的时候的解决方案。 还有一种情况是用devtools安装R的扩展包,可以用: 如果Mac电脑上没有安装X11,装package时需要指定repository(来源),比如: Docker 很多软件提供了Dockerfile。用户需要使用docker build来制作自己的容器(docker)。 如果在防火墙后面, 简单的方法是用(link): 另外一个方法是手动修改Dockerfile,加入ENV语句,例如:

Feb 192015
 

2015 春节 2015 Chinese New Year 爆竹声中一岁除,春风送暖入屠苏。 千门万户曈曈日,总把新桃换旧符。 ——王安石的《元日》 除夕夜里,稍稍反思一下2014年经历的投稿。 背景:我们写的软件速度快,功能全。比如在文章中,我们强调该软件可以处理很大规模的数据,这是已知其他软件达不到的。 审稿:除去离题万里的审稿意见以及错误的审稿意见,我们得到的反馈是对自己的软件过分称赞,对有竞争关系的其他软件过分贬低。同时这些意见往往对我们有强烈的负面情绪。 反思:比较软件从某种意义上就像是比较孩子。说自己的孩子好没问题,但要是说自己的孩子比别人的好,就易引起强烈的情绪反扑。即使是数字上的比较(比如身高说自己的孩子高于其他小孩XX厘米),有意见的人也可能会说:身高不是重要的标准,你比较的方式不对等等。这个经历告诉我指出其他人/事在科研上的不足应该非常谨慎。或许寻找“和气”、容易接受的比较方式是有必要的。 2015年我们会继续努力,写有价值的文章,做好科研,也希望能通过温和的方式,得到其他科研同行的肯定。此外,希望看到这里的各位 新春快乐!

Feb 032015
 

R函数式编程 Functional Programming in R R语言的循环很慢,同时语法有很复杂,陷阱很多。举例来说,逐行遍历一个矩阵有两种写法: (1) (2) 这两个写法那个对?或者两个都对? 熟悉R的用户自然会用(1),因为当m矩阵的行数为零时,只有(1)的结果是正确的: 为了避免手动写循环,一种折中的解决方案是用Filter, Reduce, Map等函数式编程的概念。 这些函数可以比较自然的处理列表(list)或者向量(vector)类型。举个例子: 可以看到,对于列表(list)或者向量(vector)类型,R内建的函数式函数(更多的见[1])还是很方便的。 如果想把函数式编程玩出更多花样,可以看看Hadley的新书Adv-R(参见[2])。 这里还有个问题,这种函数式的写法速度如何?下面的代码告诉我们函数式编程节省了写程序的时间,在某些情况下也许也能节省计算时间。 下面是结果: 在上面的例子中, 第一种是向量化的赋值,速度最快。 第二种是用循环,速度最慢, 第三种是函数式,速度居中。 不过据我的经验,当v的维数很小时(比如只有100个数),函数式相对于循坏的优势越来越小,甚至比循环要慢。 题外话:据说R的核心是类似于Scheme,也是一种函数式语言。但是在速度上对函数式编程的支持看来还有很大的提升空间。 [1]https://stat.ethz.ch/R-manual/R-devel/library/base/html/funprog.html [2]Advanced R http://adv-r.had.co.nz/Functional-programming.html

Dec 262014
 

拉格朗日乘数检验 Lagrange Multiplier test (Score test) 拉格朗日乘数检验,英文是Lagrange multiplier test,或者叫做Score test是一种常用的统计检验。 拉格朗日乘数检验的名称来源于这个检验用的是拉格朗日乘数的分布,见[2]。 Score test的名称则来自于Score本身。 为了写起来方便,下面都用Score Test来代替拉格朗日乘数检验。 Score Test,Likelihood Ratio Test和Wald Test的图形表示 假设似然函数 只有一个参数, 这三种检验可以在一张图里表示出来: (摘自:http://www.ats.ucla.edu/stat/mult_pkg/faq/general/nested_tests.htm) Likelihood Ratio Test计算的是, 就是两倍红色竖线的长度,这个统计量近似有自由度为1的卡方分布(假设只有一个自由变量) 这里的是最大似然估计(Maximum Likelihood Estimator): $$ \xi^R = 2 (l(\hat{a}) – l(0)) = 2 log(L(\hat{a})/L(0)) \sim \chi^2_1 $$ Score Test只需要考虑 这一点上似然函数的性质,在图像上通过蓝色部分表示(通过考察这一点的斜率和曲率)。 写成公式是: $$ \begin{align} 斜率 &= U = l'(0) […]

Dec 112014
 

加快C++编译速度 Speed up C++ compiling speed 在开发rvtests时,常常需要重新编译整个项目,但是这是个耗费时间的事。 比如,编译最耗费时间的regression库,用“make”,需要50秒。 在多CPU的机器上,打开并行编译选项后,“make -j”,需要11秒(系统时间是11秒,用户时间是64秒)。 这就是说重新编译regression库,最快也要11秒。 今天我发现有一个好用的工具,ccache, 可以把这个时间缩小到5秒 (加速编译速度一倍)。 方法是用: ccache首次使用会稍微慢一点。因为它会缓存(cache)源程序。但之后编译的速度就变得飞快。 比如上面的方法用5秒。如果还想更快,打开“-j”选项,只需要不到1秒。

Nov 172014
 

最近学习和研究有点松懈,转一篇古文激励一下自己。 小时候读这篇文章的时候没什么感受,现在读了就觉得应该立刻干活。 为学一首示子侄 彭端淑 天下事有难易乎?为之,则难者亦易矣;不为,则易者亦难矣。人之为学有难易乎?学之,则难者亦易矣;不学,则易者亦难矣。 吾资之昏,不逮人也,吾材之庸,不逮人也;旦旦而学之,久而不怠焉,迄乎成,而亦不知其昏与庸也。吾资之聪,倍人也,吾材之敏,倍人也;屏弃而不用,其与昏与庸无以异也。圣人之道,卒于鲁也传之。然则昏庸聪敏之用,岂有常哉? 蜀之鄙有二僧:其一贫,其一富。贫者语于富者曰:“吾欲之南海,何如?”富者曰:“子何恃而往?”曰:“吾一瓶一钵足矣。”富者曰:“吾数年来欲买舟而下,犹未能也。子何恃而往!”越明年,贫者自南海还,以告富者,富者有惭色。 西蜀之去南海,不知几千里也,僧富者不能至而贫者至焉。人之立志,顾不如蜀鄙之僧哉?是故聪与敏,可恃而不可恃也;自恃其聪与敏而不学者,自败者也。昏与庸,可限而不可限也;不自限其昏与庸,而力学不倦者,自力者也。 选自《白鹤堂文集》,作于乾隆九年(公元1744年)。

Sep 082014
 

拉普拉斯方法(Laplace’s Method) 拉普拉斯方法又称为拉普拉斯近似(Laplace Approximation)。它可以用来计算一元或多元积分[1]。 举例来说,假设 \(f(x)\) 是一维函数,我们要计算 . 如果\(f(x)\)形式很复杂,我们往往找不到定积分的公式(Close form)。 如果用数值方法来计算积分,计算量又很大。 所以要想个办法,得到一个比较精确的积分结果。 拉普拉斯方法可以适用于这种情况。我们先用泰勒展开(Taylor Expansion): $$ f(x) \approx f(x_0) + f'(x_0)(x-x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$ 如果选取\(x_0\)使得\(f'(x_0) = 0\),则可以进一步简化为 $$ f(x) \approx f(x_0) + \frac{1}{2}f”(x_0) (x-x_0)^2 $$ 再引入一个假设,\(\int_{-\infty}^{\infty} f(x) \mathrm{d} x \) 可以被写成 \(\int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x \) 的形式,那么就有: $$ \begin{align} \int_{-\infty}^{\infty} e^{f(x)} \mathrm{d} x &\approx \int_{-\infty}^{\infty} […]

Jul 062014
 

卡方检验 Chi-square test 卡方检验的全称是Pearson’s Chi-square test,这是统计中最重要的检验。 从用途来讲, 卡方检验可以做两件事: (1)适配度检定(Goodness of fit),检查样本是否符合某种随机分布; (2)独立性检定(Independence test),检查多个变量之间是否独立。 这篇不讲怎么用卡方检验,而是讲讲卡方检验的公式,推导,以及和似然检验的联系,目的是温故知新。 从公式来讲卡方检验: $$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} $$ 这里是观察到的次数,是期望的次数, 统计量 服从自由度的卡方分布( distribution)。 假设有一个有可能是多项分布的随机变量。 观察到类的次数分别是. 为了统计检验是否服从这个分布,列出零假设 为服从多项分布,对立假设为不服从这个多项分布,可以得到: $$ \chi^2 = \sum_i \frac{(O_i – E_i)^2}{E_i} = \sum_i \frac{(O_i – N \cdot p_i)^2}{N \cdot p_i}$$ 注意这个每一项的分母是, 但根据中心极限定理(Central Limit Theorem), $$\frac{O_i – N \cdot […]

Apr 132014
 

R package的UBSAN测试 How to perform UBSAN tests on R packages 最近需要提交一个R package (seqminer)到CRAN,我们需要保证这个包代码通过UBSAN 测试(Undefined bahavior sanitizer tests,详见这个链接). 这个测试可以有效的检测结果不确定的指令(例如打印INT_MAX+1)。 目前只有Clang(版本3.3以后)才能支持UBSAN。为了对R package进行UBSAN测试,我们首先需要一个特殊的R版本,即用打开UBSAN支持(-fsanitize=undefined)的Clang编译的R。参考了网上一篇Blog,写一下简单步骤: 1. 下载编译工具 2. 下载R源代码(devel版本) 3. 替换R-devel/config.site 4. 编译 这里–without-recommended-packages会省略一些R的packages(例如boot,nlme,Matrix) 以提升编译速度。如果不用这个选项,我发现Matrix package总是没法编译 5. 建立~/.R/Makevars 这里CXX 一定要用clang++,如果用clang的话会出现链接错误,比如找不到这个符号: _ZTVN10__cxxabiv117__class_type_infoE 这里还需要指定PKG_LIBS包括LLVM的库文件,不然就找不到UBSAN的符号。 6. 测试package UBSAN测试的结果会显示在:seqminer.Rcheck/tests/tests.Rout 如果有错误,你会看很长的一行, runtime error…. 其他技巧 检测内存还可以用valgrind:

Mar 092014
 

A pitfall in R 先看下面的代码: 最后一行输出的是115, 不是116。为什么? 答案是: altAllele是一个numeric类型的数,尽管打印出来是116,但在电脑里存储的是一个比整数116稍小的浮点数。因此rep(1, altAllele)是一个长度是115的vector。 R的程序有很多陷阱。这个陷阱是因为R不需要定义数据类型,但它的底层还是有类型的。因此写R程序一定要注意这点。解决方法如下: Python里有没有类似的陷阱?结论是一样有 (Python 2.7.5 NumPy 1.8.0): 但是Python里打印altAllele的值不是116,而是115.99999999999999。这总比R打印出16要好一点。