标签归档:统计计算

Breiman访谈实录

COS编辑部按:本文是一篇Richard Olshen对Leo Breiman的采访稿(原文发表在Statistical Science)。翻译工作已经得到作者授权。翻译: 张晔、成慧敏、李宇轩。审校:高涛、侯澄钧、丁鹏、魏太云。此外,郑重感谢施涛、丁鹏、郁彬老师为文章的翻译指导和版权沟通提供的帮助。

译者简介:张晔,毕业于华南统计科学研究中心,现严肃科技平台开发工程师,主要负责docker容器调度系统开发。成慧敏,就读于中央财经大学统计与数学学院,硕士研究生二年级,研究兴趣为复杂网络分析与深度学习。李宇轩,就读于中国人民大学统计学院,大二本科生,目测统计有关都可以是学习方向。

1928年1月28日,Leo Breiman生于纽约。5年后,他们家搬到了旧金山,然后Leo开始了他的学业。在他读初中的时候,他们家又搬去了洛杉矶。1945年,Leo从Roosevelt高中毕业后考进了加州理工学院,在那里他花了4年时间主修物理。1950年,Leo拿到了哥伦比亚大学的数学硕士学位,1954年,他又拿到了加州大学伯克利分校的数学博士学位。

Leo对科学和数学有着广泛的兴趣,包括信息论和博弈论。他曾参与汽车交通、空气质量和有毒物质识别等方向的研究。 他写过一篇著名的关于概率论的毕业论文,他是分类回归树(CART, Classification and Regression Trees)及其配套软件$CART^R$的四位作者之一,另外他还写了两本专著。Leo和Jerome Friedman一起开创了ACE(alternating conditional expectations)算法,该算法描述了因变量和自变量之间的非线性回归关系。 他开创性地提出将”bagging”和”arcing”这两种需要大量计算的方法用于分类,目前很多学者对此十分感兴趣。

Leo的职业履历包括,加州大学洛杉矶分校(UCLA)数学系教职,13年的独立咨询顾问,加州大学伯克利分校(UC Berkeley)统计系教授,同时也是该校统计计算实验室的创始人兼主任。 另外,他还是斯坦福大学和耶鲁大学的客座教授。 由于他的诸多贡献,Leo被授予数理统计研究所(Institute of Mathematical Statistics)和美国统计协会(American Statistical Association)的荣誉基金。 同时,他还是美国艺术与科学学院(American Academy of Arts and Sciences)的选举成员,并被加州大学授予Berkeley Citation荣誉奖项。

Leo Breiman是一个兴趣广泛的人,他不仅是专业的统计学家和概率学家,还在其他方面也取得了很多成就。他在Catskills当过服务员,在Merchant Marine当过洗碗工,同时他是一名探寻过热带雨林核心地带的背包客,是一群来自墨西哥农村孩子的慈爱父亲,是Santa Monica学校董事会的主席,是他美丽小屋的建筑师,还是一个技艺高超的雕刻家。Leo和他的妻子Mary Lou,居住在加州伯克利,他们育有两个女儿,Rebecca和Jessica。

采访者简介:Richard Olshen,斯坦福大学生物统计教授,生物统计方向(Division of Biostatistic)首席科学家,卫生研究和政策系(Department of Health Research and Policy)的副主任,斯坦福大学电气工程系和统计系兼职教授。该访谈于1999年2月19日在Leo和他的妻子Mary Lou的家中进行。 继续阅读Breiman访谈实录

Sweave后传:统计报告中的大规模计算与缓存

学无止境。我曾以为我明白了如何在Sweave中使用缓存加快计算和图形,但后来发现我并没有真的理解,直到读了另外一些手册才明白,因此本文作为前文“Sweave:打造一个可重复的统计研究流程”之续集,向大家介绍一下如何在Sweave的计算和图形中使用缓存,以节省不必要的重复计算和作图,让那些涉及到密集型计算的用户不再对Sweave感到难堪。

如果你还没读前文,建议先从那里开始读,了解Sweave与“可重复的统计研究”的意义。简言之,Sweave是一种从代码(R代码和LaTeX)一步生成报告的工具,我们可以把整个统计分析流程融入这个工具,让我们的报告具有可重复性。然而,就普通的Sweave而言,这样做的一个明显问题就是,所有计算和作图都被融入一个文档之后,每次运行这个文档都要重复所有的计算和作图,这在很多情况下纯粹是浪费时间;比如,我只想对新添加的部分内容运行计算,而文档中的旧内容希望保持不变。这都是很合理的需求,我们需要的实际上就是一种缓存机制,将不想重复计算的对象缓存起来,需要它的时候再从缓存库中直接调出来用。

Sweave本身做不到缓存(或者很难做到),但由于R语言的扩展性以及一些疯子和天才的存在,我们可以通过R的一些附加包来完成缓存。这里的缓存有两种(不关心细节的读者可以略过):

  • 一种是对R对象的缓存,比如数据、模型(拟合好的模型)和其它计算结果等,这个缓存容易理解,我们可以把它简单想象为把计算过程中的对象存入某些特殊的缓存库(文件),下次计算的时候先查看缓存库中是否存在我们需要的对象,如果存在,那么就直接加载进来,否则就重新计算;如(伪代码):
    ## 如果存在,就加载
    if (file.exists('x.RData')) {
        load('x.RData')
    }
    ## 如果当前工作环境中不存在x,那么就重新创建它
    if (!exists('x')) {
        x = rnorm(100000)
    }

    R对象的缓存可以通过cacheSweave包(作者Roger D. Peng)实现,它真正采用的方法比上面的伪代码当然要高级(比如用代码的MD5值作数据库名等),不过用户可以不必关心这其中的细节。这种“当需要时才加载”的方式在术语上叫“延迟加载”,即Lazy Load,这在R里面也很常见,比如很多R包的数据都是采用延迟加载的方式(加载包的时候并不立刻加载数据,而是用data(name)在需要的时候加载)。

  • 另一种是图形的缓存,这是一个不可思议的魔法。cacheSweave的帮助文档和vignette中特别提到了它只能缓存R对象,无法缓存图形和其它附属输出,而开源世界的魅力就是你总能找到一些奇妙的解决方案。R包pgfSweave借力于pgf实现了图形的缓存。pgfSweave包的两位作者Cameron Bracken和Charlie Sharpsteen都是水文学相关专业出身,却在R的世界里做了这些奇妙的工作,让我觉得有点意外,这是题外话。pgf是LaTeX世界的又一项重大发明,它提供了一套全新的绘图方式(语言),而且它的图形可以通过LaTeX直接生成PDF输出(这是pgfSweave能实现缓存的关键)。关于pgf,我们可以翻一翻它726页的手册,如果你能坚持看上10分钟,感叹超过50次,那么你一定是一位超级排版爱好者。它的图形质量之高、输出之漂亮,真的是让人(至少让我)叹为观止。同样,这里我也不详细介绍细节。pgf图形的“缓存”是通过“外置化”(externalization)来实现的,简言之,pgf图形有一种输出形式如下:
    \beginpgfgraphicnamed{graph-output-1}
    \input{graph-output-1.tikz}
    \endpgfgraphicnamed
    

    如果LaTeX在运行过程中遇到这样的代码段,那么它会把tikz图形(tikz是高层实现,pgf是基础设施)先编译为PDF图形,等下次重复运行LaTeX的时候,这段代码就不必重新运行了,LaTeX会跳过它直接引入PDF图形。这样就省去了编译tikz图形的时间,同时也能得到高质量的图形输出。

    pgfSweave从R的层面上可以输出这样的代码段,并生成一个shell脚本来编译图形为PDF;而pgfSweave也会尝试根据R代码的MD5来判断一段画图的代码是否需要重复运行——如果R代码没有变化,那么图形就不必重制。这是自然而然的“缓存”。对了,pgfSweave对Sweave的图形输出作了扩展,新增了tikz输出,这是通过另一个R包tikzDevice实现的(同样两个作者),这个包可以将R图形录制为tikz绘图代码。pgfSweave依赖于cacheSweave包,把它缓存R对象的机制也引入进来,所以我们可以不必管cacheSweave,把精力放在pgfSweave上就可以了。

介绍完血腥的细节之后,我们就可以坐享高科技的便利了。下面我以一个统计计算中的经典算法“Gibbs抽样”为例来展示pgfSweave在密集型计算中的缓存功能。

Gibbs抽样的最终目的是从一个多维分布$f(X_1,X_2,\ldots,X_n)$中生成随机数,而已知的是所有的条件分布$f(X_1|X_2,\ldots,X_n)$、$f(X_2|X_1,\ldots,X_n)$、……、$f(X_n|X_1,\ldots,X_{n-1})$,并且假设我们可以很方便从这些条件分布中抽样。剩下的事情就很简单了:任意给定初始值$x_1^{(0)},x_2^{(0)},\ldots,x_n^{(0)}$,括号中的数字表示迭代次数,接下来依次迭代,随机生成$x_1^{(k)} \sim f(X_1|x_2^{(k-1)},\ldots,X_n^{(k-1)})$、$x_2^{(k)} \sim f(X_2|x_1^{(k)},\ldots,x_n^{(k-1)})$、……、$x_n^{(k)} \sim f(X_n|x_1^{(k)},\ldots,x_{n-1}^{(k)})$(迭代次数$k=1,\cdots,N$),也就是,根据上一次的迭代结果并已知条件分布的情况下,从条件分布中生成下一次的随机数(这个算法很简单,你什么时候看明白上标下标了就明白了)。这是个串行的过程,所以无法避免循环。

接下来的计算没有那么复杂,我们仅仅用一个二维正态分布为例。当我们知道一个二维正态分布的分布函数时,我们也知道它的两个边际分布都是一维正态分布,而且分布的表达式也很明确,后文我会详细提到,读者也可以参考维基百科页面。生成一维正态分布的随机数对我们来说当然是小菜一碟——直接用rnorm()就可以了。这样,我们只需要交替从$f(X|Y)$和$f(Y|X)$生成随机数,最终我们就能得到服从二维联合分布的随机数(对)。假设我们要从下面的二维正态分布生成随机数:

$
\left[\begin{array}{c}X\\ Y\end{array}\right]\sim\mathcal{N}\left(\left[\begin{array}{c}0\\ 1\end{array}\right],\left[\begin{array}{cc} 4 & 4.2\\ 4.2 & 9\end{array}\right]\right)
$

详细过程和结果参见下面这份PDF文档(点击下载):

A Simple Demo on Caching R Objects and Graphics with pgfSweave (PDF)

二维正态分布随机数及其等高线图
二维正态分布随机数及其等高线图

我们生成了50万行随机数,并画了X与Y的散点图。由于我们设定了相关系数为0.7,所以图中自然而然显现出正相关;而等高线也体现出多维正态分布的“椭球形”特征。均值在(0, 1)附近,都和理论分布吻合。所以这个Gibbs抽样还不太糟糕。

生成上面的PDF文档和图形的LyX/Sweave源文档在这里:

A Simple Demo on Caching R Objects and Graphics with pgfSweave (LyX)

如果你已经按照我前面的文章配置好你的工具(即使当时配置过,现在也需要重新配置,因为我最近作了重大修改),这个文档应该可以让你重新生成我的结果。文档中有两处关键选项:

  • cache=TRUE 使得文档中的R对象可以被缓存;
  • external=TRUE 使得图形可以被缓存;

第一次运行文档时,也许需要等待一分钟左右的时间,因为R需要长时间的Gibbs抽样计算,并编译生成的图形为PDF,但如果下一次再运行这个文档,则速度会快很多,因为不需要再计算任何东西,只需从缓存加载就可以了。我们前一次的抽样结果已经在缓存库中了。

本文至此可以结束,但真正想享用这个结果的读者可能还要继续看几个问题:

  1. pgfSweave对图形的缓存经常有点过度,即使你改动了代码,缓存仍然不会被替换,这种情况下你可以把生成的tikz和pdf图形删掉;
  2. 中文问题:tikzDevice包好像还没有完全解决中文问题,作者正在尝试我提供的中文图形,也许在未来不久的新版本中会有完全的中文支持,但目前实际上我已经在我的自动配置脚本中解决了这个问题。我们只需要保证我们的中文LyX文档用UTF-8编码就可以了——在文档选项中使用CTeX类(以及UTF8选项),并设置语言为中文、编码为Unicode (XeLaTeX) (utf8)。这种设置至少我在Windows和MikTeX下可以成功编译中文Sweave文档。目前高亮代码功能在中文文档中还不能用,但我已经向相关作者发了补丁,等CRAN管理员回来发布新版parser包之后这个问题也会解决。所以,中文用户也不必担心,一切都(将)可以使用。
  3. 依旧是样式问题:细心的读者也许能发现上面的PDF文档中的图形和这里网页中的图形略有区别,PDF文档中图形的字体和正文字体是一致的,这也是pgf的优势,前文第5节曾提到过。
  4. 如果使用缓存,那么代码段中的对象都无法打印出来,即使用print(dat)也无法打印dat,这是因为缓存的代码段只管是否加载对象,而无法执行进一步的操作(包括打印),对于这种情况,我们应该设计好我们的Sweave文档结构,让某些代码段作计算(并缓存),某些代码段作输出(不要缓存)。

