当前位置:天才代写 > tutorial > 其他教程 > R语言的执行效率问题

R语言的执行效率问题

2017-12-04 08:00 星期一 所属: 其他教程 浏览:418

用过底层语言做计较的人转入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()执行之。面临这样一个问题,我们有若干种改造计较的大概性。首先我们看“硬”写入措施代码的要领:








1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20



## 为了使计较可反复,设定随机数种子


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检讨,那么可以这样计较:








1


2


3


4


5


6


7


8


9


10



## 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()本质上仍然是轮回,但它在某些环境下比直接用显式轮回要快:








1


2


3


4


5


6


7


8



## 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内置函数都是用底层语言写的,我们需要做的就是把一个工具作为整体传给函数去做计较,而不要自行把工具解析为一个个小部门计较,这个例子大概更能浮现向量式计较的思想:








1


2


3


4


5


6


7


8


9


10


11


12



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值:








1


2


3


4


5


6


7


8


9



## “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值:








1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30



## “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代码中,一步完成:








1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30



## “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们的编程程度,他们已经做了许多事情让用户离开底层编程

注:本文中的运算时间大概不行反复,这与计较机所处的状态有关,但概略来说,运算速度的快慢是可以必定的。本文仅仅是关于统计计较的一个微小的例子,今后若有更好的例子,大概会更新本文;也接待列位提供更多示例。



摘自统计之都 : http://cos.name/2009/12/improve-r-computation-efficiency/

 

    关键字:

天才代写-代写联系方式