R语言教程-R编程6

随机模拟示例

线性回归模拟

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

LaTex语法参考

标准差(Standard Deviation) 和 标准误差(Standard Error)

假设有样本量的一组样本， R函数lm()可以 可以得到截距, 斜率的估计, 以及相应的标准误差, 。 样本可以模拟产生。

x <- sample(1:10, size = 10, replace = T)
eps <- rnorm(10, 0, 0.5)
y <- 10 + 2*x + eps
# 计算线性回归
lm(y ~ x)
# 计算线性回归的多种统计量
summary(lm(y ~ x))
# 返回一个矩阵， 包括的估计值、标准误差、t检验统计量、检验p值
summary(lm(y ~ x))$coefficients
# 把上述矩阵的前两列拉直成一个向量返回
c(summary(lm(y ~ x))$coefficients[,1:2])

用replicate(重复次数，重复语句)执行多次模拟， 返回向量或矩阵结果， 返回矩阵时，每列是一次模拟的结果。

reg.sim <- function(){
  set.seed(1)
  resm <- replicate(1000, {
    x <- sample(1:10, size = 10, replace = T)
    eps <- rnorm(10, 0, 0.5)
    y <- 10 + 2*x + eps
    c(summary(lm(y ~ x))$coefficients[,1:2])
  })
  # 将行的数值加上行名
  row.names(resm) <- c('a', 'b', 'SE.a', 'SE.b')
  cat(1000,'次模拟的平均值：\n')
  print(apply(resm, 1, mean))
  cat(1000,'次模拟的标准差：\n')
  print(apply(resm, 1, sd))
}
set.seed(1)
reg.sim()

apply()函数参数简介：简单的说，apply函数经常用来计算矩阵中行或列的均值、和值的函数
apply(b,1,sum)
上面的指令代表对矩阵b进行行计算，分别对每一行进行求和。函数涉及了三个参数：

第一个参数是指要参与计算的矩阵；

第二个参数是指按行计算还是按列计算，1——表示按行计算，2——按列计算；

第三个参数是指具体的运算参数。

—

另一种写法：

reg.sim <- function(
  a=10, b=2, sigma=0.5, 
  n=10, B=1000){
  set.seed(1)
  resm <- replicate(B, {
      x <- sample(1:10, size=n, replace=TRUE)
      eps <- rnorm(n, 0, 0.5)
      y <- a + b * x + eps
      c(summary(lm(y ~ x))$coefficients[,1:2])
  })
  # 转置为列的数值
  resm <- t(resm)
  # 将列的数值加上列名
  colnames(resm) <- c('a', 'b', 'SE.a', 'SE.b')
  cat(B, '次模拟的平均值:\n')
  print( apply(resm, 2, mean) )
  cat(B, '次模拟的标准差:\n')
  print( apply(resm, 2, sd) )
}
set.seed(1)
reg.sim()

核密度的bootstrap置信区间

R自带的数据框faithful内保存了美国黄石国家公园Faithful火山的272次爆发持续时间和间歇时间。 为估计爆发持续时间的密度，可以用核密度估计方法， R函数density可以执行此估计， 返回NN个格子点上的密度曲线坐标:

x <- faithful$eruptions
est0 <- density(x)
plot(est0)

核密度估计是统计估计， 为了得到其置信区间（给定每个x坐标，真实密度f(x)的单点的置信区间）， 采用如下非参数bootstrap方法：

# 利用replicate()函数实现
set.seed(2)
resm <- replicate(10000,{
  x1 <- sample(x, replace = T)
  density(x1, form=min(est0$x), to=max(est0$x))$y
})

CI <- apply(resm, 1, quantile, c(0.025, 0.975))

plot(est0, ylim=range(CI), type = 'n')
polygon(c(est0$x, rev(est0$x)),
        c(CI[1,], rev(CI[2,])),
          col="red", border=F)

lines(est0, lwd=2)

quantile()函数计算B个bootstrap样本的2.5%和97.5%分位数， 循环用apply()函数表示。

polygon()函数指定一个多边形的顺序的顶点坐标用col=指定的颜色填充， 本程序中实现了置信下限与置信上限两条曲线之间的颜色填充。 lines()函数绘制了与原始样本对应的核密度估计曲线。

函数进阶

函数调用的各种形式

前面学习的部分就有涉及，这里只大概的看一下，不重复运行

嵌套定义与句法作用域(lexical scoping)

