Draw non full matrix of random numbers

  • A+

I am doing a Monte-Carlo simulation, where each repetition requires the sum or product of a random number of random variables. My problem is how to do this efficiently as the entire simulation should be as vectorized as possible.

For example, say we want to take the sum of 5, 10 and 3 random numbers, represented by the vector len = [5;10;3]. Then what I am currently doing is drawing a full matrix of random numbers:

A = randn(length(len),max(len)); 

Creating a mask of the non-needed numbers:

lenlen = repmat(len,1,max(len)); idx = repmat(1:max(len),length(len),1); mask = idx>lenlen; 

and then I can "pad", the matrix as I am interested in the sum the padding have to be zero (for the case with the product the padding had to be 1)


To obtain:

A =  1.7708   -1.4609   -1.5637   -0.0340    0.9796         0         0         0         0         0 1.8034   -1.5467    0.3938    0.8777    0.6813    1.0594   -0.3469    1.7472   -0.4697   -0.3635 1.5937   -0.1170    1.5629         0         0         0         0         0         0         0 

Whereafter I can sum them together

B = sum(A,2); 

However, I find it rather superfluous that I have to draw too many random numbers and then throw them away. In the real case, I need in the range of hundred thousands of repetitions and the vector len might vary a lot, i.e. it can easily be that I have to draw twice or three times the number of random numbers than of what is needed.

You can generate the exact amount of random numbers required, create a grouping variable with repelem, and compute the sum of each group using accumarray:

len = [5; 10; 3]; B = accumarray(repelem(1:numel(len), len).', randn(sum(len),1)); 


:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: