A dynamic program will make short work of this.
Suppose we administer all questions to the students and then randomly select a subset $\mathcal{I}$ of $k=10$ out of all $n=100$ questions. Let's define a random variable $X_i$ to compare the two students on question $i:$ set it to $1$ if student A is correct and student B not, $-1$ if student B is correct and student A not, and $0$ otherwise. The total
$$X_\mathcal{I} = \sum_{i\in\mathcal{I}} X_i$$
is the difference in scores for the questions in $\mathcal I.$ We wish to compute $\Pr(X_\mathcal{I} \gt 0).$ This probability is taken over the joint distribution of $\mathcal I$ and the $X_i.$
The distribution function of $X_i$ is readily calculated under the assumption the students respond independently:
$$\eqalign{
\Pr(X_i=1) &= P_{ai}(1-P_{bi}) \\
\Pr(X_i=-1) &= P_{bi}(1-P_{ai}) \\
\Pr(X_i=0) &= 1 - \Pr(X_i=1) - \Pr(X_i=0).
}$$
As a shorthand, let us call these probabilities $a_i,$ $b_i,$ and $d_i,$ respectively. Write
$$f_i(x) = a_i x + b_i x^{-1} + d_i.$$
This polynomial is a probability generating function for $X_i.$
Consider the rational function
$$\psi_n(x,t) = \prod_{i=1}^n \left(1 + t f_i(x)\right).$$
(Actually, $x^n\psi_n(x,t)$ is a polynomial: it's a pretty simple rational function.)
When $\psi_n$ is expanded as a polynomial in $t$, the coefficient of $t^k$ consists of the sum of all possible products of $k$ distinct $f_i(x).$ This will be a rational function with nonzero coefficients only for powers of $x$ from $x^{-k}$ through $x^k.$ Because $\mathcal{I}$ is selected uniformly at random, the coefficients of these powers of $x,$ when normalized to sum to unity, give the probability generating function for the difference in scores. The powers correspond to the size of $\mathcal{I}.$
The point of this analysis is that we may compute $\psi(x,t)$ easily and with reasonable efficiency: simply multiply the $n$ polynomials sequentially. Doing this requires retaining the coefficients of $1, t, \ldots, t^k$ in $\psi_j(x,t)$ for $j=0, 1, \ldots, n.$ (we may of course ignore all higher powers of $t$ that appear in any of these partial products). Accordingly, all the necessary information carried by $\psi_j(x,t)$ can be represented by a $2k+1\times n+1$ matrix, with rows indexed by the powers of $x$ (from $-k$ through $k$) and columns indexed by $0$ through $k$.
Each step of the computation requires work proportional to the size of this matrix, scaling as $O(k^2).$ Accounting for the number of steps, this is a $O(k^2n)$-time, $O(kn)$-space algorithm. That makes it quite fast for small $k.$ I have run it in R
(not known for excessive speed) for $k$ up to $100$ and $n$ up to $10^5,$ where it takes nine seconds (on a single core). In the setting of the question with $n=100$ and $k=10,$ the computation takes $0.03$ seconds.
Here is an example where the $P_{ai}$ are uniform random values between $0$ and $1$ and the $P_{bi}$ are their squares (which are always less than the $P_{ai}$, thereby strongly favoring student A). I simulated 100,000 examinations, as summarized by this histogram of the net scores:
The blue bars indicate those results in which student A got a better score than B. The red dots are the result of the dynamic program. They agree beautifully with the simulation ($\chi^2$ test, $p=51\%$). Summing all the positive probabilities gives the answer in this case, $0.7526\ldots.$
Note that this calculation yields more than asked for: it produces the entire probability distribution of the difference in scores for all exams of $k$ or fewer randomly selected questions.
For those who wish a working implementation to use or port, here is the R
code that produced the simulation (stored in the vector Simulation
) and executed the dynamic program (with results in the array P
). The repeat
block at the end is there only to aggregate all unusually rare outcomes so that the $\chi^2$ test becomes obviously reliable. (In most situations this doesn't matter, but it keeps the sofware from complaining.)
n <- 100
k <- 10
p <- runif(n) # Student A's chances of answering correctly
q <- p^2 # Student B's chances of answering correctly
#
# Compute the full distribution.
#
system.time({
P <- matrix(0, 2*k+1, k+1) # Indexing from (-k,0) to (k,k)
rownames(P) <- (-k):k
colnames(P) <- 0:k
P[k+1, 1] <- 1
for (i in 1:n) {
a <- p[i] * (1 - q[i])
b <- q[i] * (1 - p[i])
d <- (1 - a - b)
P[, 1:k+1] <- P[, 1:k+1] +
a * rbind(0, P[-(2*k+1), 1:k]) +
b * rbind(P[-1, 1:k], 0) +
d * P[, 1:k]
}
P <- apply(P, 2, function(x) x / sum(x))
})
#
# Simulation to check.
#
n.sim <- 1e5
set.seed(17)
system.time(
Simulation <- replicate(n.sim, {
i <- sample.int(n, k)
sum(sign((runif(k) <= p[i]) - (runif(k) <= q[i]))) # Difference in scores, A-B
})
)
#
# Test the calculation.
#
counts <- tabulate(Simulation+k+1, nbins=2*k+1)
n <- sum(counts)
k.min <- 5
repeat {
probs <- P[, k+1]
i <- probs * n.sim >= k.min
z <- sum(probs[!i])
if (z * n >= 5) break
if (k.min * (2*k+1) >= n) break
k.min <- ceiling(k.min * 3/2)
}
probs <- c(z, probs[i])
counts <- c(sum(counts[!i]), counts[i])
chisq.test(counts, p=probs)
#
# The answer.
#
sum(P[(1:k) + k+1, k+1]) # Chance that A-B is positive