9
$\begingroup$

Let $k<d$ two positive integers, and $\{p_i\}_{i=1}^d$ a series of probabilities, with $p_i \in (0,1)$ and $\sum_{i=1}^d p_i = k$.

We wish to sample exactly $k$ distinct indices $\mathcal{I}\triangleq\{i_1,...,i_k\}$ such that $$\forall i\in \{1,\dots,d\}: P(i\in \mathcal{I}) = p_i \,.$$ How can we do so efficiently?

I solved a few easy cases (e.g., $k=2$, $d=4$), but I couldn't find a general computationally feasible method.

$\endgroup$

1 Answer 1

9
$\begingroup$

A direct approach is to use a wheel of fortune with $k$ equidistant pointers at angular positions $$\zeta_j~=~\alpha+ 2\,\pi\,j/k,~~~~j=0,\ldots,k-1,$$ with equi-distributed $\alpha \in (0,2\pi/k]$ and randomly permuted sectors of angular size $2\,\pi\, p_i/k$. Choose the indices whose sectors are indicated by the pointers. To show that for each $i=1,\ldots,d$ one has $P(i\in I) ~=~p_i$ consider the angular sector belonging to the index $i$ can at most be reached by two pointers $\zeta_j, \,\zeta_{j+1 \mod k}$, and that the sector can be divided in two sub-sectors each of which can be reached by only one of these pointers, but not at the same time. From that it is easily calculated that the combined probability to reach the angular sector belonging to the index $i$ by choosing an equi-distributed $\alpha \in (0,2\pi/k]$ indeed is $p_i$.

A slightly simplified version (using intervals instead of sectors) is implemented in the following Mathematica function that expects a list indprob of $d$ pairs $(i, p_i)$ and the value of $k$ fulfilling the condition $\sum p_i =k$. It then returns a list of $k$ indices having the corresponding probabilities.

RandomIndices[indprob_, k_] := Module[{ransam = RandomSample[indprob], acc, pointer, ind, 
alpha = RandomReal[{0, 1}]}, 
acc = Accumulate[ransam[[;; , 2]]]; 
pointer = Range[k] + alpha - 1;
ind = With[{p = #}, FirstPosition[acc, _?((# > p) &)]] & /@ pointer;
ransam[[Flatten[ind]]][[;; , 1]] // Sort]

To destroy potentially remaining correlations between the indices one could further subdivide the sectors of each index into several pieces and permute the then much larger number of sectors.

$\endgroup$
3
  • $\begingroup$ Thanks, I'm verifying this now. Two questions: 1. How can a sector be reached by two pointers? The pointer difference is $2\pi / k $ is larger than the sector size $2\pi p_i / k$ since $p_i<1$. So this implies that a sector can be reached only by one pointer, right? 2. If we further sub-divide the sectors and permute (like you suggested in the last sentence), can't this cause several identical indices to be sampled? $\endgroup$ Commented Jul 23 at 9:08
  • 1
    $\begingroup$ If one of your sectors overlaps with e.g. $\zeta_1$ it is hit by pointer 1 for $\alpha$ close to 0, and by pointer 0 for $\alpha$ close to $2\pi/k$. You are right with the further subdivision, this requires care to make sure that the sectors for each index e.g. lie within a sector of size $2\pi/k$. To answer your question this is not necessary, but I guess you also like to have no unnecessary correlations between the indices. If you do not care, you can leave out the permutation. Maybe permutation is already sufficient to destroy correlation, but I could not prove this so far. $\endgroup$ Commented Jul 23 at 9:24
  • 1
    $\begingroup$ 1. Ah, I see, I thought you meant 2 hits for the same $\alpha$ sample. 2. Yeah, I'm not sure yet myself if I need this. Thanks again, I'll approve this nice answer. $\endgroup$ Commented Jul 23 at 10:31

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .