Skip to contents

Age of the Youngest Infection (AoY)

AoY.Rmd

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} \]