Age of the Youngest Infection (AoY) • pf.memory

library(deSolve)
library(pf.memory)

We want to derive a random variable describing the age of the youngest infection (AoY). The derivations are presented below:

We compute the distribution function first:

$F_Y(a) \sim \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{1-e^{-m_\tau(a)}} = \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{p_\tau(a)}$

We differentiate to get a formula for the density function:

$Y_\tau(a) \sim f_Y(\alpha; a, \tau) = f_A(\alpha, a,\tau) e^{-m_\tau(a) F_A(\alpha, a,\tau)} \frac{m_\tau(a)}{p_\tau(a)}$

The mean AoY is:

$\left< Y_\tau(a) \right> = \int_0^a \alpha \; f_Y(\alpha, a, \tau) \; d\alpha$

And the higher order moments for the AoY are:

$\left< Y_\tau(a)^n \right> = \int_0^n \alpha^n \; f_y(\alpha, a, \tau) \; d\alpha$

AoY Density

aoy = dAoY(1:1825, 5*365, foiP3)
raoy = rAoY(10^5, 5*365, foiP3)

hist(raoy, breaks=c(0:(5*365)), right=F, plot=F)->daoy
plot(daoy$mids, daoy$density, type = "l")
lines(1:1825, aoy, type = "l", col = "red")

Moments

aa = seq(5, 5*365, by = 5) 
moment1y = momentAoY(aa, foiP3)
moment2y = momentAoY(aa, foiP3, n=2)
moment3y = momentAoY(aa, foiP3, n=3)

The first three moments of the AoY plotted over time. In the top plot, we’ve also plotted the $n^{th}$ root of the $n^{th}$ moment.

par(mfrow = c(4,1), mar = c(0.5, 4, 0.5, 2))
plot(aa, moment1y, type = "l", xlab = "", ylab = expression(E(Y)), lwd=2, xaxt="n", ylim = range((moment3y)^(1/3)) ) 
lines(aa, sqrt(moment2y), col = "darkgreen")
lines(aa, (moment3y)^(1/3), col = "purple")
plot(aa, moment2y, type = "l", xlab = "", ylab = expression(E(Y^2)), lwd=2, xaxt="n", col = "darkgreen") 
plot(aa, moment3y, type = "l", xlab = "", ylab = expression(E(Y^2)), lwd=2, col = "purple") 
mtext("Age (in Days)", 1, 3)

Derivation

First, we note that the CDF for the minimum of $n$ random values of AoI in some cohort is:

$1-(1-F_A(\alpha,a,\tau))^n$ Using this, we can compute the CDF for the full distribution of AoY across the population:

$F_Y(\alpha, a, \tau) = \frac{1}{p_\tau(a)} \sum_{\zeta>0} \mbox{dpois}\left(\zeta,m_\tau(a)\right) \left( 1-(1-F_A(\alpha,a,\tau))^\zeta \right)$

or rewriting slightly:

$F_Y(\alpha, a, \tau) p_\tau(a) = \sum_{\zeta>0} \frac{m_\tau(a)^\zeta}{\zeta!} e^{-m_\tau(a)} \left(1- \left(1-F_A \left(\alpha,a, \tau\right)\right)^\zeta\right)$

pulling out $e^{-m_\tau(a)}$ and recombining terms, we get:

$F_Y (\alpha, a, \tau) p_\tau(a)= p_\tau(a) - e^{-m_\tau(a)} \sum_{\zeta>0} \frac{ \left[m_\tau(a) \left(1-F_A \left(\alpha,a, \tau\right)\right)\right]^\zeta}{\zeta!}$ if the sum were from $\zeta=0\ldots\infty,$ the expression on the inside would be simple, but the sum is for $M>0$ , so we add the term for $M=0$ into the sum, subtract it outside:

$F_Y(\alpha, a,\tau) p_\tau(a) = p_\tau(a) - e^{-m_\tau(a)} \left(-1+ \sum_{\zeta>=0} \frac{ \left[m_\tau(a) \left(1-F_A \left(\alpha,a\right)\right)\right]^\zeta}{\zeta!} \right)$

but $e^{-m_\tau(a)} = 1-p_\tau(a)$ and the infinite sum is simply an exponential, so:

$F_Y(\alpha, a, \tau) p_\tau(a) = p_\tau(a) + 1-p_\tau(a) - e^{-m_\tau(a)} e^{m_\tau(a) (1-F_A(\alpha, a, \tau))}$ so

$F_Y(\alpha, a,\tau) p_\tau(a) = 1-e^{-m_\tau(a)F_a(\alpha, a)}$

and

$F_Y(\alpha, a,\tau) = \frac{1-e^{-m_\tau(a)F_A(\alpha, a,\tau)}}{p_\tau(a)}.$

The PDF for the AoY is found by taking the derivative of $F_Y$ :

$\begin{array}{rcl} Y \sim f_Y(\alpha, a, \tau) & = & \frac{d F_Y(\alpha, a,\tau)}{d\alpha} \\ &=& \; f_A(\alpha, a) e^{-m_\tau(a) F_A(\alpha, a,\tau)} \frac{\textstyle{m_\tau(a)}}{\textstyle{p_\tau(a)}}\\ \end{array}$