Find the expected number of hours all 512 people suffer

Question

In a crowded room with 512 people, there are 9 different types of blood antigens. Each person has a blood type corresponding to a distinct subset of the antigens; the remaining of the antigens are foreign to them.

Matt the Mosquito flies around the room, picks a subset of the people uniformly at random, and bites the chosen people in a random order. After biting a person, Matt stores a bit of any antigens that person had. A person bitten while Matt had k blood antigens foreign to him or her will suffer for k hours. What is the expected total suffering of all 512 people in hours?

Let $n=9$ be the number of different antigens.

I think the linearity of expectation could be useful. Note that each person is bitten with probability 1/2. Consider a person P and an antigen a foreign to him/her. Assuming P has been bitten, it would be useful to find the probability P will suffer due to a. There are $2^{n-1}$ people with antigen a. But I'm not sure how to find the probability P will suffer due to a. Though with this given probability, it should be possible to use linearity over pairs (P,a) where P is a person and a is an antigen foreign to that person.

Answer 1

I agree: linearity of expectation is a good approach here.

Notation setup

I'm going to treat the problem size as a parameter: Let's say there are $N$ types of antigens and $2^N$ people in the room. (Your statement is the special case $N=9$.)

Let $S$ be the total suffering across all people, $S_p$ be the total suffering for person $p$, and $S_{p,a}$ be the amount of suffering of person $p$ from antigen $a$. Let $P_{p,a}$ be the probability that person $p$ ends up suffering from antigen $a$. Thus $S_{p,a}$ will always be $1$ (with probability $P_{p,a}$) or $0$ otherwise.

Reducing using linearity of expectation

You're asking for $E[S]$, and we can write $$\begin{align} E[S] &= E\left[ \sum_{p} S_p \right] \\ &= E\left[ \sum_{p,a} S_{p,a} \right] \\ &= \sum_{p,a} E[S_{p,a}] \\ &= \sum_{p,a} P_{p,a} \end{align}$$ where I've used linearity to pull the sum outside the $E$.

Solving the reduced problem

Opening comment: Biting a uniform random subset of people is the same as flipping an independent fair coin to decide whether to bite each person. I find the second interpretation much more useful in terms of intuition and I recommend readers think about it this way.

With the L.o.E. reduction in mind, the problem is reduced to finding $P_{p,a}$. If person $p$ already has antigen $a$, or if $p$ never gets bitten, then they certainly won't suffer from $a$. So let's focus on the case where person $p$ does not have antigen $a$ and does get bitten. In that case, there are $2^{N-1}$ other carriers who do have antigen $a$. Some number of those carriers will be bitten; the probability of exactly $k$ carriers being bitten is $\binom{2^{N-1}}{k} / 2^{\left(2^{N-1}\right)}$.

Given that exactly $k$ carriers are bitten, the probability of person $p$ catching disease $a$ is exactly $\frac{k}{k+1}$, since if we ignore everyone except those $k$ carriers and $p$, they'll escape the disease only if they are the first one bitten out of these $k+1$ people.

The probability that person $p$ gets bitten is $\frac 1 2$, so we can write (assuming $p$ does not start out with antigen $a$) $$\begin{align} P_{p,a} &= \frac 1 2 \left[ \frac{1}{2^{\left(2^{N-1}\right)}} \sum_{k=0}^{2^{N-1}} \binom{2^{N-1}}{k} \frac{k}{k+1} \right] \\ &= \frac 1 2 \left(\frac{1}{2^{\left(2^{N-1}\right)}} \right) \left(\frac{1 - 2^{\left(2^{N-1}\right)} + 2^{\left(2^{N-1}\right)} 2^{N-1}}{2^{N-1}+1}\right) \\ &= \frac 1 2 \left[ \frac{1}{2^{\left(2^{N-1}\right)}\left(2^{N-1} + 1\right)} + \frac{2^{N-1}-1}{2^{N-1}+1} \right] \end{align}$$

I'm omitting proof of the step where I replace the $\sum_k$ with a closed form. It can be shown by combinatorics or by induction, but easiest is by WolframAlpha. I used this query and got these results:

Anyway, now we can put the subanswers together to answer the full problem. For each $a$, exactly $2^{N-1}$ people start out with the disease (hence have no chance of catching it) and the other $2^{N-1}$ people start without the disease (hence their chance of catching it is given by the messy formula above). So our final answer is:

$$\begin{align} E[S] &= \sum_a \sum_p P_{p,a} \\ &= \sum_a 2^{N-1} \frac 1 2 \left[ \frac{1}{2^{\left(2^{N-1}\right)}\left(2^{N-1} + 1\right)} + \frac{2^{N-1}-1}{2^{N-1}+1} \right] \\ &= \boxed{N 2^{N-2} \left[ \frac{1}{2^{\left(2^{N-1}\right)}\left(2^{N-1} + 1\right)} + \frac{2^{N-1}-1}{2^{N-1}+1} \right]} \end{align}$$

