北京大学R语言教程(李东风)第40章: 随机模拟

40.1 随机数

随机模拟是统计研究的重要方法, 另外许多现代统计计算方法(如MCMC)也是基于随机模拟的。 R中提供了多种不同概率分布的随机数函数, 可以批量地产生随机数。 一些R扩展包利用了随机模拟方法,如boot包进行bootstrap估计。

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

为了得到可重现的结果, 随机模拟应该从固定不变的种子开始模拟。 用set.seed(k)指定一个编号为k的种子, 这样每次从编号k种子运行相同的模拟程序可以得到相同的结果。

还可以用set.seed()加选项kind=指定后续程序要使用的随机数发生器名称, 用normal.kind=指定要使用的正态分布随机数发生器名称。

R提供了多种分布的随机数函数,如runif(n)产生n个标准均匀分布随机数, rnorm(n)产生n个标准正态分布随机数。 例如:

round(runif(5), 2)
## [1] 0.44 0.56 0.93 0.23 0.22
round(rnorm(5), 2)
## [1] -0.20  1.10 -0.02  0.16  2.02

注意因为没有指定种子,每次运行会得到不同的结果。

在R命令行运行

?Distributions

可以查看R中提供的不同概率分布。

40.2 sample()函数

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

例如,为了设随机变量X取值于{正面,反面}, 且P(X=正面)=0.7=1−P(X=反面), 如下程序产生X的10个随机抽样值:

sample(c('正面', '反面'), size=10, 
       prob=c(0.7, 0.3), replace=TRUE)
## [1] "反面" "反面" "反面" "反面" "正面"
## [6] "正面" "正面" "正面" "反面" "反面"

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

如果要做无放回等概率的随机抽样, 可以不指定probreplace(缺省是FALSE)。 比如,下面的程序从1:10随机抽取4个:

sample(1:10, size=4)
## [1]  1  5  8 10

如果要从1:n中等概率无放回随机抽样直到每一个都被抽过,只要用如:

sample(10)
## [1]  3  5  9  2 10  7  4  1  6  8

这实际上返回了1:10的一个重排。

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

40.3 随机模拟示例

40.3.1 估计期望值

设随机变量或随机向量X有复杂的分布, 使得期望θ=Eh(X)很难计算。 如果可以得到X的N个独立同分布抽样值X1,X2,…,XN, 则可估计θ为:

θ̂ =1N∑i=1Nh(Xi),

这时Eθ̂ =θ, 估计无偏; 均方误差为

MSE==E|θ̂ −θ|2=Var(θ̂ )=SE2Var(h(X))N≈S2NN,其中

S2N=1N−1∑i=1N(h(Xi)−θ̂ )2是样本方差。 由强大数律可知N→∞时θ̂ 以概率1收敛到θ, 由中心极限定理可知N充分大时θ̂ 有近似正态分布

N(θ,SE2).

例如,考虑正方形区域Ω={(x,y):x∈[0,1],y∈[0,1]}, 以及其中的四分之一扇形A={(x,y):(x,y)∈Ω,x2+y2≤1}。 设X服从Ω上的均匀分布, 令

