Assignment 3

8.3. Assignment 3#

Same instructions as past assignments.

This assignment involves a small amount of work to modify the code provided below. You can download a “jupyter-notebook” using the download link at the top. Alternately you can use the ‘rocket’ icon at the top to open it in Binder or in your browser. Be careful opening in your browser as it will quit if you leave it idle too long.

  1. Consider an infectious disease which has recently been introduced into a very large population from wildlife. It is believed that the offspring distribution can be well-approximated by a Negative-Binomial distribution of unknown parameters. There have been five known independent introductions into the population. So far all outbreaks have died out, but there is concern that the reproduction number (the mean of the offspring distribution) may be larger than \(1\). We will use the parametrization of the Negative Binomial distribution from Example 1.11.

    (a) Find the probability of an outbreak of size \(L\) for a given \(p\) and \(r\) (\(r\) may not be an integer). Your expression should be given in terms of the Gamma function \(\Gamma(x)\) — replace all factorials \(M!\) with \(\Gamma(M+1)\). It will be helpful to observe that if \(\mu_{r,p}(x)\) is the PGF of the distribution, then \(\mu_{r,p}(x)^N\) is the PGF for another Negative Binomial distribution.

    Solution (a)

    In this model, each node in the Galton-Watson tree represents a distinct individual. So the number of infections is equal to the number of nodes in the tree. Following Size distribution of small outbreaks Section 4.5.1, the size of the outbreak is \(L\) with probability

    \[\begin{align*} \mathbb{P}[\text{outbreak of size $L$}] &= \frac{1}{L} [x^{L-1}] \mu_{r,p}(x)^L\\ &= \frac{1}{L} [x^{L-1}] \left(\left(\frac{p}{1-(1-p)x}\right)^r\right)^L\\ &= \frac{1}{L}[x^{L-1}] \mu_{rL,p}(x)\\ &= \frac{1}{L}[x^{L-1}] \frac{\Gamma(N+rL)}{N! \,\Gamma(rL)}p^{rL} (1-p)^N x^N\\ &= \frac{1}{L} \frac{\Gamma(L-1+rL)}{(L-1)! \,\Gamma(rL)}p^{rL}(1-p)^{L-1} \end{align*}\]

    (b) Find a reasonably simplified formula for the probability that \(M\) outbreaks have sizes \(L_1\), \(L_2\), \(\ldots\), \(L_M\) for given parameters \(p\) and \(r\) (expressed in terms of Gamma functions)

    Solution (b)

    As these are independent, it is just the product:

    \[\begin{align*} \mathbb{P}[L_1,L_2,\ldots, L_m] &= \prod_{i=1}^M \frac{1}{L_i} \frac{\Gamma(L_i-1+rL_i)}{(L_i-1)! \,\Gamma(rL_i)}p^{rL_i}(1-p)^{L_i-1}\\ &= \frac{1}{\prod_{i=1}^M L_i} \left(\prod_{i=1}^M \frac{\Gamma(L_i-1+rL_i)}{(L_i-1)! \, \Gamma(rL_i)}\right) p^{r\sum L_i}(1-p)^{-M + \sum L_i} \end{align*}\]

    The probability of observing data for a given set of parameters is called the likelihood of the parameters. The log-likelihood is the logarithm of the likelihood. The parameters that maximize the likelihood will also maximize the log-likelihood, but maximization is often easier with log-likelihood, and also much easier to use computationally. In Python, we can calculate \(\ln \Gamma(x)\) using

    from scipy.special import gammaln  #put this at the top

