Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
You can speed your code up by changing p = Flatten[Position[l, s]]; l = Drop[l, p] to l=DeleteCases[l,s]. –  b.gatessucks Apr 5 '13 at 20:20
    
thank you, I'll keep that in mind –  Emmet Apr 5 '13 at 20:27
add comment

2 Answers

up vote 4 down vote accepted

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

share|improve this answer
    
Thank you very much, So basically you do a RandomSample then check at what position is v, right? –  Emmet Apr 5 '13 at 20:18
add comment

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.

share|improve this answer
    
I'll sure do, much obliged. –  Emmet Apr 5 '13 at 20:34
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.