# Probability distribution of the number of infected people in a room

Suppose we have a room of $$N$$ people. The number of people who have the common cold is $$k$$, which is equal to $$1$$ at the start. Now, $$h$$ handshakes occur completely randomly. If an infected person shakes hands with a uninfected person, the uninfected person gets infected and can then go on to infect other people if they shake hands with them. Otherwise, nothing happens. What is the probability distribution of $$k$$ after all the handshakes occur?

My attempts to solve this:

The probability of an additional person being infected with the cold in 1 handshake would be represented by $$\frac{2k(N-k)}{N(N-1)}$$ , but aside from this, I don't know where to go next.

I also looked at hypergeometric distributions and compound probability distributions. Would either be related?

• We can see it as follows, though maybe not an efficient way to solve it. The probability you wrote is the one for $k$ to increase by 1 after one handshake. We have $1-{}$ that probability for $k$ to stay the same. For all $k$ this defines a transition matrix, with elements on the main diagonal and the one below. After $h$ handshakes, this matrix is multiplied by itself $h$ times. The initial probability for $k$ is concentrated on $k=1$. The probability for $k$ after $h$ handshakes is then given by the $h$th power of the transition matrix times the initial (singular) probability distribution. Aug 27, 2019 at 14:23
• When you say "completely randomly" do you mean that two people are selected at random (without replacement, so you get a pair of distinct people) from the total of $N$, they shake hands, then you put them both back into the pool of $N$ before selecting again (so that you're selecting $h$ times from the $N\choose 2$ pairs, now with replacement)? Aug 28, 2019 at 4:50
• Yeah, that's what I meant. Aug 28, 2019 at 5:04

Let $$\mathbf{K} = \{ K_h | h \in \mathbb{N}_{0+} \}$$ denote the stochastic time-series showing the number of infected people after each handshake, and let $$K_0 = 1$$ at the start of the series. This is a Markov chain that falls within the category of discrete "pure birth" processes". A single random handshake gives the transition probabilities:

$$p_{k,k+r} \equiv \mathbb{P}( K_{h+1} = k+r | K_h = k ) = \begin{cases} 1-\frac{2k(N-k)}{N(N-1)} & & \text{if } r=0, \\[6pt] \frac{2k(N-k)}{N(N-1)} & & \text{if } r=1, \\[8pt] 0 & & \text{otherwise}. \\[6pt] \end{cases}$$

Thus, the transition matrix for the chain is:

$$\mathbf{P} \equiv \begin{bmatrix} 1-\frac{2}{N} & \frac{2}{N} & \cdots & 0 & 0 & 0 \\ 0 & \frac{4(N-2)}{N(N-1)} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & \cdots & 1-\frac{4(N-2)}{N(N-1)} & \frac{4(N-2)}{N(N-1)} & 0 \\ 0 & 0 & \cdots & 0 & 1-\frac{2}{N} & \frac{2}{N} \\ 0 & 0 & \cdots & 0 & 0 & 1 \\ \end{bmatrix}.$$

After $$h$$ random handshakes, the probability that $$k$$ people are infected is:

$$\mathbb{P}(K_{h} = k) = [\mathbf{P}^{h}]_{1,k}.$$

You can compute this probability by programming the transition probability matrix into an appropriate piece of mathematical software (e.g., R) and then obtaining the first row of the appropriate power of the matrix. If you would like to try to get a closed-form expression for the probability, I would recommend deriving the eigen-decomposition or Jordan decomposition of the matrix to see if this simplifies the problem.

Computing the probability vector: It is quite simple to program this Markov chain in R. In the code below I create a function to compute the vector of probabilities (or log-probabilities) for arbitrary input values for the number of people and the number of handshakes.

#Load required library
library(expm);

#Create a function to compute the probability vector
COMPUTE_PROBS <- function(N, h, log.p = FALSE) {

#Define the transition probability matrix
P <- matrix(0, nrow = N, ncol = N);
for (k in 1:N)     { P[k,k]   <- 1 - 2*k*(N-k)/(N*(N-1)); }
for (k in 1:(N-1)) { P[k,k+1] <- 1 - P[k,k]; }

#Compute probability vector
PPP <- expm::'%^%'(P,h);
if (log.p) { PPP <- log(PPP); }
PPP[1, ]; }


We can use this function to compute the vector of probailities (or log-probabilities) for arbitrary values of N and h. Here is an example using some chosen values for the parameters.

#Compute an example of this probability vector
N <- 40;
h <- 80;
PROBS <- COMPUTE_PROBS(N,h);

#Plot the probability mass function
library(ggplot2);
DATA  <- data.frame(Infected = 1:N, Probability = PROBS);
THEME <- theme(plot.title    = element_text(hjust = 0.5, size = 14, face = 'bold'),
plot.subtitle = element_text(hjust = 0.5, face = 'bold'));
ggplot(aes(x = Infected, y = Probability), data = DATA) +
geom_bar(stat = 'identity', fill = 'blue') + THEME +
ggtitle('PMF of Number of Infected People') +
labs(subtitle = paste0('(', N, ' people and ', h, ' handshakes)')) +
xlab('Number of Infected People');


Monte-Carlo simulation: We can confirm that the above result is correct by comparing the theoretical probabilities to Monte-Carlo simulations of the process. To do this, we can program a simulation function in R. (Hat tip to user2974951 for suggesting this approach, and writing the initial code.) In the code below I create a function to simulate outcomes of the chain and take empirical estimates of the vector of probabilities for arbitrary input values for the number of people and the number of handshakes.

#Create a function to simulate the chain
SIMULATE_CHAIN <- function(N, h, times = 10^5) {

#Set the simulation vector
SIM <- rep(0, times);

#Run simulations
for (s in 1:times) {

#Compute initial vector of infected people
INFECTED <- c(1, rep(0, N-1));

#Implement random handshakes
for (i in 1:h) {
H <- sample(1:N, size = 2, replace = FALSE);
if (INFECTED[H[1]] == 0 & INFECTED[H[2]] == 1) { INFECTED[H[1]] <- 1 }
if (INFECTED[H[1]] == 1 & INFECTED[H[2]] == 0) { INFECTED[H[2]] <- 1 } }
SIM[s] <- sum(INFECTED); }

SIM; }


Using this function we can simulate the Markov chain and take empirical estimates of the probabilities of each outcome. The plot confirms the same shape we obtained in our theoretical analysis, which confirms that the calculations are correct.

#Simulate the chain
set.seed(1)
SIMS <- SIMULATE_CHAIN(N,h);

#Estimate the probability vector
PROBS_EST <- rep(0,N);
for (i in 1:N) { PROBS_EST[i] <- sum(SIMS == i)/length(SIMS); }

#Plot the probability mass function
DATA  <- data.frame(Infected = 1:N, Probability = PROBS_EST);
THEME <- theme(plot.title    = element_text(hjust = 0.5, size = 14, face = 'bold'),
plot.subtitle = element_text(hjust = 0.5, face = 'bold'));
ggplot(aes(x = Infected, y = Probability), data = DATA) +
geom_bar(stat = 'identity', fill = 'red') + THEME +
ggtitle('Monte-Carlo estimate of PMF of Number of Infected People') +
labs(subtitle = paste0('(', N, ' people and ', h, ' handshakes)')) +
xlab('Number of Infected People') + ylab('Estimated Probability');


Here is an alternative to Ben's answer using simulations in R, using his parameters.

Edit: fixed the bug.

N=40 #number of people
h=80 #handshakes
k=1 #number of infected people at the start
n=1e5 #number of simulations

result=rep(NA,n)
for (r in 1:n) {
initial=rep(0,N) #N healthy people
initial[1:k]=1 #k infected
for (t in 1:h) {
random2=sample(1:N,2) #two random people
if (initial[random2[1]]==1 | initial[random2[2]]==1) {
initial[random2[1]]=initial[random2[2]]=1 #now both infected
}
}
result[r]=sum(initial)
}


which looks like this

• This is a great attempt and a useful addition to the discussion (+1), but unfortunately it (presently) has a bug in the code. As a general piece of coding advice, try to write your code out more clearly by using variable names that are meaningful to the reader, and adding comments describing what each part of the code is doing. I have added some simulation calculations to my own answer, using the same method, but with clearer coding of the function. Have a look at my code above to get some tips on structure/naming. Great work.
– Ben
Aug 28, 2019 at 22:43