A probability question that I failed to answer in a job interview

Question

There are $N$ different sized targets, numbered $1, 2,\dots, N$. A blindfolded shooter shoots towards random directions. Because target sizes are different, the hit probabilities of each target are different, say $p_1, p_2,\dots ,p_N$. A bullet hole is left on the target each time a target is hit. The shooter does not know whether a shot hit a target or not.

Note that a shoot could be "void", i.e., $p_1+p_2+\cdots+p_N\le1$. One shot at most hits one target.

The shooter keeps shooting until $X$ out of the $N$, $X<N$, targets have bullet holes. ( each of the $X$ targets is hit at least once). Then the shooter is told to stop.

The question is: at the end of the game, what is the probability that target $k$ has NOT been hit, $1\leq k\leq N$.

Answer 1

We may as well suppose the game continues until all targets have been hit (which will happen eventually if all $p_j > 0$; we may as well remove any targets that have $p_j = 0$).

For each subset $S$ of $\{1,\ldots, N\}$, let $p_S = \sum_{s \in S} p_s$ be the probability that a shot hits a member of $S$, and let $a_{i,t}(S)$ be the probability that target $i$ is one of the first $t$ targets in set $S$ to be hit. You want $1 - a_{i,X}(\{1,\ldots,N\})$. Of course $a_{i,t}(S) = 0$ if $i \notin S$, and we also take it to be $0$ if $t = 0$. Otherwise, conditioning on the first target in $S$ to be hit, $$ a_{i,t}(S) = \dfrac{p_i}{p_S} + \sum_{j \in S \backslash \{i\}} \dfrac{p_j}{p_S} a_{i,t-1}(S \backslash \{j\}) $$ Now I claim that $$ a_{i,t}(S) = \sum_{T \subseteq S \backslash \{i\}} c(t,|T|,|S|) \frac{p_i}{p_{T \cup \{i\}}} $$ for some constants $c(t,m,n)$, $0 \le m \le n-1$. I will prove this by induction on $t$. In the case $t=1$ we have $a_{i,1}(S) = p_i/p_S$, so $c(1,m,n) = 1$ if $m = n-1$, $0$ otherwise.

If $t >1$, we have (with $|S|=n$): $$ \eqalign{a_{i,t}(S) &= \dfrac{p_i}{p_S} + \sum_{j \in S \backslash \{i\}} \dfrac{p_j}{p_S} \sum_{T \subseteq S \backslash \{i,j\}} c(t-1, |T|,n-1) \dfrac{p_i}{p_{T \cup \{i\}}}\cr &= \dfrac{p_i}{p_S} + \sum_{m=0}^{n-2} c(t-1,m,n-1) \sum_{T \subseteq S \backslash \{i\}: |T| = m} \sum_{j \in S \backslash (T \cup \{i\})} \dfrac{p_j p_i}{p_S p_{T \cup \{i\}}}\cr &= \dfrac{p_i}{p_S} + \sum_{m=0}^{n-2} c(t-1,m,n-1) \sum_{T \subseteq S \backslash \{i\}: |T| = m} \dfrac{p_{S \backslash (T \cup \{i\})} p_i}{p_S p_{T \cup \{i\}}}\cr &= \dfrac{p_i}{p_S} + \sum_{m=0}^{n-2} c(t-1,m,n-1) \sum_{T \subseteq S \backslash \{i\}: |T| = m} \dfrac{(p_S - p_{T \cup \{i\}}) p_i}{p_S p_{T \cup \{i\}}}\cr &= \dfrac{p_i}{p_S} + \sum_{m=0}^{n-2} c(t-1,m,n-1) \sum_{T \subseteq S \backslash \{i\}: |T| = m} \left(\dfrac{p_i}{p_{T \cup \{i\}}} - \dfrac{p_i}{p_S}\right)\cr &= \sum_{T \subseteq S \backslash \{i\}} c(t,|T|,n) \dfrac{p_i}{p_{T \cup \{i\}}} }$$ where $c(t, m,n) = c(t-1, m,n-1)$ if $m < n-1$ while $$c(t, n-1, n) = 1 - \sum_{m=0}^{n-2} {n-1 \choose m} c(t-1,m,n-1)$$

Hmm: it looks like $$ c(t, m, n) = \cases{1 & for $t=n,m=0$\cr (-1)^{n+m+t} {m-1 \choose {t+m-n}} & $ n-m \le t \le n$\cr 0 & otherwise\cr}$$ There ought to be an inclusion-exclusion proof for this.

Answer 2

At the risk of shooting myself in the foot, here's how I'd approach it:

The key principle is that we can ignore any action that doesn't produce a tangible result. In particular, this includes missing a target, so we can start by normalizing the probabilities: write $\tilde{p_i} =\frac{p_i}{\sum_{1\leq k\leq N}p_k}$ so that we can work with values that have $\sum\tilde{p_i}=1$.

But from here, it should be clear that the process is just drawing without replacement - after our trial, we'll have some set of $X$ items, and the probability that those items are $i_0, i_1, \ldots, i_X$ is $\tilde{p}_{i_0}\cdot\tilde{p}_{i_1}\cdots\tilde{p}_{i_X}$. So the probability that target $k$ is hit is just $\dfrac{\sum_{\{S:\left|S\right|=X \wedge k\in S\}} \mathbb{P}_S}{\sum_{\{S: \left|S\right|=X\}} \mathbb{P}_S}$, where the lower sum is over all $X$-element subsets $S$ of $1\ldots N$ and the upper sum is over all those subsets with $k\in S$, and $\mathbb{P}_S=\prod_{i\in S}\tilde{p}_i$ (and of course, the probability that $k$ isn't hit is just the complement of this value).

EDIT: as pointed out by Nate Eldredge in a comment below, the formula can't simplify to one independent of the other probabilities; let the below serve as a cautionary tale about making assumptions!

This formula itself should offer further simplification, I'm reasonably sure, but that will take a little bit of chewing. (I strongly suspect that everything else is just a red herring and the final probability is a simple expression in $\tilde{p}_k$ and $\tilde{q}_k=1-\tilde{p}_k$, but I need to work through the details and have my ducks in a row before I'm certain of that.)

Answer 3

I don't think it's easy. Say we have p1 = p2 = 0.49, and p3 .. p22 = 0.01. For the first two objects, the chance of being the first object is 0.49, the chance of being first or second is a bit above 0,98, and then it gets close to 1 very, very quickly.

For the other objects, the chance of being the first object is 0.01, the chance of being one of the first two is about 0.03, and to be among the first 2 + n the chance is about n/20.

If the probabilities are less extreme then I think it gets quite difficult, for example 20 items with 0.03 ≤ pk ≤ 0.07.

Answer 4

Guys please criticize the following answer to my own question.

Let variable $L$ be the number of shots by the end of the game. Then, clearly, the probability of the $k$-th target survives is $(1-p_k)^L$.

Then the real question is how to determine $L$ or the $E(L)$, the expectation of $L$.

Let $q_i=1-p_i$ be the probability of missing a target for one single shot.

With $L$ shots, the probability that the $i$-th target has been hit at least once is $1-q_i^L$. Then given a set of targets $S$, where $k\notin S$ and $|S|=X$, the probability that all targets in $S$ have been hit at least once, let's call it $S$ completed, is that \begin{equation} P_{\mbox{S completed}} = \prod_{i\in S} (1-q_i^L) \end{equation}

Note that there are $C_{N-1}^{X}$ target set $S$ satisfying $k\notin S$ and $|S|=X$, forming a set $\tilde{S} = \{S \mid k\notin S \mbox{ and }|S|=X \}$. As said, $|\tilde{S}|=C_{N-1}^{X}$.

Because the game is over when $L$ shots have been made, then the probability that at least one $S\in\tilde{S}$ completed equals to 1 (I'm not sure if I'm correct by saying this). Or \begin{equation} \sum_{S\in\tilde{S}} P_{\mbox{S completed}} = 1 \end{equation}

With all of this, we have \begin{equation} \sum_{S\in\tilde{S}} \prod_{i\in S} (1-(1-p_i)^L) = 1 \end{equation}

Knowing $p_i$, $1\leq i \leq N$, $L$ can be obtained using computation tools like Matlab or R. (I haven't verified the solution with the tools yet)

Answer 5

I think the right way to think about this is to realize that the game is over if and only if you've hit $X$ unique targets, so the relevant game states are all subsets of $[1, N]$ of size $X$ - the ones that contain $k$ are states where you've hit target k before the game is over.

Looking at pure combinations, the number of states where you've hit k is:

$$n_{k} = {X - 1 \choose N - 1}$$

Because you're selecting the first one in the set to be $k$ and selecting $X - 1$ from the remaining $N - 1$ combinations. The number of combinations that do not include $k$ are:

$$n_{!k} = {X \choose N - 1}$$

So in the limit where all targets are equally likely, the answer would be:

$$\mathbf{p}_{!k} = \frac{n_{!k}}{n_{k} + n_{!k}}$$

(Note that the denominator is the total number of end game states.)

When the probabilities of hitting targets is different, I think you need to use the same approach, but using the approach from this question to calculate the sum of the product of all the relevant subsets.

Using the notation from the Wikipedia article on elementary symmetric polynomials, and calling $\mathbf{T}_{!k} = \{p_i\}$ where $i \in [1, N]$ and $i \neq k$:

$$P_{k} = p_{k} \cdot e_{X-1}(\mathbf{T}_{!k})$$ $$P_{!k} = e_{X}(\mathbf{T}_{!k})$$

And the answer is: $$\mathbf{p}_{!k} = \frac{P_{!k}}{P_{k} + P_{!k}}$$

Edit To verify that this is the correct answer, I wrote a python script that runs an empirical test by actually randomly playing the game a number of times and counting the number of games in which $k$ was hit:

from __future__ import division

from copy import copy
from itertools import combinations
from operator import mul
from functools import reduce, partial

import numpy as np

# Empirical approach
def play_game(t_list, k, X):
    if k >= len(t_list):
        raise ValueError("k cannot be greater than N")

    if X >= len(t_list):
        raise ValueError("X must be <= N")

    targ_probs = np.array(t_list) / sum(np.array(t_list))
    targ_inds = np.arange(0, len(targ_probs))
    targ_set = set()

    while True:
        c_target = np.random.choice(targ_inds, p=targ_probs)

        if c_target == k:
            return False

        targ_set.add(c_target)

        if len(targ_set) >= X:
            return True

def empirical_probability(n_games, t_list, k, X):
    game_outcomes = [play_game(t_list, k, X) for ii in range(n_games)]
    return sum(game_outcomes) / n_games

# Analytical approach
def elementary_symmetric_polynomial(m, t_set):
    return sum(map(partial(reduce, mul), combinations(t_set, m)))

def analytical_probability(t_list, k, X):
    t_k = t_list[k]
    t_set = np.array(t_list)
    t_set = t_set[[ii for ii in range(0, len(t_list)) if ii != k]]

    # All subsets of size X not containing k
    p_nk = elementary_symmetric_polynomial(X, t_set)

    # All subsets of size X - 1 not containing k
    p_k = t_k * elementary_symmetric_polynomial(X-1, t_set)

    return p_nk / (p_k + p_nk)

if __name__ == "__main__":
    n_games = 5000

    # Randomly generate a game set
    N = 6
    t_list = np.random.random((N,))
    t_list /= sum(t_list)

    empirical = [['X = '] + [str(X) for X in range(2, N-1)]]
    analytical = copy(empirical)
    for k in range(0, len(t_list)):
        ana_answers = ['k = {}'.format(k)]
        emp_answers = ['k = {}'.format(k)]
        for X in range(2, N-1):
            ana_answers.append(round(analytical_probability(t_list, k, X), 3))
            emp_answers.append(round(empirical_probability(n_games, t_list, k, X), 3))

        empirical.append(emp_answers)
        analytical.append(ana_answers)

    print("Empirical answers:")
    for row in empirical:
        print('\t'.join(str(x) for x in row))

    print("Analytical answers:")
    for row in analytical:
        print('\t'.join(str(x) for x in row))

Choosing $N = 6$ and running the game 5000 times, here is one run of results:

p = [ 0.3305, 0.0350, 0.04995, 0.0831, 0.1503, 0.3509]
Empirical answers:
X =     2   3   4
k = 0   0.373   0.153   0.045
k = 1   0.916   0.833   0.721
k = 2   0.892   0.779   0.591
k = 3   0.813   0.642   0.407
k = 4   0.671   0.432   0.199
k = 5   0.348   0.141   0.041
Analytical answers:
X =     2   3   4
k = 0   0.397   0.218   0.108
k = 1   0.908   0.811   0.649
k = 2   0.871   0.741   0.542
k = 3   0.792   0.608   0.374
k = 4   0.652   0.416   0.225
k = 5   0.38    0.207   0.102

I think these differences are well within the noise of the empirical measurement. The algorithm for the elementary symmetric polynomial calculation came from this question.

Answer 6

The question is obviously too difficult to be answered analytically, at least in a job interview. What it's not so difficult is to design an algorithm to evaluate the answer numerically. (Perhaps that was the idea in the interview?)

An equivalent formulation: we have $N$ balls with weights $p_1, p_2 \cdots p_N$. We extract the balls without replacement, with probabilities equal to their normalized weights. We want $P_{j,t}$, the probability that, after $t$ extrations, ball $j$ has not been yet extracted.

In particular, $P_{j,0}=1$, $P_{j,N}=0$,$P_{j,1}=1-p_j$,

To get the idea, suppose $t=3$, $j=N$. Then to get $P_{j,t}$ we need to sum over all the probabilities of the length $3$ extraction starting subsequences that don't include $j=N$. For example, (assuming normalized initial weighs) the subsequence $(1,2,3\cdots )$ would have probability

$$ p_1 \frac{p_2}{1-p_1} \frac{p_3}{1-p_1-p_2} $$ and so on.

Then $${P_{j,t}= \sum_{\sigma \in {\mathbb P}_{N,t,j} } \prod_{i=1}^t \frac{p_{\sigma_i}}{\left(1- \sum_{k=1}^{i-1} p_{\sigma_k} \right)}}$$

where ${\mathbb P}_{N,t,j} $ denotes all the permutations of $t$ elements, taken from the set $\{1, 2 \cdots N \}$ with element $j$ excluded.

The sum over such permutations can be concisely computed by recursion. For example, in Java:

public class Balls {

    // not necessarily normalized
    final static double prob[] = { 0.3, 0.2, 0.2, 0.15, 0.1, 0.05 };
    final static int N = prob.length;

    /**
     * Computes probabiity that element j (j=0..N-1, same index as in prob array)
     * does not appear in the first pos positions

     * computeProb(element, 0) = 1
     * computeProb(element, 1) = 1-p_j 
     * computeProb(element, N) = 0 
     */
    public static double computeProb(int j, int pos) {
        normalize(prob);
        double[] p0 = removeAndNormalize(prob, j, false);
        return allPermutations(p0, pos);
    }

    static double allPermutations(double p[], int depth) {
        double ac = 0;
        if( depth <= 0 ) return 1;
        for( int i = 0; i < p.length; i++ ) 
            ac += p[i] * allPermutations(removeAndNormalize(p, i, true), depth - 1);
        return ac;
    }

    // normalizes probabities in place
    static void normalize(double[] list) {
        double ac = 0;
        for( double px : list )         ac += px;
        for( int i = 0; i < list.length; i++ )          list[i] /= ac;
    }

    // removes element of array, and optionally normalizes 
    // original (untouched) must be already normalized 
    static double[] removeAndNormalize(double p[], int index, boolean normalize) {
        double[] q = new double[p.length - 1];
        for( int i = 0, j = 0; i < q.length; i++, j++ ) {
            if( i == index ) j++;
            q[i] = normalize ? p[j] / (1 - p[index]) : p[j];
        }
        return q;
    }

For example, here's some results for $N=6$

prob=[0.3, 0.2, 0.2, 0.15, 0.1, 0.05]
t=0  1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 
t=1  0.70000 0.80000 0.80000 0.85000 0.90000 0.95000 
t=2  0.44794 0.59624 0.59624 0.68615 0.78423 0.88919 
t=3  0.24803 0.39457 0.39457 0.50573 0.64561 0.81148 
t=4  0.10545 0.20779 0.20779 0.30642 0.46950 0.70305 
t=5  0.02477 0.06330 0.06330 0.11122 0.21732 0.52009 

prob=[60.0, 5.0, 4.0, 3.0, 2.0, 0.1]
t=0  1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 
t=1  0.19028 0.93252 0.94602 0.95951 0.97301 0.99865 
t=2  0.02777 0.63673 0.70856 0.78082 0.85349 0.99264 
t=3  0.00277 0.34721 0.43412 0.54870 0.68413 0.98307 
t=4  0.00014 0.12570 0.18124 0.27503 0.45035 0.96754 
t=5  0.00000 0.00568 0.00942 0.01701 0.03533 0.93256

In IdeOne, so you can play with it.