这两篇文章都有点长,但最终的操作都极其简单。如果能理解这套工具,我想无论是写作业还是写报告,都将事半功倍(比如我就一直用它写作业)。

祝大家在新的一年能写出漂亮、利索的统计文档,告别复制粘贴的体力劳动。

也谈提高R语言的运算效率

用过底层语言做计算的人转入R语言的时候一般都会觉得R语言的运算太慢,这是一个常见的对R的误解或者对R的设计的不理解。在二三十年前Chambers等一群人在贝尔实验室设计S语言之前,统计计算主要是通过那些底层语言实现的,典型的如Fortran。当时有一个基于Fortran的统计分析库,用它的麻烦就在于无论做什么样的数据分析,都涉及到一大摞繁琐的底层代码,这让数据分析变得很没劲,统计学家有统计学家的事情,天天陷在计算机程序代码中也不是个办法,要摆脱底层语言,那就只能设计高层语言了。有所得必有所失,众所周知,高层语言一般来说比底层语言低效,但对用户来说更友好。举个简单的例子,用C语言计算均值时,我们得对一个向量(数组)从头循环到尾把每个值累加起来,得到总和,然后除以向量的长度,而均值在统计当中简直是再家常便饭不过了,要是每次都要来这么个循环,大家也都甭干活儿了,天天写循环好了。

前两天COS论坛上有个帖子提到“R语言的执行效率问题”,大意如下:

有120000行数据,用perl写成12万行R命令做t.test,然后执行,大概几分钟就算完了。如果只用R语言,把所有数据先读入,然后用循环,每一行做t.test,花了几个小时,到现在还没有算完。这说明一个问题,在R里执行单行命令要比用循环快,涉及到循环的问题,最好写成单行命令执行。

我不清楚作者是如何写这“12万行”R命令的,本文假设是把t.test(dat[i, ]), i = 1, 2, ..., 120000写入一个文件,然后用source()执行之。面对这样一个问题,我们有若干种改进计算的可能性。首先我们看“硬”写入程序代码的方法:

## 为了使计算可重复,设定随机数种子
set.seed(123)
## 生成数据,随机数,120000行 x 100列矩阵
dat = matrix(rnorm(120000 * 100), 120000)
nr = nrow(dat)
nc = ncol(dat)
## 六种方法的P值向量
p1 = p2 = p3 = p4 = p5 = p6 = numeric(nr)

## via source()
f = file("test.t")
writeLines(sprintf("p1[%d] = t.test(dat[%d, ])$p.value",
    seq(nr), seq(nr)), f)
system.time({
    source(f)
})
#   user  system elapsed
#  95.36    0.19   95.86
close(f)
unlink("test.t")

1、向量式计算:apply()以及*apply()

当我们需要对矩阵的行或者列逐一计算时,apply()系列函数可能会提高效率。本例是对矩阵的行做t检验,那么可以这样计算:

## via apply()
system.time({
    p2 = apply(dat, 1, function(x) {
        t.test(x)$p.value
    })
})
#   user  system elapsed
#  63.12    0.06   63.50
identical(p1, p2)
# [1] TRUE

结果比第一种方法快了大约半分钟,而且计算结果完全一样。apply()本质上仍然是循环,但它在某些情况下比直接用显式循环要快:

## via for-loop
system.time({
    for (i in seq(nr)) p3[i] = t.test(dat[i, ])$p.value
})
#   user  system elapsed
#  69.88    0.03   70.05
identical(p2, p3)
# [1] TRUE

不过apply()系列函数在提高运算速度上优势并不会太明显,提倡用它的原因是它和统计中的矩阵运算相似,可以简化代码,相比起$\sum_{i=1}^n x_i/n$,我们可能更愿意看$\bar{x}$这样的表达式。很多R内置函数都是用底层语言写的,我们需要做的就是把一个对象作为整体传给函数去做计算,而不要自行把对象分解为一个个小部分计算,这个例子可能更能体现向量式计算的思想:

system.time(apply(dat, 1, mean))
#   user  system elapsed
#   5.28    0.04    5.25
system.time({
    x = numeric(nr)
    for (i in 1:nr) x[i] = mean(dat[i, ])
})
#   user  system elapsed
#   4.44    0.02    4.42
system.time(rowMeans(dat))
#   user  system elapsed
#   0.11    0.00    0.13

2、明确计算的目的

很多情况下,R函数返回的不仅仅是一个数字作为结果,而是会得到一系列诸如统计量、P值、各种系数等对象,在我们调用R函数之前如果能想清楚我们究竟需要什么,也许对计算的速度提升有帮助。比如本例中,也许我们仅需要12万个双边P值,其它数字我们都用不着,那么可以考虑“手工”计算P值:

