Randu в R: заполнить матрицу без цикл

Question

Возьмем следующую реализацию плохо псевдослучайных чисел генератора (ПСЧ) 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) 
} 

Я хочу, чтобы заменить цикл. Проблема здесь в том, что записи целой строки, которые генерируются с каждым шагом итерации, необходимы для последующего шага итерации. Моя идея - применить функцию для заполнения строк пустой матрицы. Так что я пытаюсь стать знаком с функцией 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) 

..which не очень много. Я не могу понять, как перебирать каждый элемент строки и как генерировать, например. 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