Sampling from specific random distribution on sets

Question

I have a random distribution on sets in mind, that has three parameters: $n, w, k$. The goal is to sample sets of $k$ integers from $[0, n)$ (without replacement) such that the elements within each set fit in a subrange of length $w$. That is, an outcome set $S$ must have properties:

1. $S \subset \mathbb{N_0} \; \wedge\; |S| = k$
2. $0\leq \min(S) \leq \max(S) < n$
3. $\max(S) - \min(S) < w$

You can assume that $k \leq w/2 < w \ll n$.

Now there are many possible distributions possible over these sets. But I'm interested in those that have as property

$$\forall x:P[x \in S] = \frac{k}{n}\;,$$

that is each integer in $[0, n)$ has an equal chance of being in a set when sampled (or as close as possible). Beyond the above requirements, it'd be ideal if the distribution is an maximum entropy one, but this isn't as important, and something close would be fine too. As a minimum bar I do think every valid set should have a non-zero chance of occurring.

Is there a practical way of sampling from a random distribution that matches the above requirements?

I've tried various methods, rejection sampling, first picking the smallest/largest elements, but so far everything has been really biased. The only method that works that I can think of is explicitly listing all valid sets $S_i$, assigning a probability variable $p_i$ to each, and solving the linear system $$\sum_i p_i = 1 \quad\bigwedge\quad \forall_x:\frac{k}{n} - \delta \leq \sum_{x \in S_i} p_i \leq \frac{k}{n} + \delta,$$ minimizing $\delta$ first, $\epsilon $ second where $\epsilon = \max_i p_i - \min_i p_i$. However this is very much a 'brute force' approach, and is not feasible for larger $n, k, w$.

Answer 1

I have a suggestion for an algorithm that might get close to what you want:

1. Pick $m$ from the distribution $\mu$ (which I will describe below).

2. Choose $S$ uniformly at random from all sets containing exactly $k-1$ integers in $[m+1,\dots,\min(m+w,n))$.

3. Add $m$ to $S$. Output $S$.

So, how do we pick the distribution $\mu$ to make this algorithm pretty good? I'll describe how, using linear programming.

Let $\mu_i$ denote the probability of $i$ according to the distribution $\mu$. Then, we have

$$\begin{align*} \Pr[x \in S] &= \sum_m \Pr[x \in S \mid m] \Pr[m]\\ &= \mu_x + \sum_{m=x-w+1}^{x-1} {k-1 \over \min(w-1,n-m-1)} \cdot \mu_m. \end{align*}$$ Notice that this is linear in the $\mu_i$'s. Now, let's treat the $\mu_i$'s as variables to be solved for. Introduce the linear inequalities

$$\left| \Pr[x \in S] - {k \over n}\right| \le t,$$

where $t$ is another variable. Also add the equality $\mu_1 + \dots + \mu_{n-k}=1$ and $\mu_{n-k+1}=\cdots=\mu_{n-1}=0$ and the inequalities $\mu_i > 0$ for all $i$. We'll try to minimize $t$ subject to all of these inequalities. These are linear inequalities in the variables $\mu_0,\dots,\mu_{n-w},t$, so you can find a solution in polynomial time using linear programming. That gives you a distribution $\mu$ to use in the algorithm above.

This comes close to meeting your requirements. It ensures the set $S$ output by the algorithm meets your conditions 1, 2, and 3. It gets close to having $\Pr[x \in S] = k/n$; this is not exactly true, but it will be approximately true; and you can measure how large the error is when you solve the linear program, by inspecting the value of $t$. Also, due to the inequality $\mu_i > 0$, every possible set has a non-zero probability of occurring. Of course, this is only a heuristic, and it probably is not the maximum-entropy distribution, so you'll have to decide whether it is good enough for your needs.

You mention that you are dealing with a very large value of $n$. Given that, an additional heuristic would be to assume that $\mu_i$ is constant for all $i \in [3w,n-4w)$. Intuitively, the endpoints are the problem cases, and there is a loose symmetry among the middle values that makes it feel reasonable that they'll end up having about the same probability. So, we can enforce the constraint $\mu_{3w} = \mu_{3w+1} = \cdots = \mu_{n-4w-1}$, by replacing $\mu_i$ with $\mu_{3w}$ for all $i \in [3w,n-4w)$. Now there are only about $6w$ variables, namely, $\mu_1,\dots,\mu_{3w},\mu_{n-4w},\dots,\mu_{n-w}$ and $t$. Note that it is easy to express $\Pr[x \in S]$ as a linear sum of the $\mu_i$ variables; this sum can be computed in $O(w)$ time. So, we end up with a linear program in about $6w$ variables, so you can apply a LP solver, and the resulting system will have $O(w)$ variables and inequalities instead of $O(n)$ variables and inequalities, so the LP solver should be a lot faster.