Y={1,0,当X∈A,其它

则Y服从两点分布, 概率为

p=EY=14π1212=π4.

在Ω中投入N=10000个点, 得到N个Yi,i=1,2,…,N的值,估计π为

π̂ =4Y¯=4N∑i=1NYi.

Eπ̂ =π, 估计的根均方误差可估计为

RMSE=E|π̂ −π|2‾‾‾‾‾‾‾‾‾√=4Var(Y)‾‾‾‾‾‾‾√N‾‾√≈4SNN‾‾√.

程序如下:

est_pi <- function(N){
  set.seed(101)
  x1 <- runif(N, 0, 1)
  x2 <- runif(N, 0, 1)
  y <- as.numeric(x1^2 + x2^2 <= 1)
  hat_pi <- 4*mean(y)
  se <- 4 * sd(y) / sqrt(N)
  cat("N = ", N, " pi估计值 =", hat_pi, " SE =", se, "\n")
  
  invisible(list(N=N, hat_pi = hat_pi, SE = se))
}
est_pi(1E4)
## N =  10000  pi估计值 = 3.14  SE = 0.01643372 

估计的标准误差(SE)还是比较大, 提高到N=106,精度可以增加一位小数:

est_pi(1E6)
## N =  1e+06  pi估计值 = 3.145576  SE = 0.001639408 

40.3.2 线性回归模拟

考虑如下线性回归模型

y=10+2x+ε, ε∼N(0,0.52).

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

模型中的自变量x可以用随机数产生, 比如,用sample()函数从1:10中随机有放回地抽取n个。 模型中的随机误差项ε可以用rnorm()产生。 产生一组样本的程序如:

n <- 10; a <- 10; b <- 2
x <- sample(1:10, size=n, replace=TRUE)
eps <- rnorm(n, 0, 0.5)
y <- a + b * x + eps

如下程序计算线性回归:

lm(y ~ x)
## 
## Call:
## lm(formula = y ~ x)
## 
## Coefficients:
## (Intercept)            x  
##      10.133        2.036

如下程序计算线性回归的多种统计量:

summary(lm(y ~ x))
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34460 -0.17868 -0.02952  0.12797  0.49934 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.13297    0.16058   63.10 4.42e-12 ***
## x            2.03556    0.03258   62.49 4.78e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2819 on 8 degrees of freedom
## Multiple R-squared:  0.998,  Adjusted R-squared:  0.9977 
## F-statistic:  3905 on 1 and 8 DF,  p-value: 4.783e-12

如下程序返回一个矩阵, 包括a,b的估计值、标准误差、t检验统计量、检验p值:

summary(lm(y ~ x))$coefficients
##              Estimate Std. Error  t value     Pr(>|t|)
## (Intercept) 10.132970 0.16058073 63.10203 4.423167e-12
## x            2.035564 0.03257544 62.48768 4.782587e-12

如下程序把上述矩阵的前两列拉直成一个向量返回:

c(summary(lm(y ~ x))$coefficients[,1:2])
## [1] 10.13296972  2.03556360  0.16058073  0.03257544

这样得到 â ,b̂ ,SE(â ),SE(b̂ )这四个值。

replicate(, 复合语句)执行多次模拟, 返回向量或矩阵结果, 返回矩阵时,每列是一次模拟的结果。 下面是线性回归整个模拟程序,写成了一个函数。

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()
## 1000 次模拟的平均值:
##         a         b      SE.a      SE.b 
## 9.9970476 1.9994490 0.3639505 0.0592510 
## 1000 次模拟的标准差:
##          a          b       SE.a       SE.b 
## 0.37974881 0.06297733 0.11992515 0.01795926

可以看出,标准误差作为â ,b̂ 的标准差估计, 与多次模拟得到多个â ,b̂ 样本计算得到的标准差估计是比较接近的。 结果中SE(â )的平均值为0.364, 1000次模拟的â 的样本标准差为0.380,比较接近; SE(b̂ )的平均值为0.0593, 1000次模拟的b̂ 的样本标准差为0.0630,比较接近。

40.4 bootstrap

随机模拟研究中, 可以从理论模型生成许多组独立样本, 对每一组样本进行估计, 最后评估估计效果。

实际问题中, 仅有一组样本, 能否利用随机模拟方法评估对这组样本的估计效果? bootstrap就是这样一种方法, 对原始样本进行“再抽样”, 得到与原始样本分布相近的许多组样本, 称为bootstrap样本, 对每一个bootstrap样本进行估计, 用多个样本的估计结果来评估估计效果。

再抽样分为非参数和参数两种方法, 其中非参数bootstrap对原始数据集进行多次(如B=1000)重复有放回抽样, 每次得到一个同样大小的bootstrap数据集, 对每个bootstrap数据集进行相同的估计过程, 用得到的多次估计结果评估参数估计的精度、分布等。 参数方法是对原始数据集利用参数方法建模, 然后把估计出来的参数当作真实模型参数产生新的多组样本。

40.4.1 boot扩展包

基本R内嵌的boot包可以用来对数据集进行bootstrap抽样, 输入数据集和对数据集进行处理的函数, 直接输出对多个bootstrap数据集分析的结果。

作为示例, 考虑肺癌放疗研究中放疗前后的肿瘤大小(v0v1)的线性回归模型估计参数β̂ 0, β̂ 1的评估问题。 得到的回归方程是:

v1≈β̂ 0+β̂ 1v0.

library(tidyverse)
library(boot)
d.cancer <- read_csv(
  "data/cancer.csv", 
  locale=locale(encoding="GBK"),
  show_col_types=FALSE)
fa <- function(d, inds){
  dd <- d[inds,]
  lm1 <- lm(v1 ~ v0, data=dd)
  unname(coef(lm1))
}

lm0 <- lm(v1 ~ v0, data=d.cancer)
print(summary(lm0))
## 
## Call:
## lm(formula = v1 ~ v0, data = d.cancer)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.536 -10.582  -5.011  11.750  58.621 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.33079    7.26335   0.459     0.65    
## v0           0.37576    0.05376   6.990  6.4e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.56 on 32 degrees of freedom
## Multiple R-squared:  0.6042, Adjusted R-squared:  0.5919 
## F-statistic: 48.86 on 1 and 32 DF,  p-value: 6.401e-08
set.seed(101)
btr <- boot(d.cancer, fa, R=1000)
cat("Intercept: bias = ", btr$t0[1] - mean(btr$t[,1]),
    "  SE = ", sd(btr$t[,1]), "\n")
## Intercept: bias =  0.2877853   SE =  4.78422
cat("v0 slope: bias = ", btr$t0[2] - mean(btr$t[,2]),
    "  SE = ", sd(btr$t[,2]), "\n")
## v0 slope: bias =  -0.003612172   SE =  0.05450733

boot()函数的第一个自变量是data, 为原始数据集, 第二个自变量是statistic, 是统计估计函数, 该函数第一个自变量为原始数据, 第二个自变量是每次重抽样的下标向量, R是重复抽样次数, 函数输出要考察的参数估计。

boot()的结果(上面程序中的btr)是一个列表, 列表元素t0是对原始数据集的参数估计, 而t是对每一个重抽样数据集的参数估计, 每个重抽样的数据集占据矩阵t的一行, 每个估计量占据一列。

在这个例子中, bootstrap对斜率项的参数估计的标准误差估计与参数模型结果相近, 但对截距项的标准误差估计与参数模型结果相差较大。 还利用bootstrap方法估计了偏度, 但是注意线性回归模型的最小二乘估计在理论上是无偏估计。

利用得到的bootstrap样本计算β0和β1的置信区间:

ci_b0 <- quantile(btr$t[,1], c(0.025, 0.975))
ci_b1 <- quantile(btr$t[,2], c(0.025, 0.975))
cat("95% CI for beta0: (",
    ci_b0[1], ", ", ci_b0[2], ")\n")
## 95% CI for beta0: ( -5.944525 ,  12.44923 )
cat("95% CI for beta1: (",
    ci_b1[1], ", ", ci_b1[2], ")\n")
## 95% CI for beta1: ( 0.277491 ,  0.4920218 )

40.4.2 rsample扩展包

boot扩展包的功能比较固定, 另一个rsample扩展包则与函数式编程风格比较接近, 可以比较灵活地进行各种重抽样, 包括bootstrap、交叉验证、重要性抽样等。

bootstraps函数输入一个数据集, 输出多个bootstrap数据集并将其存放在一个tibble数据框的列表类型的列splits中, 每个元素是rsplit类型, 可以用as.data.frame()或者analysis()函数将其转换成数据框。

利用tibble的列表列,我们可以实现管道形式的bootstrap分析。 如:

library(rsample)
fb <- function(dd){
  lm1 <- lm(v1 ~ v0, data=dd)
  as_tibble(rbind(coef(lm1)))
}

btr2 <- d.cancer |>
  bootstraps(times=1000) |>
  pull(splits) |>
  map_dfr( \(df) fb(analysis(df)) )

cat("Intercept: bias = ", coef(lm0)[1] - mean(btr2[[1]]),
    "  SE = ", sd(btr2[[1]]), "\n")
## Intercept: bias =  0.1611144   SE =  4.905335
cat("v0 slope: bias = ", coef(lm0)[2] - mean(btr2[[2]]),
    "  SE = ", sd(btr2[[2]]), "\n")
## v0 slope: bias =  -0.001650511   SE =  0.05649049
ci_b0 <- quantile(btr2[[1]], c(0.025, 0.975))
ci_b1 <- quantile(btr2[[2]], c(0.025, 0.975))
cat("95% CI for beta0: (",
    ci_b0[1], ", ", ci_b0[2], ")\n")
## 95% CI for beta0: ( -6.250138 ,  12.52947 )
cat("95% CI for beta1: (",
    ci_b1[1], ", ", ci_b1[2], ")\n")
## 95% CI for beta1: ( 0.2667938 ,  0.4824034 )

结果的btr2是一个两列的数据框, 每一行包含对一个重抽样数据集的两个回归系数估计。

40.4.3 例:核密度的bootstrap置信区间

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

x <- faithful$eruptions
est0 <- density(x)
plot(est0)
北京大学R语言教程(李东风)第40章: 随机模拟

这个密度估计明显呈现出双峰形态。

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

重复B=10000次, 每次从原始样本中有重复地抽取与原来大小相同的一组样本, 对这组样本计算核密度估计, 结果为(xi,y(j)i),i=1,2,…,N,j=1,2,…,B, 每组样本估计N个格子点的密度曲线坐标, 横坐标不随样本改变。

对每个横坐标xi, 取bootstrap得到的B个y(j)i,j=1,2,…,B的0.025和0.975样本分位数, 作为真实密度f(xi)的bootstrap置信区间。

在R中利用replicate()函数实现:

set.seed(1)
resm <- replicate(10000, {
    x1 <- sample(x, replace=TRUE)
    density(x1, from=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='grey', border=FALSE)
lines(est0, lwd=2)
北京大学R语言教程(李东风)第40章: 随机模拟

程序中用set.seed(1)保证每次运行得到的结果是不变的, replicate()函数第一参数是重复模拟次数, 第二参数是复合语句, 这些语句是每次模拟要执行的计算。 在每次模拟中, 用带有replace=TRUE选项的sample()函数从样本中有放回地抽样得到一组bootstrap样本, 每次模拟的结果是在指定格子点上计算的核密度估计的纵坐标。 replicate()的结果为一个矩阵, 每一列是一次模拟得到的纵坐标集合。 对每个横坐标格子点,用quantile()函数计算B个bootstrap样本的2.5%和97.5%分位数, 循环用apply()函数表示。 polygon()函数指定一个多边形的顺序的顶点坐标用col=指定的颜色填充, 本程序中实现了置信下限与置信上限两条曲线之间的颜色填充。 lines()函数绘制了与原始样本对应的核密度估计曲线。

40.5 附录:完全试验方案设计函数

最简单的模拟试验设计, 是将考虑的所有因素的所有水平组合都执行, 比如,分别有2、3、4个水平的三个因素要进行完全试验, 需要2×3×4=24个方法。 下面的函数可以用来产生这样的试验方案, 输出为一个数据框, 数据框的每一行是一次试验的方法。

tidyr包的expand_grid()函数可以实现这样的功能, 基本R的expand.grid()函数有类似功能。 这两个函数输入因子,产生完全试验组合。 参见23.17

下面的程序直接实现了略强一些的功能, 支持以数据框为输入。

## 生成完全试验设计方案的函数。
## 输入:
##   varlist: 列表,元素可以是正整数、向量或数据框,
##     每个元素代表一个单纯的或者组合的因素。
##     元素取正整数时,表示仅用编号表示的因素的水平数。
##     元素取数值型或字符型向量时,表示因素的各个水平;
##     元素取多列的数据框时,每一行是一种因素组合,
##     数据框的每一列仍看作一个因素,但只允许数据框的各行出现的因素组合。
##   numbering: 真值表示输出中增加一列expNumber表示互不相同的试验方案的编号,
##     假值则不增加这一列。
##   size: 每种试验方案重复多少行。
## 输出:一个R数据框,每一行代表一次试验,每一列代表一个因素取值。
complete.design <- function(
  varlist,  
  numbering=TRUE,  # 是否为试验方案编号
  size=1){         # 每个方案的重复次数
  nf <- length(varlist)
  vl <- vector(mode="list", length=nf)
  tf <- character(nf)
  for(i in seq_along(varlist)){
    if(length(varlist[[i]])==1){
      if(varlist[[i]]==as.integer(varlist[[i]]) && varlist[[i]] > 0){
        # 用一个正整数表示的因子,就只需要用编码表示各个水平
        vl[[i]] <- factor(seq(varlist[[i]]))
        tf[i] <- "factor"
      } else {
        stop("非法输入!")
      }
    } else if(is.data.frame(varlist[[i]])){
      # 数据框,可以有一到多列,但是不可拆分,看作固定的组合
      vl[[i]] <- seq(nrow(varlist[[i]]))
      tf[i] <- "df"
    } else {
      # 普通向量,元素为因素水平
      stopifnot(!any(duplicated(varlist[[i]])))
      vl[[i]] <- factor(varlist[[i]])
      tf[i] <- "factor"
    }
  }
  
  ## 每个因素的水平数
  nlev <- vapply(vl, length, 1)
  ## 总试验次数
  ne <- prod(nlev) * size
  ## 每个因子的重复数
  nr <- rev(cumprod(rev(c(nlev, size))))[-1]
  ## 每个因子的循环次数
  nc <- ne / nr
  
  ## 产生作为数据集草稿的列表
  dl <- list()
  for(i in seq_along(vl)){
    idl <- as.integer(gl(nlev[i], nr[i], ne))
    if(tf[i]=="factor"){
      ## 单个因素
      dl[[i]] <- data.frame(
        factor(idl, labels=levels(vl[[i]])))
      names(dl[[i]]) <- names(varlist)[[i]]
    } else {
      ## 数据框,要恢复数据框
      dl[[i]] <- varlist[[i]][idl,]
    }
  }
  
  ## 合并数据框得到最终结果
  d <- dl[[1]]
  if(nf>1){
    for(i in seq(2, nf)){
      d <- cbind(d, dl[[i]])
    }
  }
  
  if(numbering){
    d <- cbind(d, data.frame(
      expNumber = rep(seq(ne/size), each=size)))
  }
  
  d
}


## 测试程序
test1 <- function(){
  complete.design(list(f1=3, f2=c("a", "b"), f3=4))
}
print(test1())

test2 <- function(){
  complete.design(list(f1=3, f2=c("a", "b"), f3=4), size=2)
}
print(test2())

test3 <- function(){
  complete.design(list(f1=3, f2=c("a", "b"), f3=data.frame(
    fa=c(0.01, 0.05, 0.1), fb = c("x", "x", "y"))), size=2)
}
print(test3())

韭菜热线原创版权所有,发布者:风生水起,转载请注明出处:https://www.9crx.com/79941.html

(0)
打赏
风生水起的头像风生水起普通用户
上一篇 2023年11月29日 23:32
下一篇 2023年11月30日 23:07

相关推荐

  • 北京大学R语言教程(李东风)第43章:基于树的方法

    43.1 树回归的简单演示 决策树方法按不同自变量的不同值, 分层地把训练集分组。 每层使用一个变量, 所以这样的分组构成一个二叉树表示。 为了预测一个观测的类归属, 找到它所属的组, 用组的类归属或大多数观测的类归属进行预测。 这样的方法称为决策树(decision tree)。 决策树方法既可以用于判别问题, 也可以用于回归问题,称为回归树。 …

    2023年12月2日
    18800
  • 新兴市场和美元贬值:全球风险资产将受益于美元疲软

    如果您不是货币交易者或去年没有出国旅行,您可能不知道美元兑其他货币已跌至 15 个月来的最低水平。随着美联储加息周期接近尾声,世界金融体系的关键货币的贬值只会加速。在这种环境下,多种资产都可以获利,其中最重要的是注入人工智能的美国科技行业。然而,没有什么地方比新兴市场更容易从这种环境中受益。 新兴市场的崛起 从历史上看,新兴市场(EM)股票与美元呈反比关系,…

    2023年9月15日
    11000
  • 北京大学R语言教程第58章:R编程例子(附:第59章:使用经验)

    因为这些例子要用作学生习题,所以这里只有问题,没有解答。 R语言 用向量作逆变换 设向量x长度为n,其中保存了1到n的正整数的一个排列。把x看成是在集合{1,2,…,n}上的一个一一变换,求向量y使得y能够表示上述变换的逆变换。即任给长度为n的向量z, z[x]表示按照x的次序重新排列z的元素,而z[x][y]则应该恢复为z。 斐波那契数列计算 设数列x0=…

    2023年12月22日
    13300
  • 什么事信用管理?信用管理的目标是什么?

    您如何管理信用可以决定您的个人财务状况。信贷可以是获得您需要和想要的东西的有用工具,但如果您不小心,它也可能导致您的财务崩溃。 明智地使用信贷可以为您提供一生的机会,但滥用信贷或积累您无法偿还的债务可能会损害您的经济利益,并关闭您可能没有考虑过的大门。无论是由于年轻、经验不足、缺乏知识还是个人财务危机,许多人随着时间的推移犯下了令人遗憾的财务失误,并发现自己…

    2023年7月10日
    13400
  • 分析疫情过后的美国经济

    作者:尤金尼奥·阿莱曼 (Eugenio Aleman)、Giampiero Fuentes,24/4/5 首席经济学家 Eugenio Alemán 和经济学家 Giampiero Fuentes 指出,虽然他们预计经济增长将放缓,但他们并不认为 2024 年会出现经济衰退。 我们通常不愿意用时髦的短语来解释我们对美国经济的好坏判断。然而,今天说“这一次不…

    2024年5月16日
    3300

发表回复

登录后才能评论
客服
客服
关注订阅号
关注订阅号
分享本页
返回顶部