R语言允许在函数体内定义函数。 比如，

x <- -1
f0 <- function(x){
  f1 <- function(){
    x + 100
  }
  f1()

其中内嵌的函数f1()称为一个closure(闭包)。

句法作用域(lexical scoping)， 即函数运行中需要使用某个变量时， 从其定义时的环境向外层逐层查找， 而不是在调用时的环境中查找。

“句法作用域”指的是函数调用时查找变量是查找其定义时的变量对应的存储空间， 而不是定义时变量所取的历史值。函数运行时在找到某个变量对应的存储空间后， 会使用该变量的当前值，而不是函数定义的时候该变量的历史值。 这种规则称为动态查找(dynamic lookup)。

# 此函数中两个x的值层级相同，-1为历史值，而10000为当前值，所以x取10000
f0 <- function(){
  x <- -1 # 此x值为历史值
  f1 <- function(){
    x + 100
  }
  x <- 1000 # 此x值为当前值
  f1()
}
f0()

对比：

# x的当前值为-1
f0 <- function(){
  x <- 10000 # 此x值为历史值
  f1 <- function(){
    x + 100
  }
  x <- -1 # 此x值为当前值
  f1()
}
f0()

句法作用域：使用变量时，从其定义的位置向外层逐层查找
动态查找：在相同层级查找某个变量时，使用其当前值，而不是历史值

当前值和历史值课可理解为：书写的顺序一般从上到下，越靠下的值越接近当前，则使用这个最靠下的值为当前值。

句法作用域与动态查找一个说的是如何查找某个变量对应的存储空间， 一个说的是使用该存储空间何时的存储值， 程序运行时两个规则需要联合使用。

辅助嵌套函数

输入一元二次方程ax2+bx+c=0ax^2+bx+c=0的三个系数，求解一元二次方程：

# 一元二次方程求解
solve.sqe <- function(x){
  fd <- function(a, b, c) b^2 - 4*a*c
  d <- fd(x[1], x[2], x[3])
  if(d >= 0){
    return( (-x[2] + c(1,-1)*sqrt(d))/(2*x[1]) )
  } 
 # 返回无解
   else return('无实数解')
# 虚数求解
#  else {
#    return( complex(real=-x[2], #imag=c(1,-1)*sqrt(-d))/(2*x[1]) )
#  }
}

solve.sqe(c(1, -2, 1))
solve.sqe(c(1, -2, 2))

在这个函数中内嵌的函数fd仅起到一个计算二次判别式的作用， 没有用到任何的闭包特性， 其中的形参变量a, b, c都是局部变量。这样的函数只能在定义它的函数内运行，不能被直接调用，可以看成是函数内的私有函数。

简单写为以下的表达式更容易理解：

# 一元二次方程求解
solve.sqe <- function(x){
  d <- x[2]^2 - 4*x[1]*x[3]
  if(d >= 0){
    return( (-x[2] + c(1,-1)*sqrt(d))/(2*x[1]) )
  } 
 # 返回无解
   else return('无实数解')
}
solve.sqe(c(1, -2, 1))
solve.sqe(c(1, -2, 2))

懒惰求值

f <- function(x, y=ifelse(x>0, TRUE, FALSE)){
  x <- -111
  if(y) x*2 else x*10
}
f(5)

最后值为-1110，虽然形参x输入的实参值为5, 但是这时形参y并没按x=5被赋值为TRUE, 而是到函数体中第二个语句需要用到y时才被求值， 这时x的值已经变成了-111， 故y的值是FALSE。

程序调试

f1 <- function(x) f2(x)
f2 <- function(x) 1/x
f1("abc")
## Error in 1/x : 二进列运算符中有非数值参数
traceback() # run selected line
## 2: f2(x) at #1
## 1: f1("abc")

结果是所谓的反向追踪(traceback)， 一般编程语言中称为调用栈(calling stack)。 这个输出是从下向上堆叠显示， 下层是调用方， 上层是被调用方。

跟踪调试

# 考虑如下错误的函数
f <- function(x){
  for(i in 1:n){
    s <- s + x[i]
  }
}
print(f(1:5))
## Error in 1:n : NA/NaN参数

R提供了一个browser()函数， 在程序中插入对browser()函数的调用， 可以进入跟踪调试状态， 可以实时地查看甚至修改运行时变量的值。

f <- function(x){
  browser()
  for(i in 1:n){
    s <- s + x[i]
  }
}

当在RStudio中调试运行时， 程序编辑窗口将显示当前要运行的程序行， 用命令行窗口(Console)的Next快捷图标可以运行到下一行。
可以点击Console的continue按钮，可以直接一次运行完代码。
最后改为：

f <- function(x){
  browser()
  n <- length(x)
  s <- 0
  for(i in 1:n){
    s <- s + x[i]
  }
  s
}

自定义函数应该用各种不同输入测试其正确性和稳定性。 比如，上面的函数当自变量x为零长度向量时应该返回0才合适， 但是上面的写法会返回一个numeric(0)结果， 这个结果表示长度为零的向量：

f <- function(x){
  n <- length(x)
  s <- 0
  for(i in 1:n){
    s <- s + x[i]
  }
  s
}
f(0)
f(numeric(0))

程序输入了零长度自变量， 我们期望其输出为零而不是numeric(0)。 在自变量x为零长度时， 函数中for(i in 1:length(x)应该一次都不进入循环。

把这里的1:length(x)改成seq_along(x)解决了问题， seq_along(x)生成x的下标序列， 如果x是零长度的则下标序列为零长度向量。

函数不需要修改后， 可以把对browser()的调用删除或注释掉:

f <- function(x){
  # browser()
  s <- 0
  for(i in seq_along(x)){
    s <- s + x[i]
  }
  s
}

这里的求向量元素和知识一个简单的示例，可以用R内置的函数sum完成。

实际上，许多我们认为需要自己编写程序做的事情， 在R网站都能找到别人已经完成的扩展包。

条件断点

用browser()函数与R的if结构配合可以制作条件断点。 比如， 在调试带有循环的程序时， 发现错误发生在循环内， 如果从循环开始就跟踪运行， 会浪费许多时间。 设已知错误发生在循环变量i等于501的时候， 就可以在循环内插入：

if(i == 501) browser()

这样就可以在更接近出错的位置进入跟踪运行。

还可以用debug(f)命令对函数f开启跟踪运行， 这时每次调用f()时都会自动进入跟踪运行状态。 用undebug(f)取消对f的这种操作。

出错调试选项
比较长的程序在调试时如果从开头就跟踪， 比较耗时。可以设置成出错后自动进入跟踪模式， 检查出错时的变量值。只要进行如下设置：

options(error=recover)

则在出错后可以选择进入出错的某一层函数内部, 在browser环境中跟踪运行。
在RStudio中调试某一源程序文件时， 可以选择“Debug–On Error”菜单， 并选择“Break in Code”， 就可以在出错时自动在出错处进入跟踪状态。

f <- function(x){
  options(error=recover)
  for(i in 1:n){
    s1 <- s1 + x[i]
  }
  s1
}
f(0)
f(numeric(0))

在console中的Selection后面输入了1，就进入了函数内部跟踪。 用Q终止运行并退出整个browser跟踪。 当函数调用函数时可以选择进入哪一个函数进行跟踪。 如果在RStudio中设置了“Break in Code”， 会自动在出错处进入跟踪运行状态。ESC键退出当前跟踪，返回命令行状态。

提前发现错误状态：

  1. stop(s)使程序运行出错停止， 其中s是一个字符型对象， 用来作为显示的出错信息。
  2. 用warning(s)提交一个警告信息， 其中s是字符型的警告信息。 警告信息的显示可能比实际运行要延迟一些。
  3. 有些警告信息实际是错误， 用options()的warn参数可以设置警告级别， 如设置warn=2则所有警告当作错误处理。warn=0是默认做法， warn=1表示不延迟显示。
  4. message()是提示性的信息输出。 message()不会像warning()那样延迟显示。 比如， 长时间的等待之前给出一个提示， 正在读写文件或者网络时给出提示， 等等。 与cat()等相比较， cat()是用户要求的输出， 而message()是程序员对用户的提示。

stop()产生错误状态， 停止运行； 用warning()产生警告， 用message()产生提示

预防性设计

函数stopifnot可以指定自变量的若干个条件， 当自变量不符合条件时自动出错停止。
函数f()需要输入两个数值型向量x, y, 需要长度相等：

f <- function(x, y){
  stopifnot(is.numeric(x),
            is.numeric(y),
            length(x)==length(y))
  ## 函数体程序语句...
}

上面的程序可以理解为python中的指定x,y为数值型输入。另外还要求x,y长度相等。