6
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – pglpm
    Aug 27, 2019 at 14:23
  • $\begingroup$ 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)? $\endgroup$
    – Glen_b
    Aug 28, 2019 at 4:50
  • $\begingroup$ Yeah, that's what I meant. $\endgroup$
    – cash999
    Aug 28, 2019 at 5:04

2 Answers 2

7
$\begingroup$

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');

enter image description here


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');

enter image description here

$\endgroup$
0
2
$\begingroup$

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

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ben
    Aug 28, 2019 at 22:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.