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 hope I won't get spanked for reposting. This is related to a previous question discussed here: Performance on Multinomial Deviate is slow.

but this version is simpler because it is a binomial rather than multinomial which is my justification for the new question.

I've kluged together a real solution for this problem that someone might find useful. I'm hoping someone will improve on the readability and/or performance. The performance is acceptable but still a bit slow compared to C. I think there might be another factor of 2 or 3 speedup to be had. My question is: does anyone have suggestions that will improve on my solution?

The problem is that I need a bunch of deviates where the parameters of the distribution are different for each deviate. Here is the simple way to do it which is unacceptably slow.

Nboxes = 8000;
p = RandomReal[{0, .01}, Nboxes];
rt = RandomInteger[{0, 10}, Nboxes];
rp = RandomInteger[{0, 10}, Nboxes];
probs = 1 - p rp;
out = MapThread[
    RandomVariate[BinomialDistribution[#1, #2]] &, {rt, probs}]; // Timing

The timing result was over 4.7 seconds on my laptop. I create one below that is about 25 times faster.

Now to generate a binomial deviate with parameters n and p you generate n uniform deviates on (0,1) and count how many are < p. The total number of uniform deviates needed then is Total[rt]. Those can be generated in a small fraction of a second. The following code implements a binomial deviate generator and runs in a total of under 0.2 seconds. The deviates are stored in "counts". If anyone can improve on this code for performance, readability, etc, I'm all ears. First I generate all the uniform deviates I need. Then I "box" those up according to what rt is and count the number of deviates in each partition that are less than that boxes p value.

unitdeviates = RandomReal[1, Total[rt]]; // Timing
nt = ConstantArray[0, Nboxes + 1];
Table[nt[[i]] = nt[[i - 1]] + rt[[i - 1]], {i, 2, Length[rt] + 1}]; // Timing

counts = ConstantArray[0, Nboxes];
Do[counts[[nboxes]] = Count[unitdeviates[[nt[[nboxes]] + 1 ;; 
   nt[[nboxes + 1]]]], _?(# < probs[[nboxes]] &)], {nboxes, 1,Nboxes}] // Timing

Here I show that the mean and variance of the two sets of deviates are consistent

N[Mean[out]]
N[Mean[counts]]
N[Variance[out]]
N[Variance[counts]]
share|improve this question
    
Sorry, but I don't see a question here. –  m_goldberg Mar 29 at 1:08
    
The question is can the provided code be improved upon either for readability or for performance. –  JEP Mar 29 at 2:13
    
These deviates can be generated in C in about 5 ms---roughly 40 times faster. I'm not sure what performance hit I should expect in Mathematica over C but a factor of 40 strikes me as excessive. A call like this: RandomVariate[BinomialDistribution[5, .2], 8000]; in which the parameters don't change takes about 2 ms which is even faster than my C code. –  JEP Mar 29 at 6:27
1  
Closers, thanks for keeping this open (and not calling us C-losers? :P). Anyway I think it is nice to close questions that are unlikely to be unanswered, but that reason to close this question should be gone now. Let's see if I can get the accept :). –  Jacob Akkerboom Mar 31 at 19:15

1 Answer 1

The main thing that I am trying to show is that you can use Accumulate and that almost all these functions are compilable. I hope it also shows when to use Table rather than Do, to avoid making unnecessary ConstantArrays. I personally find the use of Table in your code confusing. Of course it is nice to localise variables from time to time, which is also done in this code. All in all it is pretty boring code and maybe it looks a lot like what you have written in C, but I hope you learn something anyway.

In the code below, there is quite a lot of code to define a single function. There are multiple ways to avoid this, for example calling other CompiledFunctions inside the CompiledFunction. Another way is to inline definitions, which can also be done in multiple ways. For simplicity, I have just dumped everything in one definition.

cfu =
 Compile[{{nBoxes, _Integer, 0}},
  Block[
   {unifs, rt,  min, max, acc, count, p, rp, probs}
   ,
   rt = RandomInteger[{0, 10}, nBoxes];

   acc = Accumulate[rt];

   unifs = RandomReal[1., acc[[-1]]];

   p = RandomReal[{0., .01}, nBoxes];
   rp = RandomInteger[{0, 10}, nBoxes];
   probs = 1 - p rp;

   min = 1;

   Table[
    max = acc[[nb]];

    count = 0;
    Do[
     If[
      unifs[[j]] < probs[[nb]]
      ,
      count++
      ]
     ,
     {j, min, max}
     ];

    min = max + 1;
    count
    ,
    {nb, 1, nBoxes}
    ]

   ]
  ,
  CompilationTarget -> "C"
  ]

Let's also make definitions to compare this with your code

yourInit :=
 (
  Nboxes = 8000;
  p = RandomReal[{0, .01}, Nboxes];
  rt = RandomInteger[{0, 10}, Nboxes];
  rp = RandomInteger[{0, 10}, Nboxes];
  probs = 1 - p rp;
  )

yours1 :=
  (
   yourInit;
   out = MapThread[
     RandomVariate[BinomialDistribution[#1, #2]] &, {rt, probs}];
   );

yours2 :=
 (
  yourInit;
  unitdeviates = RandomReal[1, Total[rt]];
  nt = ConstantArray[0, Nboxes + 1];
  Table[nt[[i]] = nt[[i - 1]] + rt[[i - 1]], {i, 2, Length[rt] + 1}];

  counts = ConstantArray[0, Nboxes];
  Do[counts[[nboxes]] = 
    Count[unitdeviates[[nt[[nboxes]] + 1 ;; 
        nt[[nboxes + 1]]]], _?(# < probs[[nboxes]] &)], {nboxes, 1, 
    Nboxes}]
  )

We then have

yours1 // Timing // First
Mean[out] // N
Variance[out] // N
3.770396  
 4.85013  
 9.70163  
yours2 // Timing // First
Mean[counts] // N
Variance[counts] // N
0.132868  
 4.8515  
 9.53689
(res = cfu[8000]) // Timing // First
Mean@res // N
Variance@res // N
0.002822  
 4.87888  
 9.66966

So the speedup is similar to your own C code. Using LibraryLink we should be able to do a little bit better.

share|improve this answer
    
Your code is much easier to follow than mine. I didn't know about Accumulate. My C code is ranlib which I highly recommend. netlib.org/random . I hadn't understood how to compile functions with lists, arrays, etc as arguments but now I know! I moved rp, rt, etc into the argument list. I am still a little hazy as to what the function returns. I guess it returns count one element at a time? –  JEP Apr 1 at 18:42
    
@JEP thank you for the recommendation, I would like to explore some C libraries. Feel free to move rp rt etc to the argument list, that should hardly matter performance wise. The structure of my code should be quite similar to what I call yours2. What the compiled function returns is what you call counts or out, hence the comparison of Means/Variances. For each value of nb (nboxes), the number of uniform variables that is smaller than probs[[nb]] (count) is calculated. min = max + 1 is just setting things up for the next iteration, it is count that is put in the list. –  Jacob Akkerboom Apr 1 at 22:55

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.