# A fast algorithm for a probabilistic counting problem without replacement?

Consider the integers $$\{1,\dots, N\}$$ for some positive integer $$N$$. Let us suppose that for each $$\{1, \dots, N\}$$ there is an associated probability $$p_1, \dots, p_N$$. We also define an integer threshold $$1 \leq n < N$$.

We sample independently and repeatedly without replacement from $$\{1,\dots, N\}$$. For each $$i \in \{1,\dots, N\}$$ we sample the integer $$i$$ with probability $$p_i$$. As we are sampling without replacement these probabilities will have to be appropriately scaled after each new integer is sampled.

Is there a fast algorithm to compute the expected number of distinct integers less than or equal to the threshold $$n$$ in a sample of size $$x$$?

The expected value we want to compute is $$\sum_{i=1}^n q_i$$ where $$q_i = \sum_{y = 1}^x \sum_{\substack{\{j_k \neq i\text{ distinct}\}\\k=1,\dots,y-1}}\left(\prod_{k=1}^{y-1}\frac{p_{j_k}}{1 - \left(\sum_{\ell = 1}^{k-1} p_{j_{\ell}}\right)}\right)\frac{p_i}{1 - \left(\sum_{\ell = 1}^{y-1} p_{j_{\ell}}\right)}.$$

However, this is infeasible to compute for anything but the smallest problem instances.

Is there an efficient algorithm for this problem?

• Do you mean probability proportional to $p_i$ rather than equal to $p_i$? Even if initially we have $\sum_i p_i=1$ after the first sampling the sum of the remaining $p_i$ will not be 1 anymore. Apr 17 at 19:12
• @MaxAlekseyev Yes. After every integer is sampled the probability of the remaining integers will of course have to be scaled to make them add to one. Apr 17 at 19:33

I'm not sure why you ask for "distinct integers" when sampling without replacement guarantees distinctness.

Let $$q_i=1-p_i$$. The ordinary generating function $$F(u,y) = \prod_{i=1}^n (q_i+p_i uy) \prod_{i=n+1}^N (q_i+p_iy)$$ counts subsets with $$u$$ monitoring the choices in $$[1..n]$$ and $$y$$ monitoring all choices.

To extract the expectation of the number of choices in $$[1..n]$$, differentiate with respect to $$u$$. If your sample size is $$k$$ (sorry but I can't use $$x$$ for an integer parameter when playing with generating functions), extract the coefficient of $$y^k$$.

So $$\text{answer} = \frac{[y^k] F_u(1,y)}{[y^k] F(1,y)},$$ where $$[y^k]$$ means extract the coefficient of $$y^k$$.

Now to computation. Start with the denominator $$F(1,y) = \prod_{i=1}^N (q_i+p_i y).$$ I don't know of any super-clever way to extract the coefficient of $$y^k$$ but you can just start with the polynomial $$1$$ and multiply by one factor at a time discarding any powers greater than $$k$$. That will have complexity $$O(kN)$$. Probably there are options with lower theoretical complexity but the constant here is so small that it will be hard to beat for practical sizes.

The numerator is $$F_u(1,y) = \sum_{i=1}^n p_i y\prod_{j\ne i} (q_j+p_jy),$$ which you can compute in similar manner.

ADDED. I'll explain how to get a good estimate for large problems, assuming that $$k$$ is not very small. I'll consider only the denominator. The function $$F(1,y)$$ is the pgf of the sum of $$N$$ independent Bernoulli random variables. The sum satisfies a local central limit theorem (LCLT) but that won't help immediately unless $$k$$ is close to the mean of the sum. Now introduce a parameter $$\alpha\gt 0$$ and replace $$q_i+p_iy$$ by $$(q_i + \alpha p_iy)/(q_i+\alpha p_i)$$ for each $$i$$. Solve for $$\alpha$$: $$\sum_{i=1}^N \frac{\alpha p_i}{q_i+\alpha p_i}=k,$$ which should be fairly easy numerically as there is no cancellation and the left side is monotonically increasing in $$\alpha$$.

With this value of $$\alpha$$, the mean of the sum of the tilted Bernoulli random variables with pgf $$(q_i + \alpha p_iy)/(q_i+\alpha p_i)$$ is exactly $$k$$. Now you can apply a LCLT such as Theorem 6.3 in this paper to get quite a good estimate.

If you want an arbitrarily accurate value, you can apply a fourier inversion theorem to the characteristic function and evaluate it numerically to the precision you desire. I don't know if this will be easy, but for sure it will be advantageous to tilt the distribution first as I did above.

• I'm not sure why you ask for "distinct integers" ... That surely comes from a previous (and equivalent) problem statement, where one samples with replacements (leaving the probabilities unchanged) and stops after getting $n$ different values. math.stackexchange.com/questions/4096509/… May 11 at 1:45
• Yes probably, but "without replacement" and "distinct" are the same, which could have indicated a mistake in the description. May 11 at 6:09
• @leonbloy Yes it was exactly that. May 11 at 15:09
• @BrendanMcKay thank you for this answer. I would be very interested in any methods with lower theoretical complexity. In particular that might be feasible when $N \approx 10^8$ and $k \approx 10^5$. May 11 at 15:11
• Thank you again. May 22 at 12:05