I have a list of probabilities $p_k$, each corresponding to a Bernoulli distribution. These values of $p_k$ are all small (let's say 0.1% or less). I want to prepare a data structure that allows quickly sampling each of the Bernoulli distributions and returning the sparse set of them that returned 1 instead of 0.

Naively, the expected runtime cost of this sampling task would be $O(N)$ where $N$ is the number of probabilities. I want the expected runtime of producing another set of sparse samples to instead be like $\tilde{O}(pN)$ where $p$ is the maximum probability $p_k$.

I know that this sort of thing can be done, I just can't remember the name of the method. The method involves converting the probabilities into spans within an exponential distribution, creating a search tree over the spans, then using the tree to quickly identify which span a sample from that exponential distribution landed inside (which tells you the next value in the sparse sample), then recursing on the spans past the one that was landed in.

(An alternative method would be to group the probabilities into groups of size $1/p$, precomputing the probability distribution for the number of True's that appear when sampling all elements within the group, then using alias sampling to sample how many samples-without-repetition to extract [also using alias sampling].)

It would be particularly useful if a method achieving this kind of speedup was available in a common python library such as scipy.

  • $\begingroup$ The number of 0's between two consecutive 1's in your data set has a geometric distribution; you could just sample that instead, to generate a sparse vector directly. $\endgroup$ – Federico Poloni Jun 4 at 10:24
  • $\begingroup$ @FedericoPoloni Yes, that is the essence of the method I was describing. The complicating factor is that the probabilities $p_k$ are not all identical. $\endgroup$ – Craig Gidney Jun 4 at 18:15
  • 1
    $\begingroup$ This seems similar in principle to the Gillespie method in chemical kinetics, which is designed to sample reaction events. That is, if you have a number of possible reactions, this is a prescription for the next reaction that occurs and the time at which it occurs. In the language of the Gillespie method, the propensity function would be the probability $p_k$. See Eq. 10b and the procedure that follows in math.pitt.edu/~swigon/Homework/gillespie1.pdf $\endgroup$ – GnomeSort Jun 5 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.