|
|
相关文章推荐
|
大方的水煮鱼 · 2024商务区民生实事答卷· 3 月前 · |
|
|
含蓄的瀑布 · 全芳備祖集 (四庫全書本) - ...· 3 月前 · |
|
|
发财的楼房 · 心理健康 | ...· 4 月前 · |
|
|
瘦瘦的小熊猫 · 驻爱丁堡总领馆主办“情系汶川”大型赈灾义演_ ...· 4 月前 · |
|
|
强健的吐司 · 基于python+opencv+mediap ...· 1 年前 · |
evalCpp()
转换单一计算表达式
cppFunction()
转换简单的C++函数—Fibnacci例子
sourceCpp()
转换C++程序—正负交替迭代例子
sourceCpp()
转换C++源文件中的程序—正负交替迭代例子
sourceCpp()
转换C++源程序文件—卷积例子
wrap()
把C++变量返回到R中
as()
函数把R变量转换为C++类型
as()
和
wrap()
的隐含调用
//[[Rcpp::export]]
sourceCpp()
函数中直接包含C++源程序字符串
cppFunction()
函数中直接包含C++函数源程序字符串
evalCpp()
函数中直接包含C++源程序表达式字符串
depends
指定要链接的库
invisible
要求函数结果不自动显示
clone
函数
is_na
seq_along
seq_len
pmin
和
pmax
ifelse
sapply
和
lapply
sign
diff
kable()
函数制作表格
Rprof()
函数可以执行性能分析的数据收集工作，
收集到的性能数据用
summaryRprof()
函数可以显示运行最慢的函数。
如果使用RStudio软件，可以用
Profile
菜单执行性能数据收集与分析，
可以在图形界面中显示程序中哪些部分运行花费时间最多。
在改进已有程序的效率时， 第一要注意的就是不要把原来的正确算法改成一个速度更快但是结果错误的算法。 这个问题可以通过建立试验套装， 用原算法与新算法同时试验看结果是否一致来避免。 多种解决方案的正确性都可以这样保证， 也可以比较多种解决方案的效率。
本章后面部分描述常用的改善性能的方法。 对于涉及到大量迭代的算法， 如果用R实现性能太差不能满足要求， 可以改成C++编码，用Rcpp扩展包连接到R中。 Rcpp扩展包的使用将单独讲授。
R的运行效率也受到内存的影响， 占用内存过多的算法有可能受到物理内存大小限制无法运行， 过多复制也会影响效率。
如果要实现一个比较单纯的不需要利用R已有功能的算法， 发现用R计算速度很慢的时候， 也可以考虑先用Julia语言实现。 Julia语言设计比R更先进，运算速度快得多， 运算速度经常能与编译代码相比， 缺点是刚刚诞生几年的时间， 可用的软件包还比较少。
18.2 性能比较和向量化编程示例
充分利用R的向量化编程和内建函数可以提高程序效率。
R函数
system.time()
可以测量某个表达式的运行时间，
proc.time()
可以用作会话期间的计时时钟，
扩展包bench可以用来比较不同表达式的运行时间。
例18.1 假设要计算如下的统计量： w = \frac{1}{n} \sum_{i=1}^n | x_i - \hat m |, 其中 \(x_1, x_2, \dots, x_n\) 是某总体的样本， \(\hat m\) 是样本中位数。 比较不同的编程方法。
用传统的编程风格， 把这个统计量的计算变成一个R函数，可能会写成：
mad_f1 <- function(x){ n <- length(x) mhat <- median(x) s <- 0.0 for(i in 1:n){ s <- s + abs(x[i] - mhat) s <- s/n return(s)用R的向量化编程，函数体只需要一个表达式：
其中
x - median(x)利用了向量与标量运算结果是向量每个元素与标量运算的规则，abs(x - median(x))利用了abs()这样的一元函数如果以向量为输入就输出每个元素的函数值组成的向量的规则，mean(...)避免了求和再除以n的循环也不需要定义多余的变量n。显然，第二种做法的程序比第一种做法简洁的多， 如果多次重复调用， 第二种做法的计算速度比第一种要快几十倍甚至上百倍。 用
system.time()函数可以求某个表达式的计算时间， 返回结果的第3项是流逝时间。 下面对x采用10000个随机数， 并重复计算1000次，比较两个程序的效率：nrep <- 1000 x <- runif(10000) y1 <- numeric(nrep); y2 <- y1 system.time(for(i in 1:nrep) y1[i] <- mad_f1(x) )[3] ## elapsed ## 1.19 system.time(for(i in 1:nrep) y1[i] <- mad_f2(x) )[3] ## elapsed ## 0.22速度相差数倍。
用
system.time()或者proc.time()测量运行时间会受到后台正在运行的其它任务的干扰。 有一个R扩展包bench可以多次运行表达式， 从而比较两个或多个表达式的运行时间（类似的包还有microbenchmark）。x <- runif(10000) bench::mark( mad_f1(x), mad_f2(x) ## Unit: microseconds ## expr min lq mean median uq max neval cld ## f1(x) 1321.802 1519.151 1642.8350 1555.001 1590.5510 5826.800 100 b ## f2(x) 275.801 312.951 330.2541 322.501 326.9515 781.001 100 a因为测试运行受到操作系统中同时运行的其它任务的干扰， 所以我们不应该看平均时间， 而应该看中位数， 或者最小时间（最快可能）。 就运行时间中位数（单位：毫秒）来看，
f2()比f1()快大约4倍。例18.2 假设要编写函数计算 \begin{aligned} f(x) = \begin{cases} 1 & x \geq 0, \\ 0 & \text{其它} \end{cases} \end{aligned} 如何向量化编程？
利用传统思维，程序写成
pos_f1 <- function(x){ n <- length(x) y <- numeric(n) for(i in seq_along(x)){ if(x[i] >= 0) y[i] <- 1 else y[i] <- 0实际上，
y <- numeric(n)使得y的每个元素都初始化为0， 所以程序中else y[i] <- 0可以去掉。利用向量化与逻辑下标，程序可以写成:
pos_f2 <- function(x){ n <- length(x) y <- numeric(n) y[x >= 0] <- 1但是，利用R中内建函数
ifelse()， 可以把函数体压缩到仅用一个语句：例18.3 考虑一个班的学生存在生日相同的概率。 假设一共有365个生日(只考虑月、日)。 设一个班有n个人, 当n大于365时
{至少两个人有生日相同}是必然事件(概率等于1)。当\(n\)小于等于365时： \begin{aligned} & P(\text{至少有两人同生日}) \\ =& 1 - P(\text{$n$个人生日彼此不同}) \\ =& 1 - \frac{365\times 364\times\cdots\times(365 - (n-1))}{365^n} \\ =& 1 - \frac{365-0}{365} \cdot \frac{365-1}{365} \cdots \frac{365 - (n-1)}{365} \end{aligned}
对\(n=1,2,\dots, 365\)来计算对应的概率。
完全用循环（两重循环），程序写成：
birth_f1 <- function(){ ny <- 365 x <- numeric(ny) for(n in 1:ny){ s <- 1 for(j in 0:(n-1)){ s <- s * (365-j)/365 x[n] <- 1 - s 不能先计算\(365\times 364\times\cdots\times(365 - (n-1))\) 再和\(365^n\)相除， 这会造成数值溢出。用
prod()函数可以向量化内层循环：birth_f2 <- function(){ ny <- 365 x <- numeric(ny) for(n in 1:ny){ x[n] <- 1 - prod((365:(365-n+1))/365)程序利用了向量与标量的除法， 以及内建函数
prod()。把程序用
cumprod()函数改写， 可以完全避免循环：birth_f3 <- function(){ ny <- 365 x <- 1 - cumprod((ny:1)/ny)用bench包比较:
bench:: mark( birth_f1(), birth_f2(), birth_f3() ## Unit: microseconds ## expr min lq mean median uq max neval cld ## f1() 2467.000 2792.1515 2963.27806 2979.9515 3033.202 8058.001 100 c ## f2() 330.100 387.1005 674.42992 413.3505 438.301 8445.601 100 b ## f3() 1.901 2.3515 27.51003 2.8510 3.651 2322.500 100 a
birth_f2()比birth_f1()快大约6倍，birth_f3()比birth_f2()又快了大约140倍倍，birth_f3()比birth_f1()快了1000倍！18.3 减少显式循环
显式循环是R运行速度较慢的部分， 有循环的程序也比较冗长， 与R的向量化简洁风格不太匹配。 在循环内修改数据子集，例如数据框子集， 可能会先制作副本再修改， 这当然会损失很多效率。 R 3.1.0版本以后列表元素在修改时不制作副本， 但数据框还会制作副本。
前面已经指出， 利用R的向量化运算可以减少很多循环程序。
R中的有些运算可以用内建函数完成， 如
sum,prod,cumsum,cumprod,mean,var,sd等。 这些函数以编译程序的速度运行， 不存在效率损失。R的
sin,sqrt,log等函数都是向量化的， 可以直接对输入向量的每个元素进行变换。对矩阵，用
apply函数汇总矩阵每行或每列。colMeans,rowMeans可以计算矩阵列平均和行平均，colSums,rowSums可以计算矩阵列和与行和。apply类函数有多个， 包括apply, sapply, lapply, tapply, vapply, replicate等。 这些函数不一定能提高程序运行速度， 但是使用这些函数更符合R的程序设计风格， 使程序变得简洁， 程序更简洁并不等同于程序更容易理解， 要理解这样的程序， 需要更多学习与实践。 参见25。
18.3.1
replicate()函数
replicate()函数用来执行某段程序若干次， 类似于for()循环但是没有计数变量。 常用于随机模拟。replicate()的缺省设置会把重复结果尽可能整齐排列成一个多维数组输出。
replicate(重复次数, 要重复的表达式)其中的表达式可以是复合语句, 也可以是执行一次模拟的函数。
下面举一个简单模拟例子。 设总体\(X\)为\(\text{N}(0, 1)\), 取样本量\(n=5\)， 重复地生成模拟样本共\(B=6\)组， 输出每组样本的样本均值和样本标准差。 模拟可以用如下的
replicate()实现：## [,1] [,2] [,3] [,4] [,5] [,6] ## [1,] 0.1292699 0.1351357 0.03812297 0.4595670 0.08123054 -0.3485770 ## [2,] 0.9610394 0.6688342 1.49887443 0.4648177 1.20109623 0.7046822结果是一个矩阵，矩阵行数与每次模拟的结果（均值、标准差）个数相同， 这里第一行都是均值，第二行都是标准差； 矩阵每一列对应于一次模拟。 这样结果存储应该是来源于R的数组元素次序为列优先次序。
18.4 避免制作副本
类似于
x <- c(x, y)这样的累积结果每次运行都会制作一个x的副本， 在x存储量较大或者重复修改次数很多时会减慢程序。例18.4 执行10000次模拟， 每次模拟求10个U(0,1)均匀随机数的极差， 保存每次的极差， 求平均极差。
用
x <- c(x, y)这样的写法：set.seed(101) system.time({ M <- 1E5 x <- c() for(i in seq(M)){ x <- c(x, diff(range(runif(10)))) mean(x) ## 用户 系统 流逝 ## 22.81 2.58 25.49上面的程序不仅是低效率的做法， 也没有预先精心设计用来保存结果的数据结构。 数据建模或随机模拟的程序都应该事先分配好用来保存结果的数据结构， 在每次循环中填入相应结果。如：
set.seed(101) system.time({ M <- 1E5 x <- numeric(M) for(i in seq(M)){ x[[i]] <- diff(range(runif(10))) mean(x) ## 用户 系统 流逝 ## 0.96 0.01 0.96这样的程序结构更清晰， 效率更高， 而且循环次数越多， 比
x <- c(x, ...)这样的做法的优势越大。例18.5 在循环内修改数据框的值也会制作数据框副本， 当数据框很大或者循环次数很多时会使得程序很慢。 下面给出示例。
set.seed(101) m <- 2E4; n <- 100 x <- as.data.frame(matrix( runif(n*m), nrow=n, ncol=m)) system.time({ for(j in seq(m)){ x[[j]] <- x[[j]] + 1 ## 用户 系统 流逝 ## 11.53 0.03 11.55在循环内修改列表元素就不会制作副本， 从而可以避免这样的效率问题，如：
set.seed(101) m <- 2E4; n <- 100 x <- replicate(m, runif(n), simplify=FALSE) system.time({ for(j in seq(m)){ x[[j]] <- x[[j]] + 1 ## 用户 系统 流逝 ## 0.02 0.02 0.04 x <- as.data.frame(x)
replicate()函数中用simplify=FALSE使结果总是返回列表。 要注意的是， 上面第二个程序中的as.data.frame(x)也是效率较差的。 将数据保存在列表中比保存在数据框中访问效率高， 但数据框提供的功能更丰富。18.5 用性能分析工具找到程序效率的瓶颈
R的性能分析工具可以按较高频率（如每隔几毫秒）抽查正在被调用的函数， 因为函数可以嵌套调用， 所以会记录下来正在嵌套调用的各个函数。 因为抽查的速度很快， 所以会得到大量的被调用函数的抽样数据， 这样就可以进行概括， 得知哪些函数被调用最频繁， 调用的途径是怎样的。 基本R的
utils::Rprof()可以收集程序运行的性能分析用数据， 而utils::summaryRprof()则可以用文本格式提供性能分析的概括。在RStudio软件中， 借助于profvis扩展包， 我们可以从Profile菜单对选定行进行性能分析。 对比较复杂的程序， 应该将要分析的程序保存为一个R源文件， 并用
source()的方法将要运行的函数载入， 并用profvis()函数调用该函数并显示运行后性能分析结果。例18.6 以18.4的程序为例。 我们将第一个程序包装为一个函数
bad_copy()并存入文件bad_copy.R中， 然后用profvis进行性能分析。
bad_copy.R文件内容为：bad_copy <- function(){ M <- 1E5 x <- c() for(i in seq(M)){ x <- c(x, diff(range(runif(10)))) mean(x)安装好profvis扩展包后，运行如下的性能分析程序：
结果在RStudio中将显示如下的两个窗格内容：
p-prof-profiler-bc01.png
在Flame Graph窗格中， 下方显示了函数调用所占用的时间， 当鼠标光标移动到下方的某一个函数调用的色块时可以显示该函数的调用路径和耗时等信息。 可以看出是
c()函数和GC()函数占用了绝大多数时间， 而GC()是被c()调用的。 反复调用内存垃圾回收函数GC()一般说明程序有反复分配内存的问题， 从该图上半部分的Memery图表和Time图表也可以看出这一点， Memery图表中左侧为释放内存量， 右侧为分配内存量； Time图表看出程序运行集中在x被重复赋值的那一行。在Data窗格中， 点击Profvis后才展开详细的汇总统计图表。 可以看出耗时最多的是
c()函数调用， 但要注意的是， 这里的耗时也包括c()函数内部的diff(),range(),runif()函数的调用时间。18.6 R的计算函数
R中提供了大量的数学函数、统计函数和特殊函数， 可以打开R的HTML帮助页面， 进入“Search Enging & Keywords”链接， 查看其中与算术、数学、优化、线性代数等有关的专题。
这里简单列出一部分常用函数， 并对数学计算中常用函数
filter,fft,convolve进行较详细说明。18.6.1 数学函数
常用数学函数包括
abs,sign,log,log10,sqrt,exp,sin,cos,tan,asin,acos,atan,atan2,sinh,cosh,tanh。 还有gamma,lgamma(伽玛函数的自然对数)。用于取整的函数有
ceiling,floor,round,trunc,signif,as.integer等。 这些函数是向量化的一元函数。
choose(n,k)返回从\(n\)中取\(k\)的组合数。factorial(x)返回\(x!\)结果。combn(x,m)返回从集合\(x\)中每次取出不计次序的\(m\)个的所有不同取法， 结果为一个矩阵，矩阵每列为一种取法的\(m\)个元素值。18.6.2 概括函数
sum对向量求和,prod求乘积。
cumsum和cumprod计算累计， 得到和输入等长的向量结果。
diff计算前后两项的差分（后一项减去前一项）。
mean计算均值，var计算样本方差或协方差矩阵，sd计算样本标准差,median计算中位数，quantile计算样本分位数。cor计算相关系数。
colSums,colMeans,rowSums,rowMeans对矩阵的每列或每行计算总和或者平均值， 效率比用apply函数要高。
rle和inverse.rle用来计算数列中“连”长度及其逆向恢复， “连”经常用在统计学的随机性检验中。18.6.3 最值
max和min求最大和最小，cummax和cummin累进计算。
range返回最小值和最大值两个元素。对于
max,min,range， 如果有多个自变量可以把这些自变量连接起来后计算。
pmax(x1,x2,...)对若干个等长向量计算对应元素的最大值， 不等长时短的被重复使用。pmin类似。 比如，pmax(0, pmin(1,x))把x限制到\([0,1]\)内。18.6.4 排序
sort返回排序结果。 可以用decreasing=TRUE选项进行降序排序。sort可以有一个partial=选项， 这样只保证结果中partial=指定的下标位置是正确的。 比如:## [1] 2 1 3 4 5只保证结果的第三个元素正确。 可以用来计算样本分位数估计。
在
sort()中用选项na.last指定缺失值的处理， 取NA则删去缺失值， 取TRUE则把缺失值排在最后面， 取FALSE则把缺失值排在最前面。
order返回排序用的下标序列, 它可以有多个自变量， 按这些自变量的字典序排序。 可以用decreasing=TRUE选项进行降序排序。 如果指定选项method="radix", 则decreasing选项可以输入一个逻辑向量， 对每个参与排序的自变量分别指定是否降序排列。
sort.list函数与order功能相同但只支持按一个变量排序。
rank计算秩统计量，可以用ties.method指定同名次处理方法， 如ties.method=min取最小秩。
order,sort.list,rank也可以有na.last选项，只能为TRUE或FALSE。
unique()返回去掉重复元素的结果，duplicated()对每个元素用一个逻辑值表示是否与前面某个元素重复。## [1] 1 2 3## [1] FALSE FALSE TRUE FALSE TRUE
rev反转序列。18.6.5 一元定积分
integrate
integrate(f, lower, upper)对一元函数f计算从lower到upper的定积分。 使用自适应算法保证精度。函数的返回值不仅仅包含定积分数值， 还包含精度等信息。
允许使用
-Inf和Inf作为积分边界，如：18.6.6 一元函数求根
uniroot
uniroot(f, interval)对函数f在给定区间内求一个根，interval为区间的两个端点。 要求f在两个区间端点的值异号。 即使有多个根也只能给出一个。 如\(x(x-1)(x+1)\)函数在\([-2,2]\)的根：## $root ## [1] 0 ## $f.root ## [1] 0 ## $iter ## [1] 1 ## $init.it ## [1] NA ## $estim.prec ## [1] 2对于多项式， 可以用
polyroot函数求出所有的复根， 输入为升幂排列的多项式系数组成的向量。 如多项式\(x^3 - x\)和\(x^2 + 1\)的所有复根：18.6.7 离散傅立叶变换
fftR中
\begin{aligned} y_{k+1} =& \sum_{j=0}^{n-1} x_{j+1} \exp\left(-i 2\pi \frac{k j}{n} \right), \ k=0,1,\dots,n-1. \end{aligned} 注意没有除以\(n\)，结果是复数向量。fft函数使用快速傅立叶变换算法计算离散傅立叶变换。 设x为长度n的向量，y=fft(x)，则另外，若
y=fft(x),z=fft(y, inverse=TRUE), 则x == z/length(x)。快速傅立叶变换是数值计算中十分常用的工具， R软件包fftw可以进行优化的快速傅立叶变换。
18.6.8 用
filter函数作迭代R在遇到向量自身迭代时很难用向量化编程解决，
filter函数可以解决其中部分问题。filter函数可以进行卷积型或自回归型的迭代。 语法为下面用例子演示此函数的用途。
18.6.8.1 示例1：双侧滤波
对输入序列\(x_t, t=1,2,\dots,n\)， 希望进行如下滤波计算: \begin{aligned} y_{t} = \sum_{j=-k}^k a_j x_{t-j}, \ k+1 \leq t \leq n-k-1, \end{aligned} 其中\((a_{-k}, \dots, a_0, \dots, a_{k})\)是长度为\(2k+1\)的向量。 注意公式中\(a_j\)与\(x_{t-j}\)对应。
假设\(x\)保存在向量x中， \((a_{-k}, \dots, a_0, \dots, a_{k})\)保存在向量f中， \(y_{k+1}, \dots, y_{n-k}\)保存在向量y中, 无定义部分取NA, 程序可以写成
比如，设\(x = (1, 3, 7, 12, 17, 23)\), \((a_{-1}, a_0, a_{1}) = (0.1, 0.5, 0.4)\), 则 \begin{aligned} y_t = 0.1 \times x_{t+1} + 0.5 \times x_t + 0.4 \times x_{t-1}, \ t=2,3,\dots,5 \end{aligned} 用
filter()函数计算，程序为：y <- filter(c(1,3,7,12,17,23), c(0.1, 0.5, 0.4), method="convolution", sides=2) ## Time Series: ## Start = 1 ## End = 6 ## Frequency = 1 ## [1] NA 2.6 5.9 10.5 15.6 NA
filter的结果为一个时间序列对象， 可以用as.vector()函数转换为普通数值型向量。18.6.8.2 示例2: 单侧滤波
对输入序列\(x_t, t=1,2,\dots,n\)， 希望进行如下滤波计算: \begin{aligned} y_{t} = \sum_{j=0}^k a_j x_{t-j}, \ k+1 \leq t \leq n, \end{aligned} 其中\((a_0, \dots, a_{k})\)是长度为\(k+1\)的向量。 注意公式中\(a_j\)与\(x_{t-j}\)对应。
假设\(x\)保存在向量x中， \((a_0, \dots, a_{k})\)保存在向量f中， \(y_{k+1}, \dots, y_{n}\)保存在向量y中, 无定义部分取NA, 程序可以写成
比如，设\(x = (1, 3, 7, 12, 17, 23)\), \((a_{0}, a_1, a_{2}) = (0.1, 0.5, 0.4)\), 则 \begin{aligned} y_t = 0.1 \times x_{t} + 0.5 \times x_{t-1} + 0.4 \times x_{t-2}, \ t=3,4,\dots,6 \end{aligned}
y <- filter(c(1,3,7,12,17,23), c(0.1, 0.5, 0.4), method="convolution", sides=1) ## Time Series: ## Start = 1 ## End = 6 ## Frequency = 1 ## [1] NA NA 2.6 5.9 10.5 15.618.6.8.3 示例3: 自回归迭代
设输入\(e_t, t=1,2,\dots, n\)， \begin{aligned} y_t = \sum_{j=1}^k a_j y_{t-j} + e_t, \ t=1,2,\dots,n, \end{aligned} 其中\((a_1, \dots, a_k)\)是\(k\)个实数， \((y_{-k+1}, \dots, y_0)\)已知。
设\(x\)保存在向量
x中，\((a_1, \dots, a_k)\)保存在向量a中， \((y_1, \dots, y_n)\)保存在向量y中。如果\((y_{-k+1}, \dots, y_0)\)都等于零， 可以用如下程序计算\(y_1, y_2, \dots, y_n\):
如果\((y_0, \dots, y_{-k+1})\)保存在向量b中(注意与时间顺序相反)， 可以用如下程序计算\(y_1, y_2, \dots, y_n\):
设\(e = (0.1, -0.2, -0.1, 0.2, 0.3, -0.2)\), \((a_1, a_2) = (0.9, 0.1)\), \(y_{-1}=y_0=0\), 则 \begin{aligned} y_t = 0.9 \times y_{t-1} + 0.1 \times y_{t-2} + e_t, \ t=1,2,\dots,6 \end{aligned} 迭代程序和结果为y <- filter(c(0.1, -0.2, -0.1, 0.2, 0.3, -0.2), c(0.9, 0.1), method="recursive") print(y, digits=3) ## Time Series: ## Start = 1 ## End = 6 ## Frequency = 1 ## [1] 0.1000 -0.1100 -0.1890 0.0189 0.2981 0.0702这个例子中， 如果已知\(y_{0}=200, y_{-1}=100\), 迭代程序和结果为：
随机模拟要设法避免不同进程的随机数序列的重复可能， 以及同一进程中不同线程的随机数序列的重复可能。y <- filter(c(0.1, -0.2, -0.1, 0.2, 0.3, -0.2), c(0.9, 0.1), init=c(200, 100), method="recursive") print(y, digits=6) ## Time Series: ## Start = 1 ## End = 6 ## Frequency = 1 ## [1] 190.100 190.890 190.711 190.929 191.207 190.979parallel包提供了
parLapply()、parSapply()、parApply()函数， 作为lapply()、sapply()、apply()函数的并行版本， 与非并行版本相比， 需要用一个临时集群对象作为第一自变量。并行计算的工具在不断进步， 这部分内容的做法具有一定的参考价值， 但所用工具不一定是最新的， 欢迎读者批评指正。
18.7.1 例1：完全不互相依赖的并行运算
考虑如下计算问题: S_{k,n} = \sum_{i=1}^n \frac{1}{i^k}
下面的程序取
n为一百万，k为2到21，循环地用单线程计算。f10 <- function(k=2, n=1000){ s <- 0.0 for(i in seq(n)) s <- s + 1/i^k f11 <- function(n=1000000){ nk <- 20 v <- sapply(2:(nk+1), function(k) f10(k, n)) system.time(f11())[3] ## elapsed ## 1.95因为对不同的
k，f0(k)计算互相不依赖， 也不涉及到随机数序列， 所以可以简单地并行计算而没有任何风险。 先查看本计算机的虚拟核心（线程）数：用
makeCluster()建立临时的有8个节点的单机集群：用
parSapply()或者parLapply()关于\(k\)并行地循环：f12 <- function(n=1000000){ f10 <- function(k=2, n=1000){ s <- 0.0 for(i in seq(n)) s <- s + 1/i^k nk <- 20 v <- parSapply(cpucl, 2:(nk+1), function(k) f10(k, n)) system.time(f12())[3] ## elapsed ## 1.2并行版本速度提高了60%左右。
并行执行结束后， 需要解散临时的集群， 否则可能会有内存泄漏：
注意并行版本的程序还需要一些在每个计算节点上的初始化， 比如调入扩展包，定义函数， 初始化不同的随机数序列等。 parallel包的并行执行用的是不同的进程， 所以传送给每个节点的计算函数要包括所有的依赖内容。 比如，
f12()中内嵌了f10()的定义， 如果不将f10()定义在f12()内部， 就需要预先将f10()的定义传递给每个节点。parallel包的
clusterExport()函数可以用来把计算所依赖的对象预先传送到每个节点。 上面的f12()可以不包含f10()的定义， 而是用clusterExport()预先传递：cpucl <- makeCluster(nNodes) clusterExport(cpucl, c("f10")) f13 <- function(n=1000000){ nk <- 20 v <- parSapply(cpucl, 2:(nk+1), function(k) f10(k, n)) system.time(f13())[3] ## elapsed ## 1.16 stopCluster(cpucl)如果需要在每个节点预先执行一些语句， 可以用
clusterEvalQ()函数执行，如18.7.2 例2：使用相同随机数序列的并行计算
为了估计总体中某个比例\(p\)的置信区间， 调查了一组样本， 在\(n\)个受访者中选“是”的比例为\(\hat p\)。 令\(\lambda\)为标准正态分布的双侧\(\alpha\)分位数， 参数\(p\)的近似\(1-\alpha\)置信区间为 \frac{\hat p + \frac{\lambda^2}{2n}}{1 + \frac{\lambda^2}{n}} \frac{\lambda}{\sqrt{n}} \frac{\sqrt{\hat p (1 - \hat p) + \frac{\lambda^2}{4n}}} {1 + \frac{\lambda^2}{n}} 称为Wilson置信区间。
假设要估计不同\(1-\alpha\), \(n\), \(p\)情况下， 置信区间的覆盖率（即置信区间包含真实参数\(p\)的概率）。 可以将这些参数组合定义成一个列表， 列表中每一项是一种参数组合， 对每一组合分别进行随机模拟，估计覆盖率。 因为不同参数组合之间没有互相依赖的关系， 随机数序列完全可以使用同一个序列。
不并行计算的程序示例：
wilson <- function(n, x, conf){ hatp <- x/n lam <- qnorm((conf+1)/2) lam2 <- lam^2 / n p1 <- (hatp + lam2/2)/(1 + lam2) delta <- lam / sqrt(n) * sqrt(hatp*(1-hatp) + lam2/4) / (1 + lam2) c(p1-delta, p1+delta) f20 <- function(cpar){ set.seed(101) conf <- cpar[1] n <- cpar[2] p0 <- cpar[3] nsim <- 100000 cover <- 0 for(i in seq(nsim)){ x <- rbinom(1, n, p0) cf <- wilson(n, x, conf) if(p0 >= cf[1] && p0 <= cf[2]) cover <- cover+1 cover/nsim f21 <- function(){ lp <- expand.grid(c(0.8, 0.9), c(30, 100), c(0.5, 0.1)) |> as.matrix() |> t() |> as.data.frame() |> as.list() res <- sapply(lp, f20) system.time(f21())[3] ## elapsed ## 6.52约运行6.52秒。
改为并行版本：
library(parallel) nNodes <- 8 cpucl <- makeCluster(nNodes) clusterExport(cpucl, c("f20", "wilson")) f22 <- function(){ lp <- expand.grid(c(0.8, 0.9), c(30, 100), c(0.5, 0.1)) |> as.matrix() |> t() |> as.data.frame() |> as.list() res <- parSapply(cpucl, lp, f20) system.time(f22())[3] ## elapsed ## 2.83 stopCluster(cpucl)运行约2.83秒， 速度提高一倍以上。 这里模拟了8种参数组合， 每种参数组合模拟了十万次， 每种参数组合模拟所用的随机数序列是相同的。
18.7.3 例3：使用独立随机数序列的并行计算
大量的耗时的统计计算是随机模拟， 有时需要并行计算的部分必须使用独立的随机数序列。 比如，需要进行一千万次重复模拟，每次使用不同的随机数序列， 可以将其分解为10组模拟，每组模拟一百万次， 这就要求这10组模拟使用的随机数序列不重复。
R中实现了
L'Ecuyer的多步迭代复合随机数发生器， 此随机数发生器周期很长， 而且很容易将发生器的状态前进指定的步数。 parallel包的nextRNGStream()函数可以将该发生器前进到下一段的开始， 每一段都足够长， 可以用于一个节点。以Wilson置信区间的模拟为例。 设\(n=30\), \(p=0.01\), \(1-\alpha=0.95\)， 取重复模拟次数为1千万次，估计Wilson置信区间的覆盖率。 单线程版本为：
f31 <- function(nsim=1E7){ set.seed(101) n <- 30; p0 <- 0.01; conf <- 0.95 cover <- 0 for(i in seq(nsim)){ x <- rbinom(1, n, p0) cf <- wilson(n, x, conf) if(p0 >= cf[1] && p0 <= cf[2]) cover <- cover+1 cover/nsim system.time(cvg1 <- f31())[3] ## elapsed ## 78.11单线程版本运行了大约78秒。
改成并行版本。 比例2多出的部分是为每个节点分别计算一个随机数种子将不同的种子传给不同节点。 parallel包的
clusterApply()函数为临时集群的每个节点分别执行同一函数， 但对每个节点分别使用列表的不同元素作为函数的自变量。library(parallel) nNodes <- 8 cpucl <- makeCluster(nNodes) each.seed <- function(s){ assign(".Random.seed", s, envir = .GlobalEnv) RNGkind("L'Ecuyer-CMRG") set.seed(101) seed0 <- .Random.seed seeds <- as.list(1:nNodes) for(i in 1:nNodes){ # 给每个节点制作不同的种子 seed0 <- nextRNGStream(seed0) seeds[[i]] <- seed0 ## 给每个节点传送不同种子： junk <- clusterApply(cpucl, seeds, each.seed) f32 <- function(isim, nsimsub=10000){ n <- 30; p0 <- 0.01; conf <- 0.95 cover <- 0 for(i in seq(nsimsub)){ x <- rbinom (1, n, p0) cf <- wilson(n, x, conf) if(p0 >= cf[1] && p0 <= cf[2]) cover <- cover+1 cover clusterExport(cpucl, c("f32", "wilson")) f33 <- function(nsim=1E7){ nbatch <- 40 nsimsub <- nsim / nbatch cvs <- parSapply(cpucl, 1:nbatch, f32, nsimsub=nsimsub) sum(cvs)/(nsim*nbatch) system.time(cvg2 <- f33())[3] ## elapsed ## 30.72 stopCluster(cpucl)并行版本运行了大约41秒，速度提高一倍以上。 从两个版本各自一千万次重复模拟结果来看， 用随机模拟方法得到的覆盖率估计的精度大约为3位有效数字。
更大规模的随机模拟问题， 可以考虑使用多CPU的计算工作站或者服务器， 或用多台计算机通过高速局域网组成并行计算集群。
还有一种选择是租用云计算服务。