简化 R 上的仿真

Question

正如我在上一个问题中提到的。我是编程新手，之前没有任何经验，但很高兴能够学习。但是，我运行遇到了以下问题，我的教授给了我们以下内容：

sim1 <- function(n) {
  xm <- matrix(nrow=n,ncol=2)
  for (i in 1:n) {
    d <- rnorm(1)
    if (runif(1) < 0.5) {
      xm[i,1] <- 1
      xm[i,2] <- 2.5*d + 69
    } else {
      xm[i,1] <- 0
      xm[i,2] <- 2*d + 64
    }
  }
  return(xm)
}

有以下任务：尝试提高这段代码的效率。使用 speed.test 查看它是否在生成 n=1000 个观测值方面有所改进。

我终于至少能够弄清楚这段代码的作用，尽管如此，我完全不知道如何才能使这段代码更有效率。

任何帮助都意义重大。谢谢！

Answer 1

我建议的一个优化是创建默认值为 0 的矩阵。一旦以 0 值作为默认值创建了矩阵，则无需在函数中填充值 0。

修改后的代码如下所示：

sim1 <- function(n) {
#create matrix with 0 value. 
xm <- matrix(0,nrow=n,ncol=2) 
for (i in 1:n) {
d <- rnorm(1)
if (runif(1) < 0.5) {
 xm[i,1] <- 1
 xm[i,2] <- 2.5*d + 69
} else {
 #xm[i,1] <- 0    --- No longer needed
 xm[i,2] <- 2*d + 64
}
}
return(xm)
}

Answer 2

我将执行我认为最明显的步骤，即将 rnorm() 移出循环并利用其矢量化特性（如 rawr 所提到的）

sim2 <- function(n) {
    xm <- matrix(nrow=n, ncol=2)
    d <- rnorm(n)
    for (i in 1:n) {
        if (runif(1) < 0.5) {
            xm[i,1] <- 1
            xm[i,2] <- 2.5*d[i] + 69
        } else {
            xm[i,1] <- 0
            xm[i,2] <- 2*d[i] + 64
        }
    }
    return(xm)
}

n <- 1e3
set.seed(1); system.time(s1 <- sim1(n)); system.time(s2 <- sim2(n))
#  user  system elapsed 
# 0.019   0.004   0.023 
#  user  system elapsed 
# 0.010   0.000   0.009 

t.test(s1[,2], s2[,2]) # Not identical, but similar, again alluded to by rawr

这给了我们合理的改进。 runif() 也可以做类似的事情，但我会把它留给你。

如果你想阅读 material 我可以推荐 Hadley Wickhams Advanced R 和章节 Optimising code。

如果您想知道，确实可以同时消除循环和条件。

Answer 3

如果可能，不要在 R 中使用循环。rep 和 rnorm 将在一次调用中非常快速地用 5、10 或 500,000 个值填充向量。调用 rnorm(1) 500,000 次是一种浪费，而且比简单地调用 rnorm(500000) 慢得多。这就像开着法拉利兜风，走 1 英尺然后停下，走 1 英尺又停下，一遍又一遍地到达目的地。

此函数将 return 统计上与您的函数相同的结果。但是，它没有使用循环，而是以 R 方式执行操作。

sim2 <- function(n) {
    n1 <- floor(n/2)  #this is how many of the else clause we'll do
    n2 <- n - n1  #this is how many of the if clause we'll do
    col11 <- rep(0, n1) #bam! we have a vector filled with 0s
    col12 <- (rnorm(n1) * 2) + 64 #bam! vector filled with deviates
    col21 <- rep(1, n2) #bam! vector filled with 1s
    col22 <- (rnorm(n2) * 2.5) + 69 #bam! vector filled with deviates
    xm <- cbind(c(col11,col21), c(col12,col22)) #now we have a matrix, 2 cols, n rows
    return(xm[sample(nrow(xm)),]) #shuffle the rows, return matrix
}

没有循环！功能可能很明显，但如果不是，我会解释。首先，n1 & n2 只是将n 的大小适当拆分（占奇数）。

接下来，可以消除每个元素的二项式过程（即if(runif(1) < 0.5) {} else {}），因为我们知道在sim1中，矩阵的一半落入if条件，另一半在 else 中（见下面的证明）。当我们知道它是 50/50 时，我们不需要为每个元素一遍又一遍地决定采用哪条随机路径。因此，我们将首先完成所有 else 50%：我们用 n/2 0（col11）填充一个向量，另一个用 n/2 随机偏差（均值 = 0，默认情况下 sd = 1），并且对于每个偏差，乘以 2 并加上 64，结果向量 col12。完成了 50%。

接下来，我们完成第二个 50%（if 部分）。我们用 n/2 1s (col21) 填充一个向量，用随机偏差填充另一个向量，对于每个偏差，乘以 2.5 并加 69。

我们现在有 4 个向量，我们将把它们变成一个矩阵。第 1 步：我们使用 c 函数将 col11（填充有 n/2 0）和 col21（填充有 n/2 1）粘合在一起以获得向量（ n 个元素）。第 2 步：使用 c 将 col12 和 col22 粘合在一起（填充偏差）以获得向量（如 1 列 x n 行矩阵）。注意：0s/1s 与基于 64/69 公式的正确偏差相关联。第 3 步：使用 cbind 从向量中创建一个矩阵 (xm)：0/1 向量成为第 1 列，偏差向量成为第 2 列。第 4 步：获取矩阵中的行数（应该只是 n）使用 nrow。第 5 步：使用 sample 随机排列所有行号，制作一个打乱的向量。第 6 步：创建一个新的（未命名的）矩阵，根据打乱后的向量按顺序排列 xm 的行。步骤 4-6 的要点只是对行进行随机排序，因为 sim1 中的二项式过程会产生随机的行排序。

此版本运行速度提高了 866%！

> system.time({ sim1(500000)})
   user  system elapsed 
  1.341   0.179   1.527 
> system.time({ sim2(500000)})
   user  system elapsed 
  0.145   0.011   0.158

如果您担心这是否能保持二项式过程的完整性，请考虑二项式过程做两件事：1) 它将 1 与 2.5*d+69 方程相关联，将 0 与 2*d + 64 等式 - 由于行被原封不动地打乱，因此关联得以维持； 2) if 子句中有 50%，else 子句中有 50%，如下所示。

sim3 <- function(n) {
    a <- 0
    for(j in 1:n) {
        if(runif(1) < 0.5) {
            a <- a + 1
        }
    }
    return(a/n)
}
> sim3(50)
[1] 0.46
> sim3(5000)
[1] 0.4926
> sim3(10000)
[1] 0.5022
> sim3(5000000)
[1] 0.4997844

二项式过程产生 50% 的 1 和 50% 的 0（第 1 列）。

简化 R 上的仿真

Simplify Simulations on R

r

if-statement

simulation

simplification