A fast algorithm for a probabilistic counting problem without replacement?

Question

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?

Answer 1

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.