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

R高级编程技巧及Rcpp的介绍颜林林北京大学生物信息中心2011年5月28日R语言的特点—FREE：开源、免费、自由、灵活> a <- available.packages(contrib.url("http://ftp.ctex.org/mirrors/CRAN", "source"))> nrow(a)截至2011年5月28日，CRAN上总计可用软件包数：3008遵守各种不同开源协议：R语言的特点GPLGPL-2￿GPL-3￿LGPL...185768013165...—FREE：开源、免费、自由、灵活> sort(table(gsub(" .*", "", a[,"License"])), decreasing = T)—高校、科研机构、企事业单位……—金融、通讯、医药、生物、环境……R语言的广泛运用COS论坛—创办于2006年5月19日—截至2011年5月28日：—整个论坛：—主题：16,000+；回复：109,000+—“S-Plus & R语言”版：—主题：4,000+；回复：27,000+—学习R的一个捷径：—在COS上回帖！论坛中一些常见关于R的问题—如何画XX图？—如何实现XX功能？—安装时遇到XX问题如何解决？—如何提高R程序的运行效率？—……《153分钟学会R》正题：如何让R飞起来—R高级编程技巧—向量运算—软件包的使用—Rcpp的介绍—结合不同语言的优势—拓展R的运用范围向量运算“是一种习惯”—循环的思维方式：for (p in 参加R大会的人){发资料给p;发胸牌给p;发……给p;}—生活中的思维方式：发资料、胸牌等给参加R大会的人能用更少的话说清楚，就别罗嗦！> a <- 1:100> a * 2> sin(a)R是一门解释型语言> a <- 1:100> # 向量向量运算运算> b <- sin(a)> # 循循环环> b <- vector()> for (i in seq_along(a)) {b[i] <- sin(a[i])}R是一门解释型语言> system.time(replicate(10000, b <- sin(a))) user system elapsed 0.44 0.10 0.59> b <- vector()> system.time(replicate(10000, for (i in seq_along(a)) {b[i] <- sin(a[i])})) user system elapsed 6.15 0.02 6.44把苦力活交给底层，尽量避免使用循环！如何掌握好向量运算—R语言基本数据类型—vector, matrix, array, data.frame, list, ...—元素的提取方法—a[1:3], a[-2], a[a>3], a[c("x","y")], ...—数据类型之间的互相转换—as.matrix(), as.vector(), as.character(), ...—apply系列函数—apply(), sapply(), lapply(), tapply(), ...—其他各种函数—sum(), mean(), cumsum(), combn(), ...R语言元素类型—整数（integer）—实数（numeric, 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(letter[1:3])R语言基本数据类型—向量（vector）—列表（list）—矩阵（matrix）—数组（array）—数据框（data.frame）R语言基本数据类型—向量（vector）—列表（list）—矩阵（matrix）—数组（array）—数据框（data.frame）一维，同类型元素一维，不同类型元素二维，同类型元素二维或多维，同类型元素二维，每列元素同类型学习R函数时应注意apply系列函数apply操作矩阵或数组，返回数组、矩阵或向量lapply操作向量或列表，返回列表sapply操作向量或列表，返回向量或矩阵tapply操作向量，根据因子分组，返回向量mapply操作多个向量，返回列表vapply与sapply同，操作列表，返回向量或列表rapply与lapply类似，操作向量与列表，返回向量函数apply> (m <- matrix(1:8, 2, 4)) [,1] [,2] [,3] [,4][1,] 1 3 5 7[2,] 2 4 6 8> 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] [,4][1,] 0.5 1.5 2.5 3.5[2,] 1.0 2.0 3.0 4.0函数lapply和sapply> (l <- list(a = 1:5, b = 6:10))$a[1] 1 2 3 4 5$b[1] 6 7 8 9 10> lapply(l, sum) # 返回返回列表列表$a[1] 15$b[1] 40> sapply(l, sum) # 返回返回向量向量 a b 15 40函数sapply和vapply> (X <- sapply(3:5, seq))[[1]][1] 1 2 3[[2]][1] 1 2 3 4[[3]][1] 1 2 3 4 5> sapply(X, fivenum) [,1] [,2] [,3][1,] 1.0 1.0 1[2,] 1.5 1.5 2[3,] 2.0 2.5 3[4,] 2.5 3.5 4[5,] 3.0 4.0 5函数sapply和vapply> (v <- c("Min."=0, "1st Qu."=0, "Median"=0, "3rd Qu."=0, "Max."=0)) Min. 1st Qu. Median 3rd Qu. Max. 0 0 0 0 0 > vapply(X, fivenum, v) [,1] [,2] [,3]Min. 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 5函数tapply> (d <- data.frame(name = c("Alice", "Bob", "Kate", "Jim"),age = c(20, 25, 28, 24),gender = c("Female", "Male", "Female", "Male"))) name age gender1 Alice 20 Female2 Bob 25 Male3 Kate 28 Female4 Jim 24 Male> c(class(d), class(d$age), class(d$gender))[1] "data.frame" "numeric" "factor" > tapply(d$age, d$gender, mean)Female Male 24.0 24.5函数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))$a[1] 11$b[1] 1 2$c[1] -9 -8 -7函数rapply> (X <- list(list(a=pi, b=list(c=1:1)), d="a test"))[[1]][[1]]$a[1] 3.141593[[1]]$b[[1]]$b$c[1] 1$d[1] "a test"> rapply(X, function(x) x) a b.c d "3.14159265358979" "1" "a test" 举例：求DNA互补链—DNA由A、C、T、G组成DNA两条链方向相反成对碱基，A与T配对，C与G配对—已知其中一条为“ACTGAAGTGC”求另一条序列举例：求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 <- paste(result, 'G', sep = '')if (n == 'G') result <- paste(result, 'C', sep = '')}result}举例：求DNA互补链# 方法二：方法二：*apply()revcom2 <- function(DNA) {s <- strsplit(DNA, '')[[1]]s <- 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 = '')}# 方法三：方法三：names()revcom3 <- function(DNA) {tr <- c('A','C','G','T')names(tr) <- rev(tr)paste(tr[rev(strsplit(DNA, '')[[1]])],collapse = '')}> tr T G C A "A" "C" "G" "T" 举例：求DNA互补链方法比较> 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 0.163 1.0000002 revcom2(DNA) 1000 0.597 3.6625771 revcom1(DNA) 1000 0.624 3.828221R的灵活性—R是什么语言写成的？—C/C++、Fortran、R—R的外部扩展—Writing R ExtensionsR与各种语言—R + bash / Rscript—RDCOM—rJava—rpy, rpy2—Embeded—Rserve、R API、RcppR 与 C++R（解释型语言）—无需编译—不需要其它文件—弱类型语言—灵活性好—运行速度较慢其他特性—向量运算—大量统计和绘图函数C++（编译型语言）—需要事先编译—需要头文件和库文件—强类型语言—灵活性相对较差—运行速度快其他特性—面向对象—模板与泛型编程从C到C++—字符串、数组等类型—内存的管理—模板和泛型编程—STL (Standard Template Library)Rcpp相关历史—RQuantLib—RcppTemplate—Rserve—rcppbind / RAbstraction / RObjects / CXXR—R API—Rcpp / RInside / RProtoBufRcpp的运用—在R中调用C++模块（Rcpp）—用C++实现计算功能—通过C++调用其他软件库—在C++中使用R（RInside）—向量运算（STL）—使用R的统计函数—使用R的绘图函数举例：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+x}}k <- function(n, x=1) for (i in 1:n) x=1/{1+x}library(rbenchmark)N <- 1e5benchmark(f(N,1), g(N,1), h(N,1), j(N,1), k(N,1),columns=c("test", "replications", "elapsed", "relative"), order="relative", replications=10)举例：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+x}}k <- function(n, x=1) 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.511716举例：Rcpp和inline带来的速度library(inline)src <- 'int n = as(ns);double x = as(xs);for (int i=0; i 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.1481举例：Rcpp API#include RcppExport SEXP foo(SEXP a, SEXP b) {Rcpp::NumericVector xa(a);Rcpp::NumericVector xb(b);int n = std::min(a.size(), b.size());Rcpp::NumericVector xab(n);for (int i = 0; i < n; ++i) {xab[i] = xa[i] * xb[i];}return xab;}举例：RInside#include int main(int argc, char *argv[]) {RInside R(argc, argv);R["txt"] = "Hello, world!\n";R.parseEvalQ("cat(txt)");exit(0);}参考—R-introduction—R-manual—谢 谢yanlinlin82@。