Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
## 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment