Algorithm to take samples from binomial data

Question

I am faced with the following problem : I have binomial data, ie. data of the form trials/successes, with a list of features.

What I need is to sample without replacement from that data, ie. take a fraction $p$ (in my case 50%)of the examples. Since the data is aggregated, this is tricky.

So far my algorithm is the following : Denote by $h_i,c_i$ the number of trials/successes on each line, and $h'_i,c'_i$ the sampled trials/successes. I compute the total number of trials $H$ and then randomly generate $pH$ distinct integers between $1$ and $H$. I also compute the cumulative sum $(s_i = \sum_{t\leq i}h_t)_{i<n}$ of successes. Then for each sampled integer $k$, I find the first index $i$ such that $k \geq s_i$, and increment $h'_i$. If $k-s_i \leq c_i$ then I increment $c'_i$.

This algorithm works, but the problem is that when generating the random integers, I must store the whole list of integers, which means allocating space proportional to $H$, which in my case is incredibly large, and causes heap to overflow (I use JAVA with a 2Gb heap).

Any ideas to avoid storing that huge vector ? Or another simpler algorithm would be nice (I'm sure this one is largely suboptimal).

Edit : I want to add that the overall complexity of the algorithm should not be more than $O(log(n)H)$ where $n$ is the number of lines.

Answer 1

If I have understood correctly, the crux of the matter is to generate a random subset of the integers $1, 2, \ldots, H$ of predetermined size $k = pH$ without having to store the entire subset. This suggests we try to generate it as an ascending sequence $0\lt i_1\lt i_2\lt \cdots\lt i_k\le H$ by producing the $i_j$ in order, thereby reducing RAM requirements to a small constant value. To that end, note that for any $x \in H$ the chance that $i_1 \gt x$ is the proportion of $k$- subsets of $\{1, 2, \ldots, H\}$ that are actually subsets of $\{x+1, x+2, \ldots, H\}$ and thereby equals

$$\binom{H-x}{k} / \binom{H}{k} = \frac{(H-x)(H-x-1)\cdots(H-x-k+1)}{H(H-1)\cdots(H-k+1)}.$$

The probability mass function for this distribution consequently is

$$p_{H,k}(x)= -\binom{H-x}{k} / \binom{H}{k} + \binom{H-(x-1)}{k} / \binom{H}{k} = \binom{H-x}{k-1} / \binom{H}{k},$$

$\ 1 \le x \le H-k+1.$

Therefore if you can sample $i_1$ efficiently from this distribution, you can proceed recursively to sample $i_2 - i_1$ from the distribution determined by $H-i_1$ and $k-1$ and so on. Until the very end of this process, the ratios $k(j)/H(j)$ = $(k-j)/(H-(i_1+\cdots+i_j))$ will stay close to $p$.

When $p\gg 0,$ the mass of this distribution is concentrated on small $x$. (After all, the expected gaps between $i_j$ and $i_{j+1}$ are around $1/p$, so $x$ will only rarely be larger than a small multiple of $1/p$). For large $H$, this distribution is closely approximated by a geometric distribution with ratio $1-p$ ($ = 1/2$ in this case) and the approximation is excellent for small $x$. This indicates that a simple rejection-sampling procedure will work well, creating a $O(k)$ algorithm.