R语言教程-R编程5

18 R程序效率

并行计算

现代电脑一般有多个核心，通常R运行并不能利用全部的CPU能力， 仅能利用其中的一个虚拟核心。 使用特制的BLAS库（非R原有）可以并发运行多个线程， 一些R扩展包也可以利用多个线程。 利用多台计算机、多个CPU、CPU中的多核心和多线程同时完成一个计算任务称为并行计算。

R的parallel包提供了一种比较简单的利用CPU多核心的功能，如果有多个任务互相之间没有互相依赖， 就可以分解到多个计算机、多个CPU、多个虚拟核心中并行计算。统计计算中最常见耗时计算任务是随机模拟， 随机模拟要设法避免不同进程的随机数序列的重复可能， 以及同一进程中不同线程的随机数序列的重复可能。

parallel包提供了parLapply()、parSapply()、parApply()函数， 作为lapply()、sapply()、apply()函数的并行版本， 与非并行版本相比， 需要用一个临时集群对象作为第一自变量。

例1：完全不互相依赖的并行运算

$ S_{k,n} = \displaystyle\sum_{i=1}^n \frac{1}{i^k}$

f10 <- function(k=2:21, n=1000000){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
}
system.time(f10()) # 时间流逝1.34秒
f10()
f10 <- function(k=2, n=1000){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
}

f11 <- function(){
  v <- sapply(2:21, f10, n=1000000)
  v
}
system.time(f11()) # 时间流逝2.09秒
f11()
f10 <- function(k=2, n=1000){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
}
# 将n =1000000放进function中预设的速度比在sapply函数中再赋值快。
f11 <- function(n=1000000){
  v <- sapply(2:21, f10)
  v
}
system.time(f11()) # 时间流逝0.01秒
f11()

检测上面的程序计算结果是否正确：

f10 <- function(k=2, n=1000){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
}
f11 <- function(n=1000000){
  nk <- 20
  v <- sapply(2:(nk+1), function(k) f10(k, n))
  v
}
system.time(f11())[3] # 施加流逝2.09秒
f11()

对于不同的k，S~k,n~计算互相独立，也不涉及随机数序列，所以可以简单地并行计算而没有任何风险。

# 使用parallel包查看计算机的虚拟核心（线程）数
library(parallel)
detectCores()
# 用makeCluster()建立临时的有8个节点的单机集群
nNodes <- 8
cpuc <- makeCluster(nNodes)
library(parallel)

f13 <- function(n=1000000){
  f12 <- function(k=2,n=100){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
  }
  
  v <- parSapply(cpuc, 2:21, f12)
  v
}
system.time(f13()) # 时间流逝0.00秒
f13()

# 并行执行结束后，需要解散临时的集群， 否则可能会有内存泄漏
stopCluster(cpuc)

parallel包的并行执行用的是不同的进程， 所以传送给每个节点的计算函数要包括所有的依赖内容。 比如，f2()中内嵌了f0()的定义， 如果不将f0()定义在f2()内部， 就需要预先将f0()的定义传递给每个节点。

parallel包的并行执行用的是不同的进程， 所以传送给每个节点的计算函数要包括所有的依赖内容。 比如，f2()中内嵌了f0()的定义， 如果不将f0()定义在f2()内部， 就需要预先将f0()的定义传递给每个节点。

library(parallel)
cup <- makeCluster(8)
f10 <- function(k=2, n=100){
  s <- 0.0
  for(i in seq(n)) s <- s + 1/i^k
  s
}
clusterExport(cup,c("f10"))
f14 <- function(n=1000000){
  v <- parSapply(cup, 2:21, f10)
  v
}
system.time(f14()) # 时间流逝0.04秒
f14()
# 关闭集群
stopCluster(cup)

19 随机模拟

随机数

所谓随机数，实际是“伪随机数”， 是从一组起始值（称为种子）， 按照某种递推算法向前递推得到的。 所以，从同一种子出发，得到的随机数序列是相同的。

为了得到可重现的结果， 随机模拟应该从固定不变的种子开始模拟。 用set.seed(k)指定一个编号为k的种子， 这样每次从编号k种子运行相同的模拟程序可以得到相同的结果。
还可以用set.seed()加选项kind=指定后续程序要使用的随机数发生器名称， 用normal.kind=指定要使用的正态分布随机数发生器名称。
R提供了多种分布的随机数函数，如runif(n)产生n个标准均匀分布随机数， rnorm(n)产生n个标准正态分布随机数。

runif(5)
rnorm(5)
# 查看R中提供的不同概率分布。
?Distributions

sample()函数

sample()函数从一个有限集合中无放回或有放回地随机抽取， 产生随机结果。

# 模拟10次抛硬币的值
sample(c('正面', '反面'), size = 10,
       prob=c(0.5,0.5),replace = T)

sample()的选项size指定抽样个数， prob指定每个值的概率， replace=TRUE说明是有放回抽样。

如果要做无放回等概率的随机抽样， 可以不指定prob和replace(缺省是FALSE)。 比如，

# 1:10随机抽取5个
sample(1:10, size=5)
# 1:10中等概率无放回随机抽样直到每一个都被抽过
sample(10)

dplyr包的sample_n()函数与sample()类似， 但输入为数据框， 输出为随机抽取的数据框行子集。

随机模拟示例

线性回归模拟

y=10+2x+ε，ε∼N(0,0.52)y=10+2x+\varepsilon，ε\thicksim N(0,0.5^2)

LaTex语法参考