# Sampling from specific random distribution on sets

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$$.

• Just to check that I understand correctly: does solution "sample every possible subset with equal probability" has problems only for $x$ near the end of the interval? I.e. it has the same probability for each $x \in [w, n-w)$ and some other probabilities for other $x$'s, and therefore it's not acceptable?
– user114966
Oct 24, 2020 at 17:38
• @Dmitry Correct.
– orlp
Oct 24, 2020 at 18:06
• In this case, it seems to me that you can do the following: to fix the probability at $x=0$ (the same for other $x$), consider a mixture of distribution produced by $D_1=$"sample every possible subset with equal probability" and $D_2=$"sample subsets so that 0 has probability $p_1$ and any other number has probability $p_2 < p_1$". How to build the latter? As an example, consider case $n \mod k = 1$. We assigne to all subsets of $[0, k+1)$ of size $k$, containing $0$ (i.e. $\{0,2,3,...,k\}$, $\{0,1,3,...,k\}$, ..., $\{0,1,2,...,k-1\}$), weight $1$ (weight is a non-normalized probability).
– user114966
Oct 24, 2020 at 18:10
• Note that $0$ has weight $k$ and all $x \in [1,k+1)$ have weight $k-1$. Now, we assign to all sets $[k+1, 2k+1)$, $[2k+1, 3k+1)$ weight $k-1$. Therefore, each element $>0$ has weight $k-1$, while $0$ has weight $k$. It remains to find with which probability sample from $D_1$ and with which - from $D_2$, which is trivial if you know initial probabilities for each number.
– user114966
Oct 24, 2020 at 18:10
• @Dmitry Sorry, I don't think I understand your proposed scheme. On another note, when I said that you're correct I forgot one caveat: I'm not 100% sure I know how to sample from $D_1$ in an efficient way (that is, a concrete algorithm) either. There are some methods I conjecture are pretty close, but I have no convincing proof (or even a suspicion) that they are actually correct and not biased. On another note, I'd prefer it if you didn't re-use parameters $n,k$ or $w$ for other purposes because it makes the notation very confusing.
– orlp
Oct 24, 2020 at 18:31

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.

• Your approach almost works, but there's still some issues. One is trivial (you forgot to bound $\sum_i \mu_i = 1$), but a bigger issue is that you didn't explain how to handle the right boundary. If $m > n - w$ then your step $2$ might select integers outside $[0, n)$. Finally, $n$ can be quite large (> $10^8$), so ideally I'd avoid a linear programming solution.
– orlp
Oct 25, 2020 at 3:40
• @orlp, good points! See edited answer for some refinements.
– D.W.
Oct 25, 2020 at 6:18
• There's still an issue on the right boundary, and that is when $m > n - k$. It's not possible at all to choose $k-1$ other elements without replacement in step 2 in that scenario. Your assumption that $\mu_i$ is (nearly) constant in the middle is also false I believe (at least in this approach), you end up with a sort of 'ripple' effect. In order to compensate for the extra weight when $0 \leq i < w$ you find that $w \leq i < 2w$ has reduced weight, but that in turn gives extra weight to $2w \leq i < 3w$, etc. At least that was the case when I tried this for a similar problem I've had.
– orlp
Oct 25, 2020 at 6:29
• @orlp, Regarding the right boundary, I see that I was unclear. See edited answer. OK, I see, I didn't anticipate the ripple effect. Can you come up with a better one? Perhaps plot the value of $\mu_x$ as a function of $x$, for some tractable value of $n$, without assuming it is constant in the middle, and show that in the question and see if it suggests a better heuristic?
– D.W.
Oct 25, 2020 at 16:41