20200507.epidemic-attack-rate-by-variance.sage

## Plot the attack rate of an epidemic for a fixed reproduction number (set as
## k=2.5 below) in function of the square root of the variance in infectious
## contacts.  The blue plot is for a directed graph (size of the in-component
## of a test node in a directed graph whose distribution of in-degrees is as
## specified and that of out-degrees is irrelevant: this is computed as the
## non-extinction probability of a Galton-Watson process over the in-degree
## distribution); the red plot is for an undirected graph with the
## corresponding degree distribution (see reference below).

## A binomial distribution with mean k and n tries has variance k*(n-k)/n and
## generating function ((n-k+k*z)/n)^n.

## A negative binomial distribution with mean k and r failures has variance
## k*(r+k)/r and generating function (r/(r+k-k*z))^r.

## The limiting case is Poisson, which for mean k has variance k and generating
## function exp(k*(1-z)).

## The extinction probability of a Galton-Watson process associated to
## a probability distribution with generating function G(z) is the
## smallest fixed point p of G.  The epidemic attack rate is the
## non-extinction probability, namely 1-p.

## For the size of the giant connected component an undirected graph,
## following <https://arxiv.org/abs/2002.04004>, we take the smallest
## fixed point q of G1 where G1(z) = G'(z)/G'(1) (namely
## ((n-k+k*z)/n)^(n-1) for binomial, and (r/(r+k-k*z))^(r+1) for
## negative binomial), and the epidemic attack rate is 1-G(q).

k = 2.5
critval = N(-lambert_w(-k*exp(-k))/k)  ## Solution of z == exp(k*(1-z))
var('z')

tab = []
tab.append((N(sqrt(k)), N(1-critval)))
for i in range(1, 501):
    sigma = (2*k*i/500)
    v = sigma^2
    if abs(v-k)<0.1:
        pass
    elif v<k:
        n = N(k^2/(k-v))
        ## v == N(k*(n-k)/n)
        tab.append((sigma, 1-find_root(((n-k+k*z)/n)^n-z, 0, 0.999)))
    elif v>k:
        r = N(k^2/(v-k))
        ## v == N(k*(r+k)/r)
        tab.append((sigma, 1-find_root((r/(r+k-k*z))^r-z, 0, 0.999)))
tab.sort()

tab2 = []
tab2.append((N(sqrt(k)), N(1-critval)))
for i in range(1, 501):
    sigma = (2*k*i/500)
    v = sigma^2
    if abs(v-k)<0.1:
        pass
    elif v<(2-k)*k:
        tab2.append((sigma, 0))
    elif v<k:
        n = N(k^2/(k-v))
        ## v == N(k*(n-k)/n)
        q = find_root(((n-k+k*z)/n)^(n-1)-z, 0, 0.999)
        tab2.append((sigma, 1-(((n-k+k*z)/n)^n).subs({z:q})))
    elif v>k:
        r = N(k^2/(v-k))
        ## v == N(k*(r+k)/r)
        q = find_root((r/(r+k-k*z))^(r+1)-z, 0, 0.999)
        tab2.append((sigma, 1-((r/(r+k-k*z))^r).subs({z:q})))
tab2.sort()

list_plot(tab, plotjoined=True, ymin=0, ymax=1) + list_plot(tab2, plotjoined=True, ymin=0, ymax=1, color="red")