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