Performance on Binomial Deviate is slow

Question

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]]

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.