模拟simulation
➢概率函数
概率函数通常用来生成特征已知的模拟数据，以及在统计函数中计算概率值。
对于任意分布有四种基本函数：
正态分布：rnorm( )、dnorm( )、pnorm( )、qnorm( )； Gamma分布：rgamma( )、dgamma( )、pgamma( )和 qgamma( )； 泊松分布：rpois( )、dpois( )、ppois( ) 和 qpois( )。
dnorm(x, mean = 0, sd = 1, log = FALSE) pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) rnorm(n, mean = 0, sd = 1)
正态分布概率函数默认为标准正态分布。
・dnorm( ) 的log选项能求密度的对数值 (log of density) ・pnorm( ) 和 qnorm( ) 中log.p选项也可以对概率求对数 ・lower.tail选项默认计算分布的左尾 (lower tail) ・若Φ是标准正态分布的累积分布函数，则pnorm(q)=Φ(q)且qnorm(p)=Φ-1(p) (qnorm() 等于Φ的反函数)
set.seed( ) 设置随机数生成种子 (seed)，可以重复生成相同的随机数。
set.seed(1)rnorm(5) #生成5个服从标准正态分布的随机数[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078rnorm(5)[1] -0.8204684 0.4874291 0.7383247 0.5757814 -0.3053884set.seed(1)rnorm(5)[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078
生成服从泊松分布的随机数
x10, x[1] 1 1 1 0 2 1 2 0 1 1summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 0 1 1 1 1 2#泊松分布中，均值应该等于平均发生率
估计泊松分布的累积分布函数
ppois(2, 2)[1] 0.6766764#在平均发生率为2的泊松分布中出现小于等于2的随机变量的概率，Pr(x <= 2)ppois(4, 2)[1] 0.947347# Pr(x <= 4)ppois(6, 2)[1] 0.9954662# Pr(x <= 6)
➢从线性模型模拟随机值
假设线性模型y=β0+β1x+ε，其中β0等于0.5，β1等于2，随机噪声ε服从标准差为2的正态分布。
当x服从正态分布时：
set.seed(20)x <- rnorm(100) #生成服从正态分布的随机变量xe <- rnorm(100, 0, 2) #生成随机噪声y<-0.5+2*x+e summary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. -6.4084 -1.5402 0.6789 0.6893 2.9303 6.5052plot(x, y)
用回归模型模拟x和y，两者之间存在明显线性关系，符合回归模型。
当x服从二项分布时：
set.seed(10)x <- rbinom(100, 1, 0.5) e <- rnorm(100, 0, 2)y<-0.5+2*x+esummary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.4936 -0.1409 1.5767 1.4322 2.8397 6.9410 plot(x, y)
x是二元变量时，y依旧是连续的且服从正态分布，x和y之间存在线性趋势。
➢从广义线性模型模拟随机值
假设y服从泊松分布，Y ~ Poisson(μ)；log(μ)服从一个截距为β0，斜率为β1的线性函数，logμ=β0 +β1x；β0为0.5，β1为0.3。
set.seed(1)x <- rnorm(100) #生成服从正态分布的随机变量xlog.mu<-0.5+0.3*xy <- rpois(100, exp(log.mu)) #模拟生成100个服从泊松分布的随机变量summary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 1.00 1.00 1.55 2.00 6.00plot(x, y)
可以看到y不是连续变量，而是计数变量，但x和y之间存在线性关系，当x越大，y也增加。
➢随机抽样
sample函数可以从给定的一组对象中随机抽取样本。
set.seed(1)sample(1:10, 4)[1] 9 4 7 1 sample(1:10, 4)[1] 2 7 3 6 sample(letters, 5)[1] "r" "s" "a" "u" "w"
默认为无放回抽样，不会得到重复的元素。
sample(1:10)[1] 10 6 9 2 1 5 8 4 3 7
如果不指定其他参数，返回这些整数重新排列后的结果，向量里1到10的顺序被随机重新排列。
sample(1:10, replace = TRUE)[1] 5 5 2 10 9 1 4 3 6 10
指定replace=TRUE，函数进行10次有放回抽样。
分析器profiler
在开发规模较大的程序或者做大型数据分析时，运行代码可能耗费很长时间，分析器可以分析代码的运行时间并提供一些建议解决问题，优化代码。
➢system.time( )
system.time( )取出并分析任意一条R语句，返回运行语句所需时间。如果执行分析时代码有错误，返回代码运行直到错误发生所用的时间。
●用户时间 (user time)：运行这行语句需要占用中央处理器的时间。
●运行时间 (elapsed time)：所经历的时间，也称时钟时间 (wall clock time)。
・通常情况下，用户时间和运行时间是比较接近的，计算机执行函数或者语句所需时间大致上等于等待的时间。
・运行时间比用户时间长，说明中央处理器实际可能耗费很多时间在其它后台事项上。
・运行时间比用户时间短，通常是因为机器有多核或多处理器且都处于可利用状态。
#运行时间比用户时间长system.time(readLines("")) #readLines() 从远程服务器读取网页 user system elapsed 0.002 0.002 0.305
中央处理器实际上没有耗费很多时间运行代码，大部分时间用在等待网络，等待数据经过网络传递到计算机里，等待网络对数据的来回处理上。
#运行时间比用户时间短hilbert function(n) { i<-1:n 1/outer(i-1,i,"+")}x 1000) system.time(svd(x)) #svd()计算这个矩阵的奇异值分解 用户 系统 流逝 1.605 0.094 0.742
svd函数使用苹果机的 Accelerate 框架，一个多线程线性代数库，把计算任务分给两个核进行。
#运行时间与用户时间大致相等system.time({n r for(i in 1:n){xr[i] } }) 用户 系统 流逝 0.099 0.006 0.109
不需要多线程或网络活动的简单表达式，用户时间和运行时间基本相同。
➢Rprof( )
Rprof( )启动R中的分析器，可以在规律的样本区间内记录追踪函数调用栈，即每过0.02秒打印输出函数的调用栈。(当函数运行时间少于0.02秒时无效)
system.time函数适用于较小型、较简单、很明确问题发生位置的程序。如果不明确问题在哪里，可以用Rprof函数。(注意，不可以同时使用Rprof( )和sysyem.time( )。)
summaryRprof( )总结分析器的原始输出结果，计算每个函数各自耗费的时间。
library(ISLR)Wage Rprof(tmp lm(wage~age + jobclass + education,data=Wage) #运行lm函数Call:lm(formula = wage ~ age + jobclass + education, data = Wage)Coefficients:(Intercept)age 59.5800 0.5554 jobclass2. Informationeducation2. HS Grad 4.5107 11.education3. Some Collegeeducation4. College Grad 23.3304 38.3862 education5. Advanced Degree 62.9115 Rprof()summaryRprof(tmp)$by.self #真正耗时的函数 self.time self.pct total.time total.pct"model.matrix.default"0.02 25 0.04 50""0.02 25 0.02 25"conditionCall" 0.02 25 0.02 25"stop" 0.02 25 0.02 25$by.total #顶层函数为lm，耗费0.06s total.time total.pct self.time self.pct"lm"0.06 750.00 0"model.matrix.default" 0.04 500.02 25"model.matrix"0.04 500.00 0""0.02 250.02 25"conditionCall" 0.02 250.02 25"stop" 0.02 250.02 25".External2" 0.02 250.00 0"[[.data.frame" 0.02 250.00 0"[["0.02 250.00 0"base::try" 0.02 250.00 0"contrasts0.02 250.00 0"eval" 0.02 250.00 0"model.frame.default" 0.02 250.00 0……$sample.interval[1] 0.02$sampling.time[1] 0.08unlink(tmp)
sysyem.time( )将从R分析器中得到的数据规范化的方式：
・by total：把每个函数依照总运行时间排序，第一个通常是顶层函数
因为lm位于顶层，耗费75%的时间，但它只是调用真正进行运算的辅助函数，因此知道顶层函数使用的时间用处不大
・by self：先减去调用栈中低级函数花的时间再进行排序，即每个函数真正的耗费时间
编辑：李雪纯 冯文清 校审：张健 罗鹏
