模拟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("http://www.google.cn")) #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){     x      r[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. Information          education2. HS Grad                       4.5107                      11.2009     education3. Some College     education4. 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        75      0.00        0"model.matrix.default"       0.04        50      0.02       25"model.matrix"               0.04        50      0.00        0""                0.02        25      0.02       25"conditionCall"              0.02        25      0.02       25"stop"                       0.02        25      0.02       25".External2"                 0.02        25      0.00        0"[[.data.frame"              0.02        25      0.00        0"[["                         0.02        25      0.00        0"base::try"                  0.02        25      0.00        0"contrasts                0.02        25      0.00        0"eval"                       0.02        25      0.00        0"model.frame.default"        0.02        25      0.00        0……$sample.interval[1] 0.02$sampling.time[1] 0.08unlink(tmp)

sysyem.time( )将从R分析器中得到的数据规范化的方式：

・by total：把每个函数依照总运行时间排序，第一个通常是顶层函数

因为lm位于顶层，耗费75%的时间，但它只是调用真正进行运算的辅助函数，因此知道顶层函数使用的时间用处不大

・by self：先减去调用栈中低级函数花的时间再进行排序，即每个函数真正的耗费时间

编辑：李雪纯 冯文清 校审：张健 罗鹏