Intuition and double-checking

Intuitively, unless $N$ is tiny, if person $p$ starts out without antigen $a$ and gets bitten then I would pretty much expect them to catch $a$. Like if $N = 100$, there are just so so many infected people that had a chance to get bitten first; it seems pretty unlikely that person $p$ will avoid infection. This is borne out by the boxed formula above: all the stuff inside $\left[ \cdots \right]$ will be pretty close to $1$ for large $N$ (or really for any $N$ bigger than like $2$ or $3$), so the answer will be pretty close to $N 2^{N-2}$ which is what we'd get if every bitten person catches every disease they don't start with.

We can also double-check the exact formula using a Python simulation. I used this code:

# use the formula we derived
def exact(N):
    return N*2**(N - 2)*(1/(2**(2**(N - 1))*(2**(N - 1) + 1)) + (2**(N - 1) - 1)/(2**(N - 1) + 1))

import random
# simulate a bunch of times to get an approximate answer
def sim(N, NUM_TRIALS):
    total_suffering = 0
    for trial_num in range(NUM_TRIALS):
        current_suffering = 0

        # decide who gets bitten
        bitten = []
        for p in range(2**N):
            if random.randint(0,1) == 1: bitten.append(p)
        # random bite order
        random.shuffle(bitten)
        
        diseases = 0  # track which diseases the fly is carrying
        for p in bitten:
            diseases |= p  # fly catches all of p's diseases
            current_suffering += diseases.bit_count() - p.bit_count()  # suffer 1 hour per disease that the fly has but p didn't have before
        total_suffering += current_suffering
    return total_suffering / NUM_TRIALS

for n in range(0, 10):
    ex = exact(n)
    si = sim(n, 10**5)
    print(n, ex, si, sep='\t')

and I got these printouts:

0   0.0 0.0
1   0.125   0.12479
2   0.8333333333333333  0.83008
3   3.675   3.67416
4   12.45138888888889   12.44043
5   35.29415355009191   35.31057
6   90.1818181824955    90.21302
7   217.1076923076923   217.13148
8   504.06201550387595  504.04763
9   1143.035019455253   1142.75689

That's enough accuracy to make me feel mostly convinced that the formula is right. I also tried running with more trials for some of the smaller $N$ values and confirmed that the simulation results got a little closer to the exact answers.

Additional confirmation: Full brute force for small $N$

A skeptical commenter suggested using a full brute force solver to confirm the exact answer for small $N$. Here's some Python code:

import itertools
import math
import fractions
def count_suffering(bitten):
    suffering = 0
    diseases = 0
    for p in bitten:
        diseases |= p  # fly catches all of p's diseases
        suffering += diseases.bit_count() - p.bit_count()  # suffer 1 hour per disease that the fly has but p didn't have before
    return suffering

def bruteforce(N):
    total_suffering = fractions.Fraction(0)
    for num_bitten in range(2**N+1):
        for bitten in itertools.combinations(range(2**N), num_bitten):
            current_suffering = fractions.Fraction(0)
            for bitten_ordered in itertools.permutations(bitten):
                current_suffering += count_suffering(bitten_ordered)
            total_suffering += current_suffering / math.factorial(num_bitten)
    bitten_options = 2**(2**N)  # total number of subsets we could have picked
    expected_suffering = total_suffering / bitten_options
    return expected_suffering

for n in range(4): print(bruteforce(n))

and the output:

0
1/8
5/6
147/40

I've confirmed that these are the same as my formula's output for $N=0, 1, 2, 3$. If I really cared I could code this in C++ and then the $N=4$ case would also finish in a reasonable amount of time. I'm not bothering because I honestly just feel sufficiently convinced that the formula is correct already, and I think Python is a better language to post here since more people can read it.

Bonus notes

The exact answer to the specific version you asked ($N = 9$) is: $$\frac{265742844799640668497095410594938748523254614807645094470555155298160632523653129}{232488804171798923623888618337756189986643641086481444985473430390888080605184}$$

Viewing the answer as a function $S(N)$, it turns out that $\tilde S(N) = N 2^{N-2} - N$ is an extremely close approximation. The first several values of $S(N) - \tilde S(N)$ are:

0   0.
1   0.625
2   0.833333
3   0.675
4   0.451389
5   0.294154
6   0.181818
7   0.107692
8   0.0620155
9   0.0350195
10  0.0194932

So for most practical purposes you'd be fine using that much simpler formula instead.