Speeding up a simple simulation with random extractions

Question

I wonder if anybody can help, I need a hand with a simple simulation I am writing with Mathematica (I am using version 8). Basically it creates a list (called l) with numbers from 1 to n and a random number v in the same interval. Then it draws a random number s from l and checks if it is equal to v. If s is not equal to v, that number is dropped from l and another s is drawn and so on, the cycle stops when s=v. The number i is used to count the draws and the list k to show all the i's. Here is the code with n=5000 and 5 runs

n = 5000;
runs = 5;
(* initialization *)
k = ConstantArray[0, runs];
Table[
 {v = RandomInteger[{1, n}];
  s = 0;
  i = 0;
  l = Range[1, n];
  While[s != v, 
   s = RandomChoice[l];
   p = Flatten[Position[l, s]];
   l = Drop[l, p];
   i = i + 1;];
  k[[j]] = i}
 , {j, 1, runs}]
k

The problem is I would like to run thousands of simulations with n raging from 30000 to 50000, but this code is way too slow. Is there a way to speed it up?

You can speed your code up by changing p = Flatten[Position[l, s]]; l = Drop[l, p] to l=DeleteCases[l,s]. — b.gates.you.know.what, Apr 05 '13 at 20:20

score 5 · Accepted Answer · edited Apr 13 '17 at 12:56

If I understand the problem correctly, sampling l without replacement can be done all at once with RandomSample[l]. It is equivalent to sampling the entire list, past when you would quit, but it does it very quickly. Position tells you how many times it would take to sample v on that particular run.

n = 5000;
runs = 5;
l = Range[n];
k = Table[v = RandomInteger[{1, n}];
  s = RandomSample[l];
  Position[s, v, 1][[1, 1]], {j, 1, runs}]

It runs about 1000 times faster than your code. For n = 50000 on my machine, it's about 0.01 sec. per run.

P.S. You might be interested in this answer, which collects a lot links to programming tips for Mathematica

Thank you very much, So basically you do a RandomSample then check at what position is v, right? — Emmet, Apr 05 '13 at 20:18

ciao · Answer 2 · 2014-05-08T00:36:40.123

Late answer, but perhaps useful to future visitors wanting the same kind of simulation. No need to go through machinations of code and/or compile to get performance. Use e.g.:

n = 50000;
runs = 5000;

results=RandomVariate[BetaBinomialDistribution[1, 1, n - 1], runs] + 1;

This can be hundreds of times faster than the above answers.

If the end result desired is just the statistics of the results, no need to even generate results, just operate directly on the distribution, e.g.:

Mean[BetaBinomialDistribution[1, 1, n - 1]]+1

If, on the other hand, the results are all you want (perhaps for some machinations beyond basic statistics of the distribution), we can note:

FullSimplify[PDF[BetaBinomialDistribution[1, 1, n - 1], x] ==
             PDF[DiscreteUniformDistribution[{0, n - 1}], x]]

(* True *)

In other words, if just generating a list of results of number of tries until the target was "found" is desired,

results=RandomInteger[{1,n},runs]

will suffice, and runs as quickly as is probably possible (n=50K, runs=5K is below timing resolution on a netbook...)

score 2 · Answer 3 · answered Apr 05 '13 at 20:34

Another way to get a speedup is to compile. Your own code with a minimal change would look like :

makeRuns = 
Compile[{{n, _Integer}, {runs, _Integer}},
   Module[{k = Array[0 &, runs], v, s, i, l}, 
     Table[{v = RandomInteger[{1, n}];
     s = 0;
     i = 0;
     l = Range[1, n];
     While[s != v, s = RandomChoice[l];
      l = DeleteCases[l, s];
      i = i + 1;];
      k[[j]] = i}, {j, 1, runs}];
      k
   ], CompilationTarget :> "C", RuntimeOptions -> "Speed"
]

and it gives you a 10x boost. You can go ahead and try that with the other answers you have.

Speeding up a simple simulation with random extractions

3 Answers3