It took two days and required me to think which is, of course, disgusting, but I found a solution. The original pseudo-code loop (see above) was a loop over things that I think of as "boxes" (not mentioned in the original pseudo-code). The number of boxes has typically been 8,000 but I can envision problems for which it might be some millions. What I needed was essentially the number of zeros in each returned multinomial variate. For example if the variate was (1,2,0,0,5) the result to be stored would be "2" because there are two zeros in that variate. (THis is slightly neutered for simplicity.)
Now instead of running a single loop over all boxes I run a double loop over all (np, nt) pairs. The variable "boxes" contains a list of the indices for the current (np,nt) pair. The variates (called "out" above) are stored in er, "variates" and the random number generator returns Length[boxes] results instead of just 1 as in the original version. The results are stored in results[[boxes]] . I haven't done any timing yet but the original loop (over boxes) took about 5-6 seconds to run and the double loop version runs in a small fraction of a second.
Create lists of all np and nt values and store in "nps" and "nts".
Create function vec[np] which returns the probability vector (called tmp in the original).
Initialize result=ConstantArray[0, Nboxes]
Loop over all (np,nt) pairs (up to Max[nps] and Max[nts]):
boxes = Flatten[Position[Transpose[{nps, nts}], {np, nt}]];
variates =
RandomVariate[MultinomialDistribution[nt, vec[np]],
Length[boxes]];
result[[boxes]] = Map[Count[#, 0] &, variates];
Arg. I don't know if these things get bumped when they're edited. In any case I've realized that I made an error in my previous solution. I am more or less back where I started. I will rewrite a clearer version of the original pseudo-code below.
Loop over boxes from 1 to Nboxes;
Get new np and nt values
vec=vec[np,boxes];
out = RandomVariate[MultinomialDistribution[nt, vec[np,boxes]]];
Do stuff with out
End Loop
The new wrinkle here is that vec depends on the box index. The solution above assumed that vec=vec[np] but not vec=vec[np,box]. The box-dependence of vec kills that solution. Is there a way to compile the random variate code that will improve performance? Most of my tests are on cases where Nboxes=8000 and it takes 5-6 seconds to run through the loop.
I'd like to be able to increase Nboxes to millions. Anyway, if this code were in C there wouldn't be a penalty for the repeated calls to the RandomVariate thing. I don't know how to get this to even compile currently.
RandomVariate
for results one at a time is a bit brain-dead. If you're going to need a set of variates, generate them at one go and use that result. E.g., for your example of 10K variates, there's a ~500 to 1 speed difference. $\endgroup$For
. It's generally better to use functional constructs likeTable
. When you can't, useDo
. It's more readable thanFor
, less error prone (iterator is localized), and marginally faster. Instead ofTimeUsed
, use the functionTiming
orAbsoluteTiming
. $\endgroup$