R高级编程技巧及Rcpp的介绍.PPT

1、R高级编程技巧及Rcpp的介绍颜林林北京大学生物信息中心2011年5月28日2R语言的特点FREE：开源、免费、自由、灵活 a nrow(a)截至2011年5月28日，CRAN上总计可用软件包数：30083遵守各种不同开源协议：R语言的特点FREE：开源、免费、自由、灵活 sort(table(gsub( .*, , a,License), decreasing = T)4高校、科研机构、企事业单位金融、通讯、医药、生物、环境R语言的广泛运用5COS论坛创办于2006年5月19日截至2011年5月28日：整个论坛：主题：16,000+；回复：109,000+“S-Plus & R语言”

2、版：主题：4,000+；回复：27,000+学习R的一个捷径：在在COS上回帖！上回帖！6论坛中一些常见关于R的问题如何画XX图？如何实现XX功能？安装时遇到XX问题如何解决？如何提高R程序的运行效率？153分钟学会R7正题：如何让R飞起来R高级编程技巧向量运算软件包的使用Rcpp的介绍结合不同语言的优势拓展R的运用范围8向量运算“是一种习惯”循环的思维方式：for (p in 参加R大会的人)发资料给p;发胸牌给p;发给p;生活中的思维方式：发资料、胸牌等给参加R大会的人能用更少的话说清楚，就别罗嗦！ a a * 2 sin(a)9R是一门解释型语言 a # 向量运算向量运算 b # 循环循

3、环 b for (i in seq_along(a) bi system.time(replicate(10000, b b system.time(replicate(10000, for (i in seq_along(a) bi 3, ac(x,y), .数据类型之间的互相转换as.matrix(), as.vector(), as.character(), .apply系列函数apply(), sapply(), lapply(), tapply(), .其他各种函数sum(), mean(), cumsum(), combn(), .12R语言元素类型整数（integer）实数（nu

4、meric, double）复数（complex）字符（character）逻辑（logical）原始数据（raw）因子（factor）1:10c(1.1, 3.14, 10)c(1+i, 3-2i)c(a, b, COS)c(TRUE, FALSE)as.raw(48)as.factor(letter1:3)13R语言基本数据类型向量（vector）列表（list）矩阵（matrix）数组（array）数据框（data.frame）14R语言基本数据类型向量（vector）列表（list）矩阵（matrix）数组（array）数据框（data.frame）一维，同类型元素一维，不同类型元素二

5、维，同类型元素二维或多维，同类型元素二维，每列元素同类型15学习R函数时应注意16apply系列函数17函数apply (m apply(m, 1, sum) # 相当于相当于 rowSums(m)1 16 20 apply(m, 2, sum) # 相当于相当于 colSums(m)1 3 7 11 15 apply(m, 1:2, function(x) x / 2) ,1 ,2 ,3 ,41, 0.5 1.5 2.5 3.52, 1.0 2.0 3.0 4.018函数lapply和sapply (l lapply(l, sum) # 返回列表返回列表$a1 15$b1 40 sapply

6、(l, sum) # 返回向量返回向量 a b 15 4019函数sapply和vapply (X sapply(X, fivenum) ,1 ,2 ,31, 1.0 1.0 12, 1.5 1.5 23, 2.0 2.5 34, 2.5 3.5 45, 3.0 4.0 520函数sapply和vapply (v vapply(X, fivenum, v) ,1 ,2 ,3Min. 1.0 1.0 11st Qu. 1.5 1.5 2Median 2.0 2.5 33rd Qu. 2.5 3.5 4Max. 3.0 4.0 521函数tapply (d c(class(d), class(d$

7、age), class(d$gender)1 data.frame numeric factor tapply(d$age, d$gender, mean)Female Male 24.0 24.522函数mapply mapply(function(x, y) seq_len(x) + y, c(a = 1, b = 2, c = 3), # names from first c(A = 10, B = 0, C = -10)$a1 11$b1 1 2$c1 -9 -8 -723函数rapply (X rapply(X, function(x) x) a b.c d 3.1415926535

8、8979 1 a test 24举例：求DNA互补链DNA由A、C、T、G组成DNA两条链方向相反成对碱基，A与T配对，C与G配对已知其中一条为“ACTGAAGTGC”求另一条序列25举例：求DNA互补链# 方法一：循环方法一：循环revcom1 - function(DNA) result - for (i in nchar(DNA):1) n - substr(DNA, i, i)if (n = A) result - paste(result, T, sep = )if (n = T) result - paste(result, A, sep = )if (n = C) result

9、- paste(result, G, sep = )if (n = G) result - paste(result, C, sep = )result26举例：求DNA互补链# 方法二：方法二：*apply()revcom2 - function(DNA) s - strsplit(DNA, )1s - sapply(s, function(x) if (x = A) return(T)if (x = T) return(A)if (x = G) return(C)if (x = C) return(G)paste(rev(s), collapse = )27# 方法三：方法三：names(

10、)revcom3 - function(DNA) tr - c(A,C,G,T)names(tr) tr T G C A A C G T 举例：求DNA互补链28方法比较 library(rbenchmark) benchmark(revcom1(DNA), revcom2(DNA), revcom3(DNA), columns = c(test, replications, elapsed,relative), order = relative,replications = 1000) test replications elapsed relative3 revcom3(DNA) 1000

11、 0.163 1.0000002 revcom2(DNA) 1000 0.597 3.6625771 revcom1(DNA) 1000 0.624 3.82822129R的灵活性R是什么语言写成的？C/C+、Fortran、RR的外部扩展Writing R Extensions30R与各种语言R + bash / RscriptRDCOMrJavarpy, rpy2EmbededRserve、R API、Rcpp31R 与 C+R（解释型语言）无需编译不需要其它文件弱类型语言灵活性好运行速度较慢其他特性向量运算大量统计和绘图函数C+（编译型语言）需要事先编译需要头文件和库文件强类型语言灵活性

12、相对较差运行速度快其他特性面向对象模板与泛型编程32从C到C+字符串、数组等类型内存的管理模板和泛型编程STL (Standard Template Library)33Rcpp相关历史RQuantLibRcppTemplateRservercppbind / RAbstraction / RObjects / CXXRR APIRcpp / RInside / RProtoBuf34Rcpp的运用在R中调用C+模块（Rcpp）用C+实现计算功能通过C+调用其他软件库在C+中使用R（RInside）向量运算（STL）使用R的统计函数使用R的绘图函数35举例：Rcpp和inline带来的速度f

13、- function(n, x=1) for (i in 1:n) x=1/(1+x)g - function(n, x=1) for (i in 1:n) x=(1/(1+x)h - function(n, x=1) for (i in 1:n) x=(1+x)(-1)j - function(n, x=1) for (i in 1:n) x=1/1+xk - function(n, x=1) for (i in 1:n) x=1/1+xlibrary(rbenchmark)N - 1e5benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1),col

14、umns=c(test, replications, elapsed, relative), order=relative, replications=10)36举例：Rcpp和inline带来的速度f - function(n, x=1) for (i in 1:n) x=1/(1+x)g - function(n, x=1) for (i in 1:n) x=(1/(1+x)h - function(n, x=1) for (i in 1:n) x=(1+x)(-1)j - function(n, x=1) for (i in 1:n) x=1/1+xk - function(n, x=1

15、) for (i in 1:n) x=1/1+x test replications elapsed relative5 k(N, 1) 10 6.103 1.0000001 f(N, 1) 10 6.426 1.0529254 j(N, 1) 10 6.835 1.1199412 g(N, 1) 10 7.339 1.2025233 h(N, 1) 10 9.226 1.51171637举例：Rcpp和inline带来的速度library(inline)src - int n = as(ns);double x = as(xs);for (int i=0; in; i+) x=1/(1+x)

16、;return wrap(x); l benchmark( f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), l(N,1), columns=c(test, replications, elapsed, relative), order=relative, eplications=10) test replications elapsed relative6 l(N, 1) 10 0.027 1.00005 k(N, 1) 10 6.190 229.25931 f(N, 1) 10 6.379 236.25934 j(N, 1) 10 6.808 252.14812 g(N, 1) 10 7.267 269.14813 h(N, 1) 10 9.184 340.148139举例：Rcpp API#include RcppExport SEXP foo(SEXP a, SEXP b) Rcpp:NumericVector xa(a);Rcpp:NumericVector xb(b);int n = std:m