Fast way to generate random value from a distribution in R based on different parameters inputted from very large matrices

user_19261284 0
Viewed 40

I need to generate a LOT of random numbers from a negative binomial distribution, all with different distribution parameters (the parameters I am using are mu and size for rnbinom). I need a faster way of doing this.

So right now I have 3 very large matrices: mean, exp and size. All of them have the same dimension and distinct numbers. I need to return a matrix M, where M(i, j) being a number randomly generated with NegativeBinomial(mu=mean(i, j)^exp(i, j), size=size(i, j)). Here is my current way of doing it:

f <- function(mean, exp, size) rnbinom(1, mu=mean^exp, size=l2fc)
M <- structure(mapply(f, mean, exp, size), dim=dim(mean))

From what I can tell mapply is already about as fast as it can be with R's internal vector optimization. Is there anyway to optimize this even further? The dim of mean/exp/size is roughly about 10000*25000, and it takes as long as an hour to finish running the entire code.

r functional-programming
1 Answers

rnbinom() is vectorised over the distributional parameters. You can get the draws with one call:

set.seed(42)

mu <- matrix(1:4, 2, 2)
size <- matrix(1:4, 2, 2)

n <- prod(dim(mu))
M <- rnbinom(n, size = size, mu = mu)
dim(M) <- dim(mu)

M
#>      [,1] [,2]
#> [1,]    1    4
#> [2,]    3    6

A quick performance comparison with more parameters shows ~30x improvement:

mu <- matrix(1:4, 2000, 2000)
size <- matrix(1:4, 2000, 2000)

set.seed(42)
system.time({
  f <- function(mu, size) rnbinom(1, mu=mu, size=size)
  M1 <- structure(mapply(f, mu, size), dim=dim(mu))
})
#>    user  system elapsed 
#>   29.36    0.30   30.11

set.seed(42)
system.time({
  n <- prod(dim(mu))
  M2 <- rnbinom(n, size = size, mu = mu)
  dim(M2) <- dim(mu)
})
#>    user  system elapsed 
#>    1.08    0.03    1.13

identical(M1, M2)
#> [1] TRUE
user_45506950
Related