gammaln(x)

    (c) Find the log-likelihood of \(p\) and \(r\) given \(L_1, \ldots, L_M\).

    Solution (c)

    We just take logarithms, recalling that the logarithm of a product is the sum of the logarithms.

    We use \(\log\Gamma(x)\) for \(\ln(\Gamma(x))\)

    \[\begin{align*} \text{Log-likelihood}(r,p) &= + r \ln(p) \sum L_i + \ln(1-p)\left(-M+ \sum L_i\right)\\ &\quad + \sum_{i=1}^M {\Large(}-\ln L_i + \log\Gamma (L_i-1 + r L_i) - \ln ((L_i-1)!) - \log\Gamma(rL_i){\Large)} \end{align*}\]

    (d) Using \(A = \sum L_i\), find a function \(\hat{p}(r)\) that maximizes the likelihood given \(r\).

    Solution (d)

    Holding \(r\) fixed, we differentiate with respect to \(p\).

    \[\begin{align*} \frac{\partial}{\partial p} \text{Log-likelihood}(r,p) &= Ar \frac{1}{p} - \frac{A-M}{1-p}\\ \Rightarrow 0 &=Ar \frac{1}{\hat{p}} - \frac{A-M}{1-\hat{p}}\\ &= (1-\hat{p})Ar - (A-M)\hat{p}\\ \Rightarrow \hat{p} &= \frac{Ar}{A(r+1)-M} \end{align*}\]

    Note - for any \(r\) this is a maximum.

    (e) Write down a function in terms of only \(r\) to maximize to find \((p^*,r^*)\) that maximizes the likelihood. Remove any terms that do not depend on \(r\).

    Solution (e)

    We substitute this in for \(p\) in the log-likelihood and then differentiate with respect to \(r\):

    \[\begin{align*} \text{Log-likelihood}(r,\hat{p}(r)) &= + Ar \ln\left(\frac{Ar}{A(r+1)-M}\right) + (A-M)\ln\left(1-\frac{Ar}{A(r+1)-M}\right)\\ &\quad + \sum_{i=1}^M {\Large(}-\ln L_i + \log\Gamma (L_i-1 + r L_i) - \ln ((L_i-1)!) - \log\Gamma(rL_i){\Large)} \end{align*}\]

    So to maximize this, we can maximize the function

    \[\begin{align*} f(r) &= Ar \ln\left(\frac{Ar}{A(r+1)-M}\right) + (A-M)\ln\left(1-\frac{Ar}{A(r+1)-M}\right)\\ &\quad + \sum_{i=1}^M {\Large(}\log\Gamma (L_i-1 + r L_i)- \log\Gamma(rL_i){\Large)} \end{align*}\]

    (f) One way to estimate parameters is “Maximum Likelihood Estimation” (MLE). We simply find the parameters that have the highest likelihood.

    The following code simulates numsims\(=5\) outbreaks with a negative binomial offspring distribution having \(r=1\), \(p=0.6\). Then it uses maximum likelihood estimation to estimate \(r\) and \(p\).

    • Using your answers above, modify Function_to_maximize so that it performs the correct calculation. Also modify the expression for p_est. Run the code.

    • Play around with numsims to get an idea of how many simulations are needed for this to perform well at estimating the parameters. Also play around with \(p\) and \(r\) (the code will insist that the simulations be guaranteed to finish, so it will refuse to run if \(\mathcal{R}_0=\mu'(1)=(1-p)r/p>1\).).

    • Looking at the plot of \(p_k\), the offspring distribution, comment on whether \(p\) and \(r\) need to be known exactly to make good predictions.

    • Reflect on your observations. What happens if numsims is small? What if \(\mathcal{R}_0\) is close to 1?

    Solution (f)

    See code below.

    We need a large number of simulations to have good estimates for \(p\) and \(r\). However, well before \(p^*\) and \(r^*\) are a very close match the offspring distribution fits well. Thus the outcomes will be accurately predicted event if \(p\) and \(r\) are not well-identified.

    If \(\mathcal{R}_0\) is close to \(1\), we need higher certainty about its value to make good predictions for the future (the policy implications of being unsure if \(\mathcal{R}_0 =0.5\) or \(0.95\) are small compared to the implications of uncertainty about whether \(\mathcal{R}_0 = 0.95\) or \(1.05\)).

Hide code cell source

 
import numpy as np
from scipy.special import gammaln
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt

r = 0.8
p = 0.45
numsims = 5000

if    (1-p)*r/p >= 1:
    raise ValueError("Warning: mean must be less than 1 for extinction to be certain")


L=[]
for counter in range(numsims):
    X = 1
    cumulative = X
    while X>0:
        newX = 0
        for i in range(X):
            U = np.random.negative_binomial(r,p)
            cumulative += U
            newX += U
        X = newX
    L.append(cumulative)


print("observed sizes: ", L)
L = np.array(L)
A = L.sum()



#print("student-created function is needed here, using default")
def F_to_maximize(r):
    f = A*r*np.log(A*r/(A*(r+1)-numsims)) + (A-numsims)*np.log(1-A*r/(A*(r+1)-numsims))
    for Li in L:
        f += gammaln(Li - 1 + r*Li) - gammaln(r*Li)
    return f
def F_to_minimize(r):
    return -F_to_maximize(r)


res = minimize_scalar(F_to_minimize, bounds=(1e-8, 1e6), method="bounded")

r_est = res.x

#print("student-created function is needed here, using default")
p_est = A*r_est/(A*(r_est+1)-numsims)#\frac{Ar}{A(r+1)-M}

print("Estimated r:", r_est)
print("Estimated p:", p_est)

from scipy.stats import nbinom
plt.figure()
n = np.arange(0, 40)
plt.plot(n, nbinom.pmf(n, r_est, p_est ), label="Estimated: (r, p, R_0) = (%.2f, %.2f, %.6f)"%(r_est,p_est, (1-p_est)*r_est/p_est))
plt.plot(n, nbinom.pmf(n, r, p ), label="True: (r, p, R_0) = (%.2f, %.2f, %.6f)"%(r,p, (1-p)*r/p))
plt.legend()
plt.xlabel('offspring number')
plt.ylabel('probability mass function')
plt.show()

observed sizes:  [3, 4, 1, 1, 1, 1, 1, 2, 9, 1, 1, 1, 1, 1, 2, 3, 1, 1152, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 16, 1, 6, 2, 257, 70, 1, 14, 1, 11, 1, 1, 5, 1, 1, 1, 1, 1, 1, 1, 17, 2, 40, 1, 1, 1, 2, 1, 5172, 4, 1, 1, 2, 2, 4, 1, 1, 5, 1, 1, 1, 7, 176, 1, 1, 4, 1, 2, 1, 2, 1, 2, 1, 2, 1, 1, 1, 4, 1, 5, 2, 1, 1, 1, 1, 26, 3, 1, 1, 2, 1, 2, 24, 1, 312, 25, 3, 1, 1, 1, 6, 2, 1, 11, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 2, 23, 4, 1, 1, 8, 8, 1, 1, 2, 1, 11, 1, 3, 1, 12, 1, 3, 37, 10, 124, 1, 2, 2, 16, 1, 2, 1, 1, 1, 6, 1, 53, 146, 1, 10, 64, 1, 1, 1, 1, 1, 1, 1, 1, 38, 1, 1, 1, 1, 34, 4, 15, 25, 2, 5, 1, 108, 2, 1, 5, 2, 174, 6, 54, 9, 1, 1, 1, 1, 1, 35, 1, 1, 1, 1, 368, 1, 1, 2, 1, 1, 1, 3, 8, 1, 603, 1, 3, 1, 1, 3, 1, 5, 19, 1, 1, 3, 4, 1, 1, 3, 32, 2, 1, 1, 1, 2, 12, 1, 39, 41, 1, 1, 16, 3, 1, 8, 2, 1, 10, 1, 1, 58, 9, 14, 1, 2, 1, 9, 1, 62, 13, 1, 7, 2, 4, 5, 1, 1, 4, 1, 1, 4, 1, 285, 10, 10, 6, 5, 4, 1, 1, 1, 2, 1, 1, 2, 1, 14, 21, 2, 2, 1, 3, 2, 1, 5, 2, 1, 1, 2, 1, 1, 1, 1, 558, 1, 1, 1, 1, 62, 1, 1, 6, 32, 11, 34, 2, 9, 1, 1, 1, 1, 2, 1, 1, 12, 1, 74, 1, 2, 1, 3, 10, 1, 1, 1, 1, 1, 19, 1, 1, 1, 1, 2, 2, 1, 2, 10, 1, 184, 1, 3, 1, 1, 1, 1, 1, 1, 1, 2, 1, 48, 19, 2, 3, 1, 1, 5, 1, 1, 10, 1, 1, 1, 4, 3, 3, 8, 1, 1, 2, 3, 13, 1, 18, 1, 46, 6, 2, 1, 242, 3, 2, 1, 3, 1, 18, 1, 1, 2, 1, 1, 37, 3, 1, 3, 1, 1, 2, 1, 2, 442, 1, 30, 17, 3, 1, 7, 2, 1, 1, 1, 1, 23, 1, 1, 4, 5, 1, 1, 1, 1, 1, 1, 5, 2, 134, 2, 1, 1, 1, 2, 2, 1, 2, 2, 1, 31, 1, 4, 51, 2, 1, 1, 2, 1, 1, 1, 2, 1, 4192, 1, 5, 1, 4, 1, 1, 1, 49, 1, 5, 3, 406, 3, 1, 1, 4, 1, 1, 2, 8, 1864, 1, 1, 1, 1, 4, 1, 1, 13, 1, 17, 1, 2, 1, 1, 345, 1, 1, 2, 1, 5, 59, 22, 1, 3, 2, 1, 1, 2, 3, 1, 1, 2, 6, 1, 1, 1, 87, 9, 1, 1, 1, 12, 24, 1, 89, 1, 2, 1, 2, 15, 1, 16, 1, 1, 3, 2, 25, 23, 1, 1, 2, 3, 2, 2, 13, 1, 1, 2, 5, 1, 1, 16, 2, 2, 1, 4, 32, 1, 2, 18, 2, 29, 1, 1, 1, 79, 2, 1, 2, 7, 7, 1, 1, 111, 1, 1, 5, 1, 1, 51, 1, 2, 18, 3, 1, 2, 1, 2, 4, 171, 32, 1, 1, 3, 3, 1, 1, 1, 1, 29, 2, 6, 4, 1, 1, 8, 3, 1, 1, 1, 36, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 4, 1, 5, 1, 4, 1348, 1, 1, 1, 1, 1, 1, 1, 4, 2, 71, 167, 4, 1, 1, 3, 13, 1, 1, 39, 1, 1, 4, 5, 1, 53, 1, 10, 1, 1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 1, 2, 40, 2, 1, 1, 2, 4, 1, 12, 1, 3, 1, 3, 4, 4, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 13, 1, 2, 1, 6, 2, 6, 1, 29, 1, 1, 16, 6, 1, 53, 3, 1, 1, 26, 2, 20, 67, 40, 1, 1, 5, 2, 2, 1, 3, 1, 5, 23, 9, 1, 2, 15, 1, 309, 3, 302, 1, 1, 1, 2, 1, 2, 2, 3, 33, 20, 6, 1, 1, 8, 8, 1, 1, 17, 1, 2, 1, 13, 1, 40, 6, 1, 1, 1, 9, 1, 1, 3, 1, 1, 1, 1, 50, 3, 1, 1, 1, 6, 2, 4, 1, 2, 17, 26, 2, 9, 1, 55, 1, 2, 1, 1, 2, 4, 20, 1, 8, 3, 1, 2, 1, 847, 1, 8, 1, 14, 1, 103, 1, 1, 2, 1, 1, 130, 5, 1, 1, 10, 2, 8, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1129, 1, 1, 1, 1, 2, 3, 1, 1, 43, 1, 3, 1, 1, 1, 121, 37, 2370, 2, 1, 1, 3, 1, 1, 32, 13, 1, 97, 1, 2, 1, 3, 1, 1, 1, 1, 2, 1, 1, 4, 1, 5, 3, 1, 1, 1, 1, 1, 1, 177, 1, 1, 1, 43, 5, 2, 2, 2, 1, 1, 1, 1, 3, 1, 1, 2, 5, 2, 1, 1, 2, 21, 1, 109, 5, 2, 1, 3, 2, 13, 5, 1, 6, 2, 2, 1, 4, 2, 1, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 1, 6, 1, 1, 1, 5, 1, 1, 23, 863, 1, 2, 4, 1, 7, 1, 1, 1, 3, 1, 1, 20, 10, 1, 5, 4, 2, 1, 1, 1, 1, 4, 1, 1, 3, 1, 1, 3, 2, 1, 32, 1, 1, 1, 24, 1, 1, 1, 4, 16, 242, 2, 23, 1, 1, 1, 1, 271, 1, 2, 1, 1, 1, 1, 1, 2, 8, 1, 1, 1, 1, 1, 1, 1, 6, 2, 1, 1, 2, 1, 1, 3, 1, 118, 3, 1, 1021, 2, 22, 1, 2, 4, 1, 3, 1, 1, 14, 1, 1, 1, 2, 2, 1, 12, 4, 1, 1, 1627, 10, 10, 1, 1, 15973, 1, 9, 2, 1, 1, 1, 1, 10, 2, 1, 1, 1, 29, 8, 3, 2, 19, 1, 1, 3, 5, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 2, 1, 1, 10, 1, 1, 1, 3, 1, 1, 9, 1, 1, 1, 1, 2, 2, 1, 1, 12, 1, 1, 1, 236, 2863, 2, 4, 2, 3, 11, 1, 4, 3, 23, 6, 1, 1, 1, 10, 1, 3, 1, 2, 2, 4, 5, 1, 8, 1711, 1, 1, 42, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1287, 3, 1, 1, 89, 1, 1, 1, 1, 1, 148, 1, 2, 553, 8, 10, 164, 2, 1, 1, 1, 1, 19, 1, 1, 1, 1, 3, 14, 1, 5, 1, 19, 1, 14, 1, 1, 149, 2, 1, 8, 3, 1, 1, 1, 521, 4, 1, 4, 1, 1, 2, 6, 3, 2, 22, 1, 1, 2, 9, 1, 247, 2, 1, 1, 2, 1, 5, 1, 1, 1, 45, 1, 119, 10, 4, 1, 3, 1, 22, 2, 40, 2, 1, 1, 1, 1, 373, 1, 2, 1, 1, 1, 1, 4, 1, 1, 1, 1, 4, 6, 3, 2, 5, 2, 6, 1, 2, 1, 4, 15, 2, 22, 1, 1, 1, 3, 1, 1, 1, 1, 1, 7, 2, 1, 1, 1, 1, 5, 1, 1, 443, 1, 1, 60, 45, 1, 1, 1, 1, 52, 1, 2, 1, 3, 1, 1, 1, 2, 1, 1, 1, 1, 2, 19, 3, 1, 1, 1, 2, 4, 1, 1, 9, 8, 1, 1, 1, 2, 14, 4, 7, 9, 6, 1, 1, 1, 1, 12, 3, 1, 34, 1, 2, 1, 1, 72, 1, 2, 326, 1, 2, 82, 1, 1, 1, 1, 1, 2, 1, 15, 8, 16, 1, 9, 2, 1, 141, 1, 2, 7, 4, 2, 3, 1, 2, 1, 13, 2, 1, 293, 6, 2, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 275, 1, 1, 1, 1, 1, 3, 16, 1, 1, 1, 3, 1, 1, 2, 3478, 38, 1, 1, 1, 16, 1, 1, 7, 1, 1, 27, 1, 1, 1, 1, 35, 2, 1, 1, 2, 1, 65, 12, 1, 1, 2, 4, 4, 1, 1, 3, 2, 3, 2, 1, 2, 2, 120, 1, 3, 61, 1, 8, 1, 1, 16, 28, 1, 2, 10, 64, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 6, 1, 1, 1, 6, 42, 1, 4, 1, 1, 2, 1, 1, 60, 4, 1, 1, 2, 16, 167, 2, 3, 42, 8, 1, 1, 6, 2, 1, 1, 2, 1, 4, 7, 2, 2, 43, 1, 1, 1, 2, 1, 4, 8, 5, 1, 1, 1, 1, 14, 5, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, 1, 4, 41, 2, 1, 1, 2, 2, 1, 3, 6, 53, 1, 1, 268, 9, 13, 1, 1, 2, 49, 7, 1, 4, 1, 1, 7, 1, 2, 22, 1, 17, 7, 1, 1, 639, 68, 2, 1, 5, 1, 54, 1, 2, 52, 2, 1, 700, 3, 12, 1, 1, 1, 1, 1961, 1, 1, 3, 60, 2, 155, 1, 1, 1, 1, 3, 1, 1, 1, 113, 116, 218, 4, 5, 377, 1, 1, 3, 3, 1, 2, 1, 5, 1, 1, 2, 2, 1, 1, 2, 2, 17, 2, 2, 1, 1, 1, 38, 3, 15, 1, 3, 1, 1, 8, 13, 3, 1, 1, 1, 59, 1, 1, 4, 4, 1, 1, 1, 1, 6, 1, 42, 7, 1, 11, 4, 1, 2, 2, 1, 1, 16, 11, 1, 1, 2, 1, 1, 1, 4, 1, 14, 2, 9, 105, 1, 57, 29, 4, 3, 29, 1, 1, 1, 5, 1, 5, 11, 2, 1, 1, 2, 70, 3, 1, 4, 2, 1, 1, 1, 1, 9, 513, 1, 2, 27, 2, 2, 1, 1, 1, 44, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 7, 2, 53, 25, 1, 1, 2, 1, 363, 4, 26, 8, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 1, 83, 1, 1, 1, 4, 1, 67, 1, 1, 1, 1, 1, 6, 5, 174, 1, 3, 72, 2, 213, 4, 1, 4, 18, 1, 1, 2, 1, 1, 17, 1, 1, 11, 21, 1, 1, 1, 35, 12, 3325, 4, 2, 1, 1, 3, 1, 5, 2, 2, 1, 1, 8, 6, 1, 1, 1, 1, 3, 1, 1, 1, 1, 14, 22, 1, 1, 3, 2, 4, 1, 2, 4, 2, 1, 2, 257, 1, 1, 2, 1, 1, 36, 1, 1, 4, 1, 1753, 7, 12, 1, 2, 2, 1, 1, 1, 13, 2, 3, 119, 1, 7, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 5109, 1, 1174, 1, 3, 1, 9, 7, 1, 1, 2, 1, 1, 1, 5, 1, 1, 1, 1, 2, 1, 1, 1, 1, 14, 1, 1, 2, 1, 1, 3, 6, 3, 1, 1, 5, 1, 3, 1, 1, 22, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 4, 5, 1, 2, 1, 1, 1, 3, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1, 1, 1, 1, 45, 5, 1842, 1, 4, 12, 11, 15, 1, 43, 1, 1, 1, 1, 4, 1, 1, 1, 9, 7, 1, 1, 2, 3, 3, 5, 15, 1, 1, 1, 5, 1, 5, 1, 4, 6, 5, 1, 2247, 1, 1, 3, 3, 21, 18, 2, 1, 1, 1, 1, 66, 1, 15, 4, 1, 56, 4, 1, 1, 1, 2, 12, 1, 11, 1, 1, 3, 3, 5, 8, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 40, 1, 230, 1, 1, 1, 1, 5, 1610, 89, 1, 8, 1, 3, 1, 1, 134, 1, 1, 3, 1, 5, 2, 1, 2, 1, 1, 1, 21, 1, 14, 2, 41, 1, 1, 1, 3, 1, 6, 1, 5, 5, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 3, 42, 4, 1, 2, 3, 1, 2, 2, 1, 1, 18, 10, 1, 1, 1, 1, 6, 2, 3, 1, 9, 1, 1, 2, 154, 1, 1, 1, 1, 1, 1, 5, 628, 37, 2, 3, 25, 1, 1, 1, 5, 1, 1, 1, 1, 1, 1, 3, 1, 3, 1, 1, 474, 1, 5, 10, 3, 2, 6, 1, 1, 123, 4, 1, 1, 1, 5, 3, 3, 7, 2, 1, 1, 7, 1, 6, 4, 1, 22, 7, 1, 1, 1, 29, 5, 49, 1, 5, 1, 1, 1, 41, 103, 7, 1, 75, 1, 3, 1, 1, 1, 1, 1, 2, 77, 2, 1, 13, 3, 1, 4, 5, 1, 1, 1, 3, 1, 1, 2, 1, 1, 1, 4, 19, 2, 1, 10, 1, 6, 2, 3, 7, 1, 6, 2, 3, 1, 2, 1, 2, 2, 12, 4, 3, 1, 1, 38, 2, 1, 1, 1, 1, 22, 6, 2, 12, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 79, 1, 1, 3, 2, 1, 1, 1, 1, 1, 2, 37, 185, 1, 1, 1, 20, 1, 1, 52, 1, 2, 1, 1, 4, 2, 2, 6, 7, 309, 1, 103, 1, 2, 7, 8, 1, 1, 1, 1, 1, 1, 24, 3, 17, 92, 6, 6, 1, 1, 5, 78, 1, 1, 1, 1, 1, 20, 12, 6, 56, 2, 2, 3, 1, 4, 1, 24, 6, 2, 1, 1, 1, 1, 1, 1, 1, 9, 1, 1, 9987, 1, 23, 2, 1, 9, 1, 1, 3, 2, 1471, 2, 1, 1, 12, 1, 1, 2, 110, 1, 1, 10, 1, 4, 2, 1, 3, 4, 915, 1, 1, 5, 1, 6, 1, 6, 1, 105, 1, 13, 7394, 1, 22, 8, 1, 16, 31, 4, 6, 1034, 2, 1, 1, 4137, 1, 1, 5, 2, 549, 1, 2, 3, 1, 10, 1, 1, 456, 1, 1, 5, 1, 8, 9, 17, 2, 1, 3, 4, 3, 1, 1, 2, 1, 2, 2, 1, 278, 1, 2, 2, 1, 1, 56, 16, 5, 1, 1, 63, 6, 1, 1, 1, 1, 2, 79, 1, 1, 1, 2, 1, 24, 1, 3, 1, 2, 10, 2, 1, 3, 1, 1, 1, 14, 1, 1, 1, 36, 6, 1, 35, 1, 1, 21, 1, 2, 41, 6, 16, 1, 3, 1, 5, 63, 1, 1, 1, 1, 1, 1, 1, 2, 1, 26, 21, 1, 24, 1, 7, 1, 1, 35, 1, 10, 3, 1, 1, 1, 5, 1, 1, 3, 106, 3, 1, 2, 1, 5, 2, 37, 1, 1, 1, 1, 6, 4, 1, 3, 1, 1, 34, 5, 1, 3, 1, 1, 6, 1, 1, 7, 1, 1, 1, 1, 1, 1, 1, 1, 1417, 1, 1, 1, 3, 4, 1, 25, 1, 17, 1, 2, 5, 2, 2, 1, 1, 2, 1, 110, 1, 3, 2, 2, 1, 7, 1, 8, 1, 2, 1, 1, 7, 44, 1, 5, 29, 1, 6, 2, 35, 1, 1, 1, 5, 3, 2, 2, 1, 2, 1, 1, 1, 21, 13, 3, 30, 1, 1, 1, 1, 1, 13, 5, 1, 129, 3, 1, 1, 3, 1, 1, 1, 3, 1, 1, 3, 1, 2, 1, 1, 3, 1, 1541, 3, 1, 1, 1, 1, 1, 1, 10, 1, 30, 8, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 145, 229, 1, 3, 1, 2, 74, 1, 4, 1, 22, 1, 1, 19, 1, 1, 29, 156, 520, 22, 2, 15, 1, 9, 2, 692, 4, 7, 4, 16, 9, 1, 1, 3, 1, 2, 1, 1, 1, 6, 1, 1, 1, 1, 1, 1, 4886, 1, 1, 82, 1, 94, 4, 4, 1, 1, 11, 1, 1, 21, 1, 1, 1, 9, 1, 7, 12, 1, 2, 1, 137, 3, 2, 15, 1, 1, 2, 1, 8, 1, 3, 1, 7, 1, 2, 1, 20, 1, 1, 2, 1, 51, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 13, 4, 1, 2, 1, 1, 1, 1, 1, 13, 1, 1, 2, 67, 50, 9, 1, 5, 1, 5, 1, 76, 1, 69, 27, 1, 2, 36, 1, 1, 1, 1, 4, 3, 1, 1, 1, 1, 4, 1, 2, 1, 1, 5, 1, 4, 1, 1, 1, 47, 1, 3, 1, 6, 1, 1, 1, 52, 1, 2, 1, 2, 32, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 131, 1, 17, 3, 9, 1, 1, 5, 4, 1, 2, 2, 1, 1, 1119, 7, 3, 1, 3, 1, 2, 1, 1, 1, 6, 62, 1, 2, 2, 8, 1, 3, 1, 10, 159, 1, 1, 4, 6, 1, 2, 1, 333, 1, 1, 4, 1, 2, 71, 1, 1, 5, 1, 1, 1, 155, 1, 9, 1, 21, 1, 3, 6, 39, 1, 2, 2, 1, 14, 1, 1, 1, 8, 5, 1, 3, 37, 1, 1, 1, 1, 1, 2, 1, 1, 6, 2, 33, 2, 2, 2, 5, 1, 1, 1, 1, 57, 7, 1, 53, 1, 3, 1, 1, 3, 1, 1, 1, 1, 1, 2, 1, 1, 7, 170, 12, 91, 4, 10, 1, 1, 1, 1, 7, 8, 1, 14, 1, 2, 31, 1, 3, 7, 10, 1, 1, 3, 24, 1, 31, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 4, 9, 2, 1, 1, 19, 1, 1, 2, 1, 3, 1, 1, 17, 2, 1, 1, 1, 1, 3, 5, 1, 1, 1, 1, 39, 1, 1, 1, 8, 3, 1, 16, 19, 6, 1, 1, 1, 2, 2, 1, 1, 1, 56, 1, 6, 1, 7, 7, 1, 2, 1, 1, 4, 2, 1, 1, 2, 4, 2, 1, 44, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 46, 4, 7, 1, 2, 2, 1, 3, 12, 3, 3, 2, 1, 1, 2, 1, 3, 2, 1, 1, 1, 3, 1, 5, 1, 1, 1, 2, 1, 3, 35, 1, 6, 77, 3, 2, 1, 2, 1, 1, 3, 2, 1, 262, 1, 8, 1, 16, 5, 10, 1, 21, 2, 26, 2, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 2, 28, 3, 11, 1, 2, 5, 3, 1, 1, 13, 1, 1, 1, 1, 7, 1, 5, 35, 2, 1, 2, 1, 4, 3, 1, 1, 2, 1, 4, 1, 8, 1, 1, 127, 1, 4, 1, 96, 81, 1, 6, 2, 1, 1, 1, 5, 3, 1, 13, 28, 1, 1, 1, 1, 1, 3, 1, 1, 5, 1, 1, 1, 4, 1, 1, 3, 2, 2, 1, 1, 2, 1, 10, 1, 102, 1, 2, 2, 1, 1, 2, 9, 142, 604, 1, 4, 1, 1, 28, 1, 6, 1, 1, 24, 2, 1, 4, 2, 1, 1, 2, 7, 3, 1, 15, 1, 159, 1, 9, 2, 3, 1, 3, 9, 1, 1, 1, 8, 3, 7, 1, 6, 1, 1, 1, 11, 46, 1, 1, 1, 3, 2, 1, 1, 1, 1, 1, 4, 48, 1, 1, 766, 106, 1, 5, 1, 1, 16, 2, 1, 1, 16, 1, 1, 3, 5, 6, 2, 1, 6, 1, 1, 1046, 1, 1, 3, 1, 1, 33, 110, 1, 2, 2, 1, 9, 1, 3, 72, 11, 3, 1, 7, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 9, 1, 9, 1047, 2, 2, 2, 2, 1, 24, 4, 1, 1, 2, 5, 1, 270, 1, 1, 1, 1, 2, 11, 85, 1, 3, 1, 93, 1, 1, 6, 1, 1, 1, 2, 2, 2, 77, 1, 1, 3, 5, 3, 1, 1, 2, 3, 2, 6, 121, 59, 6, 45, 3, 1, 40, 1, 1, 4, 1, 1, 46, 2, 3, 1, 1, 3, 3, 1, 2, 3, 2, 1, 1, 1, 1, 1, 1, 3, 1, 3, 10, 4, 383, 1, 151, 1, 1, 1, 2, 1, 1, 1, 1, 3, 11, 11, 1, 2, 1, 1, 3, 5, 9, 2, 8, 191, 4, 1, 7, 33, 1, 1, 102, 1, 4, 1158, 89, 1, 5, 3, 1, 1, 73, 4, 1, 1, 1, 1, 1, 1, 2, 3, 62, 5, 1, 2, 1, 1, 1, 3, 5, 1, 1, 2, 1, 128, 1, 1, 1, 3, 1, 2, 1, 2, 6, 19, 3, 9, 1, 1, 1, 15, 1, 1, 1, 1, 1, 1, 4, 7, 59, 1, 2, 1, 1, 103, 1, 1, 7, 1, 745, 23, 5, 4, 1, 2, 1, 2, 1, 1, 4, 3, 2, 4, 79, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 3, 1, 1, 12, 1, 17, 1, 19, 1, 1, 1, 1, 2, 1, 1, 7, 6, 1, 1, 1, 1, 15, 1, 1, 1, 1, 1, 1, 102, 1, 93, 11, 1, 1, 32, 1, 1, 1, 5, 10, 1, 1, 77, 1, 3, 2, 92, 1, 1, 1, 2, 192, 1, 1, 10, 196, 4, 17, 141, 1, 19, 1, 1, 30, 1, 1, 1, 1, 3, 1, 10, 1, 3, 2718, 3, 12, 63, 5, 2, 1, 13, 1, 66, 1, 1, 1, 1, 1, 2, 20, 1, 2, 5, 3, 1, 1, 1, 2, 1, 1, 3, 1, 4, 1, 2, 2, 5, 8, 1, 10, 1, 6, 2, 1, 1, 3, 1, 11, 3, 1, 1184, 63, 1, 1, 1, 1, 1, 1, 5, 15, 1, 2, 1, 1, 34, 1, 1, 4, 1, 1, 1, 14, 1, 1, 1, 1, 1, 1, 1, 3, 1, 40, 3, 1, 1, 3, 1, 1, 1, 1, 2, 1, 1, 1, 3, 26, 1, 1, 5, 2, 24, 2, 1, 1, 4, 1, 6, 12, 3, 2, 1, 2, 1, 2, 70, 1, 1, 2, 3, 1, 2, 1, 1, 3, 5, 6, 1, 2, 1, 2, 8, 1, 21, 1, 2, 1, 1, 1, 1, 284, 1, 1, 1, 19, 6, 1, 122, 1, 24, 9, 3, 7, 1, 1, 1, 1, 5, 3, 5, 2, 4, 2, 30, 75, 2, 2, 5, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 19, 2, 5, 1, 3, 4, 3, 5, 1, 1, 1, 1, 1, 6, 1, 4, 7, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 9, 1, 19, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 57, 19, 5, 1, 10, 1, 1, 220, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 10, 89, 52, 1, 91, 1, 2, 6, 1, 67, 1, 2, 2, 1, 2, 3, 3, 1, 1, 1, 6, 1, 1, 1, 5, 1, 1, 5, 2, 56, 6, 1, 1, 9, 1, 7, 1, 5, 2, 14, 2190, 1, 1, 1, 36, 1, 1, 1, 4, 8, 1, 1, 7, 1, 1, 1, 1, 6, 1, 1, 1, 14, 1, 2, 1, 2, 1, 1, 2, 1, 3, 4, 1, 1, 1, 1, 12, 4, 11, 2, 7, 1, 28, 2, 1, 1, 1338, 1, 1, 116, 1, 2, 1, 2, 259, 2, 198, 1, 1, 5, 7, 8, 18, 3, 1, 1, 16, 1, 1, 1, 1, 1, 3, 1, 1, 517, 102, 5, 190, 10, 1, 2, 1, 1, 2, 5, 8, 1, 17, 13, 96, 1, 1, 1, 1, 4, 1, 1, 2, 1, 2, 11, 78, 1, 1, 1, 1, 26, 2, 2, 1, 1, 1, 7, 44, 1, 160, 1, 1, 1, 1, 3, 1, 2, 1, 1, 1, 1, 1, 2, 12, 39, 3, 10, 36, 1, 9, 1, 1, 1, 40, 1, 1, 1, 138, 1, 1, 5, 81, 3, 3, 2, 1, 1, 8, 7, 7, 1, 1, 1, 3, 27, 657, 1, 1, 7, 1, 1, 1, 17, 1, 1, 13, 1, 16, 3, 1, 1, 3, 2, 21, 8, 1, 1, 2, 1, 2, 2, 1, 40, 1, 2, 124, 1, 3, 1, 1, 1, 4, 1, 1, 2, 30, 1, 1, 4, 4, 1, 1, 2, 1, 1, 4, 1, 43, 1, 1, 4, 1, 14, 1, 1, 4, 1, 2, 2, 8, 1, 7, 1, 1, 2829, 1, 3, 1, 10, 4, 1, 2, 1, 1, 3, 1, 1, 31, 2, 10, 1, 1, 1, 1, 1, 1, 1, 7, 247, 1, 2, 147, 1, 1, 10, 1, 1, 1, 2, 1, 2, 74, 1, 1, 48, 3, 19, 1, 2754, 11, 3, 1, 1, 2, 1, 6, 20, 63, 1, 1, 1, 1, 1, 3, 3, 16, 2578, 1, 1, 1, 32, 12, 1, 6, 11, 1, 182, 1, 4, 7, 1, 1, 1, 1, 4, 1, 97, 17, 1, 2, 1, 1, 1, 1, 6, 3, 5, 1, 1, 1, 1, 16, 93, 51, 1, 1, 1, 6, 2, 1, 1, 1, 78, 1, 2, 1, 3, 1, 2, 2366, 1, 18, 1, 1, 1, 1, 3, 2, 45, 8, 1, 2, 1, 16, 6, 1, 1, 140, 2, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 3, 1, 3, 1, 1, 2, 1, 6, 1, 1, 118, 1, 1, 4, 1, 1, 5, 1, 39, 4, 1, 1, 1, 1, 1, 1, 1, 1, 15, 1, 5, 1, 14, 1, 2, 8, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 3, 1, 3, 77, 1, 858, 34, 2, 848, 169, 1, 1, 2, 8, 1, 10, 44, 1, 2, 1, 1, 1, 2, 2, 201, 10, 4, 1, 10, 1, 1, 1, 1, 1, 1, 1, 541, 3, 1, 5, 9, 1, 1, 1, 14, 2, 1, 1, 1, 1, 1, 601, 1, 7, 145, 3, 1, 1, 25, 5, 5, 3, 6, 1, 1, 27, 1, 1, 23, 1, 34, 3, 76, 1, 2, 1, 6, 10, 3, 1, 1, 2, 1, 3, 38, 3, 1, 13, 1, 85, 1, 1, 24, 4, 18, 1193, 1, 1, 1, 1, 2, 1, 2, 12, 3, 1, 1, 3, 2, 3, 3, 8, 18, 1, 4, 1, 16, 1, 7, 1, 1, 2, 1, 2, 1, 8, 1, 1, 1, 1, 1, 1, 1, 7, 4, 37, 4, 1, 1, 2, 1, 1, 1, 1, 195, 2, 1, 1, 3, 1, 1, 2, 2, 4, 2, 9, 6, 1, 1, 1, 1, 4, 1, 13, 5, 2, 1, 2, 1, 6, 2, 11, 3, 334, 639, 1, 2, 1, 1, 3, 1, 2, 1, 3, 27, 9, 1, 4, 1, 2, 11, 14, 2, 10, 1, 195, 1, 2, 1, 2, 37, 1, 2, 1, 1, 1, 1, 1, 4, 1, 5, 1, 2, 143, 6, 1, 2, 1, 1, 7, 1, 1, 1, 1, 22, 1, 1, 1, 2, 1, 1, 1, 1, 1, 9, 90, 1, 3, 1, 1, 1, 1, 2, 106, 2, 3, 2, 1, 1, 1, 7, 1, 5, 8, 1, 1, 67, 12, 1, 1, 2, 1, 3, 1, 295, 1, 1, 1, 2, 984, 1, 1, 1, 1, 2, 1, 218, 1, 1, 3, 5, 2, 3, 2, 1, 1, 1, 2, 1, 6, 1, 1, 4, 1, 1, 1, 4, 1, 1, 26, 1, 1, 7, 1, 7, 3, 1, 8, 2, 2, 1, 1, 1, 1, 1, 8, 1, 1, 1, 10, 1, 24, 14, 1, 6, 4, 1, 1, 1, 107, 1, 1, 2, 1, 3, 1, 1, 6, 6, 5, 1, 1, 3, 1, 4, 8, 1, 1, 1, 2, 1, 1, 1, 28, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 15, 1, 64, 2, 30, 2, 1, 1, 163, 1, 1, 1, 2, 1, 1, 40, 1, 1, 5, 510, 1, 2, 1, 1, 1, 9, 1, 4, 3, 1, 1, 2, 1, 309, 1, 5, 21, 1, 1, 1, 1, 1, 4, 1, 1, 35, 1, 1, 24, 1, 7, 2, 1, 1, 2, 2, 1, 1, 16, 35, 1, 1, 3, 20, 1, 2, 5, 1, 9, 1, 1, 362, 9, 7, 1, 42, 1, 2, 1, 1, 1, 2, 1, 2, 1, 5, 5, 1, 3, 1, 1, 1, 2, 1, 1, 1, 5, 2646, 221, 1, 27, 9, 2, 1, 1, 3, 1, 13, 1, 26, 1, 1, 1, 26, 7, 1, 1, 13, 1, 1, 1, 1, 1, 2, 8, 1, 2, 1, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 22, 1, 1, 1, 1, 1, 1, 1, 1, 8, 1, 1, 10579, 1, 2, 48, 135, 1, 1, 1, 2, 1, 1, 4, 1, 1, 1, 1744, 6, 1, 2, 129, 1, 21, 1, 8, 4, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 383, 1, 6, 1, 1, 1, 1, 1, 1, 1, 2, 1, 28, 66, 9, 29, 15, 2, 2, 16, 1, 1, 202, 6, 2, 2, 1, 1, 1, 1, 1, 1, 27, 1, 1, 1, 1, 1, 1, 3, 4, 1, 126, 1, 1, 24, 1, 1, 1, 3, 1, 1, 1, 8, 2, 3, 1, 1, 1, 9, 1, 1, 2, 1, 1, 3, 1, 4, 1, 1, 1, 1, 2, 1, 21, 447, 1, 11, 1, 1, 14, 16, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 4, 1, 2, 1, 1, 2, 52, 1, 1, 2, 1, 2, 4, 1, 1, 1, 3426, 1, 1, 6, 1, 1, 1, 2, 2, 4, 1, 6, 48, 1, 2, 1, 1, 4, 3, 1, 1, 2, 1, 3, 10, 1, 1, 37, 2, 364, 5, 1, 1, 1, 1, 2, 1, 1, 1, 9, 1, 2, 2, 1, 13, 3, 1, 39, 1, 1, 1, 4, 7, 3, 1, 1, 1, 16, 8, 2, 1, 40, 1, 1, 535, 1, 1, 1, 1, 1, 1, 1, 4, 1, 2, 3, 1]
Estimated r: 0.7761158801426653
Estimated p: 0.4429037981403129
../../_images/da0ab5c767fe7cc98af63bc3831e4c3c607fcb8c30284fb529a5c783b0578168.png

  1. The code below performs numsims simulations of SIR disease in a population of \(N=10\) individuals, assuming a negative-binomial offspring distribution with \(r=1\), \(p=0.3\). It plots a histogram of the resulting outbreak sizes. Let \(q_\ell\) be the probability of size \(\ell\).

    (a) The values \(q_1\), \(\ldots\), \(q_N\) satisfy a system of \(N\) equations of the form

    \[ 1 = \sum_{\ell=1}^M c_{M,\ell} q_\ell \]

    For this offspring distribution, write down what \(c_{M,\ell}\) is.

    Solution (a)
    \[\begin{split} C_{M,\ell} = \begin{cases} \mu\left(\frac{M-1}{N-1}\right)^{-\ell}\prod_{j=1}^\ell \frac{M-j}{N-j} & \ell\leq M\\ 0 & \ell>M \end{cases} \end{split}\]

    (b) The example in the text used a relatively sophisticated method to calculate \(q_\ell\) because of error accumulation for very small values of \(q_\ell\). However, this smaller population does not have such issues. A simpler matrix solver is used here. Enter your expression into the code below into c_formula, which will be used to create the matrix. Then it uses the linalg basic solve to solve the matrix equation \(C \vec{q} = \vec{1}\). Confirm that the solution matches the simulations.

    Solution (b)

    See code and plot below.

    (c) Briefly discuss how one could use this approach and Maximum Likelihood Estimation to infer \(p\) and \(r\) given some observed outbreaks in several small closed populations.

    Solution (c)

    Following the previous question, given some observed outbreaks and an assumed form of \(\mu(x)\) with some unknown parameters, this gives us a way to calculate the probability of the observations for each set of parameters. That is, this will let us calculate the likelihood of those parameters. Then we just find the parameters that maximize the likelihood.

Hide code cell source

 
import numpy as np
from scipy.special import gammaln
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt

rng = np.random.default_rng()
r = 1
p = 0.3
N=10

numsims = 5000

L=[]
for counter in range(numsims):
    X = 1
    S = N-X
    cumulative = X
    while X>0:
        newX = 0
        for i in range(X):
            U = rng.negative_binomial(r,p)
            for _ in range(U):
                if rng.random() < S/(N-1):
                    S -= 1
                    newX += 1
        cumulative += newX
        X = newX
    L.append(cumulative)

L = np.array(L)
A = L.sum()
bins = np.arange(0.5, N+1, 1)
plt.hist(L, bins = bins, density=True,)
plt.xticks(np.arange(0, N+1, 1))
plt.xlabel('total outbreak size')
plt.ylabel('probability mass function')




def mu(x):
    return (p/(1-(1-p)*x))**r


def c_formula(N, M, ell,):
    """
    Student-defined formula for c_{M, ell}.
    Replace this with the desired expression.
    """
    if ell <= M:
        prod = 1
        for j in range(1,ell):
            prod *= (M-j)/(N-j)
        return prod/mu((M-1)/(N-1))**(ell)
    else: 
        return 0



# Build lower-triangular matrix
C = np.zeros((N, N), dtype=float)

for M in range(1, N+1):
    for ell in range(1, M+1):
        C[M-1, ell-1] = c_formula(N, M, ell)

# Right-hand side
b = np.ones(N)

# Solve C q = b
q = np.linalg.solve(C, b)


x = np.arange(1, N + 1)
#plt.plot(ell, q)
plt.figure()
plt.bar(x, q, width=1.0, align= 'center')
plt.xlabel('total outbreak size - predicted')
plt.ylabel('probability mass function - predicted')
plt.show()
../../_images/9cb369562d35b23c4857eeaf874c712b394ded449163a1e58ee63d305e41000e.png ../../_images/17eb3db4a9a048f4e4b29b3e12e90eb5ab02f017e19b191e8dbd00acae6fb11c.png

  1. Consider the Till model of stem cells with some \(p_0\) and \(p_2=1-p_0\).

    (a) Compare this to the infectious-disease model in continuous-time. Explain why this is analagous and use the results we have about the continuous-time infectious disease model to calculate the probability that the Till model comes to completion with exactly \(\ell\) differentiated cells (and zero stem cells) in terms of \(p_0\) and \(p_2\). It will help if you make the offspring distribution PGFs identical and use this to determine \(\mathcal{R}_0=\beta/\gamma\) and \(1+\mathcal{R}_0\) in terms of \(p_0\) and \(p_2\).

    Solution (a)

    In the Till model, a stem cell either converts to an end cell or it duplicates itself. Converting to an end cell is mathematically equivalent to recovering, while duplicating is mathematically equivalent to causing an infection in the disease models.

    Then \(p_0\) corresponds to \(\gamma/(\beta+\gamma)\) and \(p_2\) corresponds to \(\beta/(\beta+\gamma)\). So \(\mathcal{R}_0 = \beta/\gamma\) corresponds to \(p_2/p_0\) and \(1+\mathcal{R}_0\) corresponds to \(1 + (p_2/p_0) = (p_0+p_2)/p_0 = 1/p_0\).

    Using the formula for the extinction probability with exactly \(\ell\) infections from Theorem 4.2, we find that

    \[\begin{align*} \mathbb{P}[\text{extinction with $\ell$ end cells}] &= \frac{1}{\ell}\frac{\mathcal{R}_0^{\ell-1}}{(\mathcal{R}_0+1)^{2\ell-1}} \binom{2\ell-2}{\ell-1}\\ &= \frac{1}{\ell}\frac{(p_2/p_0)^{\ell-1}}{(1/p_0)^{2\ell-1}} \binom{2\ell-2}{\ell-1} \\ &= \frac{1}{\ell}p_2^{\ell-1} p_0^\ell \binom{2\ell-2}{\ell-1} \end{align*}\]

    Now we will explore an alternate way to derive the final size that does not depend on the Cycle Lemma proof.

    Let \(\Phi_g(x,y)\) be the PGF for the number of stem cells and differentiated cells at generation \(g\). Assume that differentiated cells persist for each generation (or equivalently, that at each generation they are replaced by an identical differentiated cell).

    (b) Find the PGF of the offspring distribution for both stem cells and differentiated cells.

    Solution (b)
    \[ \mu(x,y) = p_0 y + p_2 x^2 \]

    (c) Using the backward Kolmogorov approach, derive a recursive formula for \(\Phi_g(x,y)\).

    Solution (c)

    We first find \(\Phi_1(x,y)=\mu(x,y)\). We can define \(\Phi_g(x,y|(1,0)) = \Phi_g(x,y)\) to be the PGF if the first cell is a stem cell and \(\Phi_g(x,y|(0,1)) = y\) to be the PGF if the first cell is an end cell. Then

    \[\begin{align*} \Phi_{g+1}(x,y) &= \mu(\Phi_g(x,y|(1,0)), \Phi_g(x,y|(0,1)))\\ &= \mu(\Phi_g(x,y), y)\\ &= p_0y + p_2 \Phi_g(x,y)^2 \end{align*}\]

    (d) What information is contained in the coefficients of the function \(\omega_g(y) = \Phi_g(0,y)\)?

    Solution (d)

    The coefficient of \(y^n\) is the probability of \(n\) end cells and \(0\) stem cells.

    (e) What recursive relation does \(\omega_g(y)\) solve? Assume \(\omega_g(y) \to \omega(y)\) converges as \(g \to \infty\). What quadratic (in \(\omega\)) equation does the recursive relation become in this limit?

    Solution (e)

    We have

    \[\begin{align*} \omega_{g+1}(y) &= \Phi_{g+1}(0,y)\\ &= p_0 y + p_2 \Phi_{g}(0,y)^2\\ &= p_0 y + p_2 \omega_g(y)^2 \end{align*}\]

    As \(g \to \infty\), we conclude \(\omega(y) = p_0y + p_2 \omega(y)^2\)

    (f) Solve the quadratic equation for \(\omega(y)\). There should be two options. Based on our knowledge of the system, what is \(\omega(0)\)? Use this to choose the correct result. Your result should have a \(\sqrt{1-4p_0p_2y}\) in it,

    Solution (f)

    We have

    \[ 0 = p_2\omega(y)^2 +(-1) \omega(y) + p_0y \]

    Using the quadratic formula for \(\omega(y)\) we find

    \[\begin{align*} \omega(y) &= \frac{- (-1) \pm \sqrt{(-1)^2- 4 p_2 p_0y}}{2p_2}\\ &= \frac{1 \pm \sqrt{1-4p_0p_2y}}{2p_2} \end{align*}\]

    So

    \[\begin{align*} \omega(0) &= \frac{1\pm\sqrt{1}}{2p_2}\\ &= 0 \text{ or } 1/p_2 \end{align*}\]

    But a positive value of \(\omega(0)\) corresponds to a scenario with no cells at all. This is not possible, so we must have \(\omega(0)=0\). Thus we conclude

    \[ \omega(y) = \frac{1 - \sqrt{1-4p_0p_2y}}{2p_2} \]

    (g) Write out the first 4 terms of \((1-4p_0p_2y)^{1/2}\) using the binomial theorem (for non-integer exponents):

    \[ (a+b)^r = a^r + ra^{r-1}b + \frac{r(r-1)}{2!} a^{r-2}b^2 + \frac{r(r-1)(r-2)}{3!} a^{r-3}b^3 + \cdots \]

    I don’t see a quick way for me to guide you to the form of the general coefficient, but if you can find it, that would be useful.

    Solution (g)
    \[\begin{align*} \left(1 - 4 p_0p_2\right)^{1/2} &= 1^{1/2} + \frac{1}{2} 1^{-1/2} (-4p_0p_2y) + \frac{(1/2)(-1/2)}{2} 1^{-3/2} (-4p_0p_2y)^2 + \frac{(1/2)(-1/2)(-3/2)}{6} 1^{-5/2} (-4p_0p_2y)^3 + \cdots\\ &= 1 - 2 p_0p_2y - 2p_0^2p_2^2y^2 - 4 p_0^3p_2^3y^3 + \cdots \end{align*}\]

    (h) Use this to write out the first 3 terms of the expansion for \(\omega(y)\). Verify that they agree with the result calculated in (a).

    Solution (h)
    \[\begin{align*} \omega(y) &= \frac{1 - \sqrt{1-4p_0p_2y}}{2p_2}\\ &= \frac{1 - (1 - 2 p_0p_2y - 2p_0^2p_2^2y^2 - 4 p_0^3p_2^3y^3 + \cdots )}{2p_2}\\ &= \frac{+ 2p_0p_2y + 2p_0^2p_2^2y^2 + 4 p_0^3p_2^3y^3 + \cdots}{2p_2}\\ &= p_0y + p_0^2p_2y^2 + 2 p_0^3p_2^2 y^3 + \cdots \end{align*}\]

    Previously we calculated that the probability is \(\frac{1}{\ell}p_2^{\ell-1} p_0^\ell \binom{2\ell-2}{\ell-1}\). Setting \(\ell =1,2,3\) we find \(1p_0 \binom{0}{0} = p_0\), \(p_2p_0^2\binom{2}{1}/2 = p_2p_0^2\), and \(p_2^2p_0^3\binom{4}{2}/3= 2 p_2^2p_0^3\). These correspond to the coefficients we just found using the binomial expansion.

  2. Consider a continuous-time version of the Macken-Perelson model. Assume stem cells divide with rate \(r_s\) and the mast cells divide with rate \(r_m\). Assume the system begins with a single stem cell. Let \(\Phi(x,y,z,t)\) be the PGF for the number of each type at time \(t\).

    (a) Derive a backward Kolmogorov equation for \(\Phi(x,y,z,t)\)

    Solution (a)

    We introduce \(\Phi_s(x,y,z,t)\), \(\Phi_m(x,y,z,t)\), and \(\Phi_e(x,y,z,t)\) to be the PGFs assuming an initial stem, mast or end cell.

    We can quickly see \(\Phi_e(x,y,z,t) = z\).

    The offspring distribution of a stem cell is given by \(\xi_s(x,y,z) = (p_sx + (1-p_s) y)^2\) and the offspring distribution of a mast cell is given by \(\xi_m(x,y,z) = (p_my+(1-p_m)z)^2\).

    At time \(\Delta t\), we can conclude

    \[\begin{align*} \Phi_s(x,y,z,\Delta t) &= (1-r_s\Delta t) x + r_s\Delta t \xi_s(x,y,z)\\ \Phi_m(x,y,z,\Delta t) &= (1-r_m\Delta t) y + r_m \Delta t \xi_m(x,y,z)\\ \Phi_e(x,y,z,\Delta t) &= z \end{align*}\]

    Following our standard argument, it follows that

    \[\begin{align*} \Phi_s(x,y,z,t+\Delta t) &= (1-r_s\Delta t) \Phi_s(x,y,z,t) + r_s\Delta t \xi_s(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t))\\ \Phi_m(x,y,z,t+\Delta t) &= (1-r_m\Delta t) \Phi_m(x,y,z,t) + r_m \Delta t \xi_m(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t))\\ \Phi_e(x,y,z,t+\Delta t) &= \Phi_e(x,y,z,t) \end{align*}\]

    So

    \[\begin{align*} \frac{\partial}{\partial t}\Phi_s(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_s(x,y,z,t+\Delta t)-\Phi_s(x,y,z,t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{(-r_s\Delta t) \Phi_s(x,y,z,t) + r_s\Delta t \xi_s(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t))}{\Delta t}\\ &= r_s\xi_s(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t)) - r_s \Phi_s(x,y,z,t)\\ \frac{\partial}{\partial t} \Phi_n(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_m(x,y,z,t+\Delta t)-\Phi_m(x,y,z,t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{(-r_m\Delta t) \Phi_m(x,y,z,t) + r_m \Delta t \xi_m(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t))}{\Delta t}\\ &= r_m\xi_m(\Phi_s(x,y,z,t),\Phi_m(x,y,z,t),\Phi_e(x,y,z,t)) - r_m \Phi_m(x,y,z,t)\\ \frac{\partial}{\partial t} \Phi_e(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_e(x,y,z,t+\Delta t)-\Phi_e(x,y,z,t)}{\Delta t}\\ &= 0 \end{align*}\]

    (b) Derive a forward Kolmogorov equation for \(\Phi(x,y,z,t)\).

    Solution (b)

    We have \(\Phi_s(x,y,z,t)\), \(\Phi_m(x,y,z,t)\) and \(\Phi_e(x,y,z,t)\). To find them at a time \(\Delta t\) later, we substitute \((1-r_s\Delta t)x + r_s\Delta t \xi_s(x,y,z)\) for \(x\), \((1-r_m\Delta t)y + r_m\Delta t \xi_s(x,y,z)\) for \(y\), and \(z\) for \(z\). We have

    \[\begin{align*} \Phi_s(x,y,z,t+\Delta t) &= \Phi_s((1-r_s\Delta t)x + r_s\Delta t \xi_s(x,y,z),(1-r_m\Delta t)y + r_m\Delta t \xi_s(x,y,z),z,t, z)\\ &= \Phi_s(x,y,z,t) + (r_s \Delta t)(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_s(x,y,z,t) + (r_m \Delta t)(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_s(x,y,z,t)\\ \Phi_m(x,y,z,t+\Delta t) &= \Phi_e((1-r_s\Delta t)x + r_s\Delta t \xi_s(x,y,z),(1-r_m\Delta t)y + r_m\Delta t \xi_s(x,y,z),z,t, z)\\ &= \Phi_m(x,y,z,t) + (r_s \Delta t)(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_m(x,y,z,t) + (r_m \Delta t)(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_m(x,y,z,t)\\ \Phi_e(x,y,z,t+\Delta t)&= z \end{align*}\]

    So

    \[\begin{align*} \frac{\partial}{\partial t}\Phi_s(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_s(x,y,z,t+\Delta t)-\Phi_s(x,y,z,t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{(r_s \Delta t)(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_s(x,y,z,t) + (r_m \Delta t)(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_s(x,y,z,t)}{\Delta t}\\ &= r_s(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_s(x,y,z,t) + r_m(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_s(x,y,z,t)\\ \frac{\partial}{\partial t} \Phi_m(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_m(x,y,z,t+\Delta t)-\Phi_m(x,y,z,t)}{\Delta t}\\ &= \lim_{\Delta t \to 0} \frac{(r_s \Delta t)(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_m(x,y,z,t) + (r_m \Delta t)(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_m(x,y,z,t)}{\Delta t}\\ &= r_s(\xi_s(x,y,z)-x)\frac{\partial}{\partial x} \Phi_m(x,y,z,t) + r_m(\xi_m(x,y,z)-y)\frac{\partial}{\partial y} \Phi_m(x,y,z,t)\\ \frac{\partial}{\partial t} \Phi_e(x,y,z,t) &= \lim_{\Delta t \to 0} \frac{\Phi_e(x,y,z,t+\Delta t)-\Phi_e(x,y,z,t)}{\Delta t}\\ &= 0 \end{align*}\]

    (c) How can we interpret \(\Phi(0,0,1,t)\)? Use one of the Kolmogorov equations to write down an ordinary differential equation (that is, just \(t\) derivatives) that we can use to solve for this. What is the appropriate initial condition?

    Solution (c)

    This is the probability that at time \(t\) the only cells are end cells.

    Since we only allow \(t\) derivatives, this needs to use the Backward Kolmogorov Equations:

    \[\begin{align*} \frac{\partial}{\partial t} \Phi_s&= r_s(\xi_s(\Phi_s,\Phi_m,\Phi_e) - \Phi_s)\\ \frac{\partial}{\partial t} \Phi_m&= r_m(\xi_m(\Phi_s,\Phi_m,\Phi_e) - \Phi_m)\\ \frac{\partial}{\partial t} \Phi_e&= 0 \end{align*}\]

    The appropriate initial condition is \(\Phi_s=0\), \(\Phi_m=0\), and \(\Phi_e=1\).

  3. Write a 1+ page summary of one of the following papers.

  • Explain the basic problem the researchers are studying.

  • What did the PGF in the paper represent?

  • Explain the equations governing their PGF (how did they determine how it changes in time? What is its initial condition?)

  • How did they use it to analyze their problem?

  • What were their conclusions?

Here are the papers to choose from:

  • Fixation Probabilities When Generation Times Are Variable: The Burst–Death Model

  • Gestational mutations and carcinogenesis

  • Stochasticity, invasions, and branching random walks

  • Quantifying uncertainty in a predictive model for popularity dynamics

  • Superspreading and the effect of individual variation on disease emergence

  • Second look at the spread of epidemics on networks

  • https://arxiv.org/pdf/0707.0080

  • A generating function perspective on the transmission forest

If there is a different paper you would like to study, please let me know.

Many of these define a PGF in terms of \(\mathbb{E}[x^S]\) where \(S\) is a random (non-negative valued variable). This may look unusual to you as I have not used that notation, but this is simply \(\mathbb{E}[x^S] = \sum_k \mathbb{P}[S=k] x^k = \sum_k p_k x^k\).