R中的randu：在没有for循环的情况下填充矩阵

Question

采取以下bad伪随机数生成器（PRNG）randu的实现：

n = 100000
randu = matrix(NA, ncol=3, nrow=n)
new_z = 1

for(i in 1:n) {
    new_x = (65539*new_z) %% 2^31
    new_y = (65539*new_x) %% 2^31
    new_z = (65539*new_y) %% 2^31
    randu[i,] = c(x=new_x/2^31, y=new_y/2^31,z=new_z/2^31)
}   

我想替换 for 循环。这里的问题是，在每个迭代步骤中生成的整行条目对于后续迭代步骤是必需的。我的想法是应用一个函数来填充空矩阵的行。因此，我正在尝试熟悉 apply 函数，而我已经走到了这一步：

n = 100000
randu = matrix(NA, ncol=3, nrow=n)
randu[1,3] <- 1 # seed
randu.fct <- function() {
  randu[,1] <- (65539 * randu[,3]) %% 2 ^ 31
  randu[,2] <- (65539 * randu[,1]) %% 2 ^ 31
  randu[,3] <- (65539 * randu[,2]) %% 2 ^ 31 
}
apply(randu[,1:3],1,randu.fct)

..这不是很多。我无法理解如何迭代一行的每个元素以及如何生成例如100000 行。

Answer 1

你可以用replicate替换（隐藏）循环：

n = 100000
x = 1
matrix(replicate(3*n, {x <<- (65539*x) %% 2^31})/2^31, ncol = 3, byrow = TRUE)

如果你需要速度，你可能应该看看Rcpp，下面是一个简单的实现：

library(inline)
library(Rcpp)

cppFunction(
    'NumericMatrix randU(int n) {
        NumericMatrix X(n, 3);
        int x = 1;
        unsigned int d = 2147483648;
        for (int i = 0; i < n; ++i) {
          for (int j = 0; j < 3; ++j) {
               x = (65539*x) % d;
               X(i,j) = x  / double(d);
          }
        }
        return X;
    }'
)

> all.equal(randu, randU(100000))
[1] TRUE

一个小的速度比较：

f1 <- function(){
    n = 100000
    randu = matrix(NA, ncol=3, nrow=n)
    new_z = 1
    for(i in 1:n) {
        new_x = (65539*new_z) %% 2^31
        new_y = (65539*new_x) %% 2^31
        new_z = (65539*new_y) %% 2^31
        randu[i,] = c(x=new_x/2^31, y=new_y/2^31,z=new_z/2^31)
    }  
    randu
}

f2 <- function(){
    n = 100000
    x = 1
    matrix(replicate(3*n, {x <<- (65539*x) %% 2^31})/2^31, ncol = 3, byrow = TRUE)
}

f3 <- function(){
    randU(100000)
}


Unit: milliseconds
 expr         min          lq        mean      median          uq       max neval
 f1() 1170.166889 1245.987545 1331.918328 1320.593903 1356.121828 1593.0860    10
 f2() 1194.103998 1449.195295 1499.362126 1514.794140 1562.868296 1798.1218    10
 f3()    2.041235    2.055671    3.386515    2.207969    2.676895   13.1357    10