1、R高級(jí)編程技巧及Rcpp的介紹,顏林林 北京大學(xué)生物信息中心 2011年5月28日,R語言的特點(diǎn),FREE：開源、免費(fèi)、自由、靈活, a nrow(a),截至2011年5月28日，CRAN上總計(jì)可用軟件包數(shù)：,3008,遵守各種不同開源協(xié)議：,R語言的特點(diǎn),FREE：開源、免費(fèi)、自由、靈活, sort(table(gsub( .*, , a,License), decreasing = T),高校、科研機(jī)構(gòu)、企事業(yè)單位 金融、通訊、醫(yī)藥、生物、環(huán)境,R語言的廣泛運(yùn)用,COS論壇,創(chuàng)辦于2006年5月19日 截至2011年5月28日： 整個(gè)論壇： 主題：16,000+；回復(fù)：109,000+ “

2、S-Plus 發(fā)胸牌給p; 發(fā)給p; ,生活中的思維方式： 發(fā)資料、胸牌等 給參加R大會(huì)的人,能用更少的話說清楚，就別羅嗦！, a a * 2 sin(a),R是一門解釋型語言, a # 向量運(yùn)算 b # 循環(huán) b for (i in seq_along(a) bi - sin(ai) ,R是一門解釋型語言, system.time(replicate(10000, b b system.time(replicate(10000, for (i in seq_along(a) bi - sin(ai) user system elapsed 6.15 0.02 6.44,把苦力活交給底層，盡量

3、避免使用循環(huán)！,如何掌握好向量運(yùn)算,R語言基本數(shù)據(jù)類型 vector, matrix, array, data.frame, list, . 元素的提取方法 a1:3, a-2, aa3, ac(x,y), . 數(shù)據(jù)類型之間的互相轉(zhuǎn)換 as.matrix(), as.vector(), as.character(), . apply系列函數(shù) apply(), sapply(), lapply(), tapply(), . 其他各種函數(shù) sum(), mean(), cumsum(), combn(), .,R語言元素類型,整數(shù)（integer） 實(shí)數(shù)（numeric, double） 復(fù)數(shù)（c

4、omplex） 字符（character） 邏輯（logical） 原始數(shù)據(jù)（raw） 因子（factor）,1:10 c(1.1, 3.14, 10) c(1+i, 3-2i) c(a, b, COS) c(TRUE, FALSE) as.raw(48) as.factor(letter1:3),R語言基本數(shù)據(jù)類型,向量（vector） 列表（list） 矩陣（matrix） 數(shù)組（array） 數(shù)據(jù)框（data.frame）,R語言基本數(shù)據(jù)類型,向量（vector） 列表（list） 矩陣（matrix） 數(shù)組（array） 數(shù)據(jù)框（data.frame）,一維，同類型元素 一維，不同類型

5、元素 二維，同類型元素 二維或多維，同類型元素 二維，每列元素同類型,學(xué)習(xí)R函數(shù)時(shí)應(yīng)注意,apply系列函數(shù),函數(shù)apply, (m apply(m, 1, sum) # 相當(dāng)于 rowSums(m) 1 16 20 apply(m, 2, sum) # 相當(dāng)于 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,函數(shù)lapply和sapply, (l lapply(l, sum) # 返回列表 $a 1 15 $b 1 40 sapp

6、ly(l, sum) # 返回向量 a b 15 40,函數(shù)sapply和vapply, (X 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,函數(shù)sapply和vapply, (v vapply(X, fivenum, v) ,1 ,2 ,3 Min. 1.0 1.0 1 1st Qu. 1.5 1.5 2 Median 2.0 2.5 3 3rd Qu. 2.5 3.5 4 Max. 3.0 4.0 5,函數(shù)tapply, (d c(class(d),

7、class(d$age), class(d$gender) 1 data.frame numeric factor tapply(d$age, d$gender, mean) Female Male 24.0 24.5,函數(shù)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,函數(shù)rapply, (X rapply(X, function(x) x) a b.

8、c d 3.14159265358979 1 a test,舉例：求DNA互補(bǔ)鏈,DNA由A、C、T、G組成DNA兩條鏈方向相反成對(duì)堿基，A與T配對(duì)，C與G配對(duì) 已知其中一條為“ACTGAAGTGC”求另一條序列,舉例：求DNA互補(bǔ)鏈,# 方法一：循環(huán) 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 (

9、n = C) result - paste(result, G, sep = ) if (n = G) result - paste(result, C, sep = ) result ,舉例：求DNA互補(bǔ)鏈,# 方法二：*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

10、= ) ,# 方法三：names() revcom3 - function(DNA) tr - c(A,C,G,T) names(tr) - rev(tr) paste(trrev(strsplit(DNA, )1), collapse = ) , tr T G C A A C G T,舉例：求DNA互補(bǔ)鏈,方法比較, library(rbenchmark) benchmark(revcom1(DNA), revcom2(DNA), revcom3(DNA), columns = c(test, replications, elapsed,relative), order = relative

11、,replications = 1000) test replications elapsed relative 3 revcom3(DNA) 1000 0.163 1.000000 2 revcom2(DNA) 1000 0.597 3.662577 1 revcom1(DNA) 1000 0.624 3.828221,R的靈活性,R是什么語言寫成的？ C/C+、Fortran、R R的外部擴(kuò)展 Writing R Extensions,R與各種語言,R + bash / Rscript RDCOM rJava rpy, rpy2 Embeded Rserve、R API、Rcpp,R 與

12、C+,R（解釋型語言） 無需編譯 不需要其它文件 弱類型語言 靈活性好 運(yùn)行速度較慢 其他特性 向量運(yùn)算 大量統(tǒng)計(jì)和繪圖函數(shù),C+（編譯型語言） 需要事先編譯 需要頭文件和庫文件 強(qiáng)類型語言 靈活性相對(duì)較差 運(yùn)行速度快 其他特性 面向?qū)ο?模板與泛型編程,從C到C+,字符串、數(shù)組等類型 內(nèi)存的管理 模板和泛型編程 STL (Standard Template Library),Rcpp相關(guān)歷史,RQuantLib RcppTemplate Rserve rcppbind / RAbstraction / RObjects / CXXR R API Rcpp / RInside / RProto

13、Buf,Rcpp的運(yùn)用,在R中調(diào)用C+模塊（Rcpp） 用C+實(shí)現(xiàn)計(jì)算功能 通過C+調(diào)用其他軟件庫 在C+中使用R（RInside） 向量運(yùn)算（STL） 使用R的統(tǒng)計(jì)函數(shù) 使用R的繪圖函數(shù),舉例：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

14、k - function(n, x=1) for (i in 1:n) x=1/1+x library(rbenchmark) N - 1e5 benchmark(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 i

15、n 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 relative 5 k(N, 1) 10 6.103 1.000000 1 f(N, 1) 10 6.426 1.052925 4 j(N, 1) 10 6.835 1.119941 2 g(N, 1) 10 7.339 1.202523

16、 3 h(N, 1) 10 9.226 1.511716,舉例：Rcpp和inline帶來的速度,library(inline) src (ns); double x = as(xs); for (int i=0; in; i+) x=1/(1+x); return wrap(x); l - cxxfunction(signature(ns=integer, xs=numeric), body=src, plugin=Rcpp),舉例：Rcpp和inline帶來的速度, benchmark( f(N,1), g(N,1), h(N,1), j(N,1), k(N,1), l(N,1), col

17、umns=c(test, replications, elapsed, relative), order=relative, eplications=10) test replications elapsed relative 6 l(N, 1) 10 0.027 1.0000 5 k(N, 1) 10 6.190 229.2593 1 f(N, 1) 10 6.379 236.2593 4 j(N, 1) 10 6.808 252.1481 2 g(N, 1) 10 7.267 269.1481 3 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(),