Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Create vector of random numbers (size of vector not known at run-time)

Tags:

r

I'm programming a Monte-Carlo Simulation that should give the user quite some flexibility. Consequently, I want the user to be able to specify the concrete probability distribution of the random numbers before the simulation is run. However, at this time it is not known to the user how many random numbers are needed.

My idea is now to get a call-object from user that creates ONE random number and then create as many of those random numbers as needed internally. However, except for a loop I can't get any other solution to work, but have the feeling that this is because of me missing something. So basically, I have two questions:

1) Is the idea with the call-object a good one? I'm still working on the project so I could still change the set-up, but I need a very intuitive solution for the user.

2) If it is a good idea, is there are more elegant way to expand the random number to a vector of size nrMCS?

Let's make an example:

#That's what I would get from the user with my current set-up:
rnd_call <- call("rnorm", 1, mean=0.1, sd=0.01)
#To create nrMCS random numbers, that's my best shot so far:
nrMCS <- 100
rnd_vec <- as.numeric(nrMCS)
for (i in 1:nrMCS){rnd_vec[i] <- eval(rnd_call)}
rnd_vec
[1] 0.09695170 0.11752132 0.11548925 0.11205948 0.10657986 0.12017120 0.09518435
...
#Question: Is there are more elegant way?
#I tried the following, but it fails for certain reasons
rep(eval(rnd_call), nrMCS) #DOES NOT WORK: Repeats ONE random number
[1] 0.1105464 0.1105464 0.1105464 0.1105464 0.1105464 0.1105464 0.1105464 0.1105464 
...
eval(rep(rnd_call, nrMCS)) #DOES NOT WORK
Error in rnorm(1, mean = 0.1, sd = 0.01, rnorm, 1, mean = 0.1, sd = 0.01,  : 
formal argument "mean" matched by multiple actual arguments
like image 591
Christoph_J Avatar asked Oct 24 '11 19:10

Christoph_J


1 Answers

I think a more idiomatic way to do this would be to take an r* function and a list of arguments. Whenever you can avoid calling eval, you should. Something like this:

rnd_fun <- rnorm
rnd_args <- list(mean=0.1,sd=0.01)
nrMCS <- 100
rnd_vec <- do.call(rnd_fun,c(list(n=nrMCS),rnd_args))

(this relies on the convention in R that the first argument to an r* (random deviate generator) function is always n, the number of deviates required ...)

Furthermore, calling rnd_fun once with n=nrMCS is generally much more efficient than calling it nrMCS times ...

library(rbenchmark)
nrMCS <- 10000
benchmark(single_call=do.call(rnd_fun,c(list(n=nrMCS),rnd_args)),
           mult_call=replicate(nrMCS,do.call(rnd_fun,c(list(n=1),rnd_args))))
         test replications elapsed relative user.self sys.self 
2   mult_call          100  11.135 91.27049    11.084    0.004 
1 single_call          100   0.122  1.00000     0.080    0.036
like image 86
Ben Bolker Avatar answered Sep 30 '22 13:09

Ben Bolker