## "hand" calculation in R
system.time({
    p4 = 2 * pt(apply(dat, 1, function(x, mu = 0) -abs((mean(x) -
        mu)/sqrt(var(x)/nc))), nc - 1)
})
#   user  system elapsed
#  15.97    0.07   16.08
identical(p3, p4)
# [1] TRUE

上面的计算更“纯净”,所以计算速度有了本质的提升,而且计算的结果和前面毫无差异。在做计算之前,人的脑子多思考一分钟,也许计算机的“脑子”会少转一个小时。

3、把四则运算交给底层语言

R是高层语言,把它拿来做简单的四则运算是很不划算的,而且容易导致程序低效。加加减减的活儿是C和Fortran等底层语言的强项,所以可以交给它们去做。以下我们用一段C程序来计算t统计量,然后用R CMD SHLIB将它编译为dll(Windows)或so(Linux)文件,并加载到R中,用.C()调用,最终用R函数pt()计算P值:

## "hand" calculation in C for t-statistic
writeLines("#include <math.h>

void calc_tstat(double *x, int *nr, int *nc, double *mu, double *tstat)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        tstat[i] = (mean - *mu) / sqrt(var / (*nc - 1));
    }
}", "calc_tstat.c")
system("R CMD SHLIB calc_tstat.c")
dyn.load(sprintf("calc_tstat%s", .Platform$dynlib.ext))
system.time({
    p5 = 2 * pt(-abs(.C("calc_tstat", as.double(dat), nr, nc,
        0, double(nrow(dat)))[[5]]), nc - 1)
})
#   user  system elapsed
#   0.86    0.06    0.92
dyn.unload(sprintf("calc_tstat%s", .Platform$dynlib.ext))

因为R可以用system()调用系统命令,所以整个过程全都可以用R完成,Windows用户需要安装Rtools并设置系统环境变量PATH才能使用R CMD SHLIB

更进一步,因为R自身的一些C程序也是可供用户的C程序调用的,所以我们可以把整个P值的计算过程全都扔进C代码中,一步完成:

## "hand" calculation in C calling Rmath.h
writeLines("#include <Rmath.h>
void calc_pvalue(double *x, int *nr, int *nc, double *mu, double *pval)
{
    int i, j;
    double sum = 0.0, sum2 = 0.0, mean, var;
    for (i = 0; i < *nr; i++) {
        for (j = 0; j < *nc; j++) {
            sum += x[i + j * *nr];
        }
        mean = sum / (double) *nc;
        sum = 0.0;
        for (j = 0; j < *nc; j++) {
            sum2 += (x[i + j * *nr] - mean) * (x[i + j * *nr] - mean);
        }
        var = sum2 / (double) (*nc - 1);
        sum2 = 0.0;
        pval[i] = 2 * pt(-fabs((mean - *mu) / sqrt(var / (*nc - 1))),
                      (double) (*nc - 1), 1, 0);
    }
}", "calc_pvalue.c")
system("R CMD SHLIB calc_pvalue.c")
dyn.load(sprintf("calc_pvalue%s", .Platform$dynlib.ext))
system.time({
    p6 = .C("calc_pvalue", as.double(dat), nr, nc, as.double(0),
        double(nrow(dat)))[[5]]
})
#   user  system elapsed
#   0.83    0.07    0.91
dyn.unload(sprintf("calc_pvalue%s", .Platform$dynlib.ext))

头文件Rmath.h的引入使得我们可以调用很多基于C程序的R函数,详情参考手册Writing R Extensions。通过C计算出来的P值和前面用R算的略有差异,下面画出p6 - p1 vs p1以及p6 - p5 vs p5的图:

P值的差异
P值的差异

导致差异的原因此处不细究,感兴趣的读者可以帮忙检查一下。

小结

  1. 若你熟悉底层语言,计算又不太复杂,那么可用底层语言写,然后用R调之;
  2. 否则把R对象当做整体去计算,能做x + 1就不要做for (i in length(x)) x[i] + 1
  3. 不要低估R core们的编程水平,他们已经做了很多工作让用户脱离底层编程

注:本文中的运算时间可能不可重复,这与计算机所处的状态有关,但大体来说,运算速度的快慢是可以肯定的。本文仅仅是关于统计计算的一个微小的例子,以后若有更好的例子,可能会更新本文;也欢迎各位提供更多示例。