20200511.epidemic-attack-rate-by-variance-2.sage

## Plot the attack rate and (some form of) herd immunity threshold for a very
## trivial epidemic model, namely for a population divided in an infinite
## number of subpopulations with varying reproduction number and density given
## by a gamma distribution around a specified mean.  That is, we're just
## averaging the SIR attack rate and herd immunity threshold by a gamma
## distribution around a specified mean.

k=2.5 ## Mean reproduction number
tab=[]
for i in range(501):
    sigma = N(2*k*i/500)
    v = sigma^2
    if i==0:
        tab.append((0, (1+lambert_w(-x*exp(-x))/x).subs({x:k})))
    elif i<10:
        pass
    else:
        scal = v/k
        shap = k^2/v
        igrd = (-lambert_w(-x*exp(-x))/x)
        dstr = (x^(shap-1) * exp(-x/scal) / (gamma(shap)*(scal^shap)))
        expr = igrd * dstr
        tbas = N(gamma(shap,1/scal)/gamma(shap))  ## From 1 to infinity
        tbs1 = N(gamma(shap,shap)/gamma(shap))    ## From k to infinity
        tbs0 = tbas - tbs1                        ## From 1 to k
        ubnd = k
        igr1 = 0
        ibs1 = 0
        while abs(tbs1-ibs1)>1.e-4 and dstr.subs({x:ubnd}) > 1.e-8:
            lbnd = ubnd
            fctr = N(2)
            while True:
                try:
                    ubnd = fctr*lbnd
                    # print "DEBUG: integrating from",lbnd,"to",ubnd,"; fctr=",fctr,"; ibs1=",ibs1,"; igr1=",igr1
                    igr1 += expr.nintegral(x,lbnd,ubnd)[0]
                    ibs1 += dstr.nintegral(x,lbnd,ubnd)[0]
                    break
                except TypeError:
                    fctr = sqrt(fctr)
        lbnd = k
        igr0 = 0
        ibs0 = 0
        while abs(tbs0-ibs0)>1.e-4 and dstr.subs({x:lbnd}) > 1.e-8 and lbnd > 1:
            ubnd = lbnd
            fctr = N(2)
            while True:
                try:
                    lbnd = max(1,ubnd/fctr)
                    # print "DEBUG: integrating from",lbnd,"to",ubnd,"; fctr=",fctr,"; ibs0=",ibs0,"; igr0=",igr0
                    igr0 += expr.nintegral(x,lbnd,ubnd)[0]
                    ibs0 += dstr.nintegral(x,lbnd,ubnd)[0]
                    break
                except TypeError:
                    fctr = sqrt(fctr)
        ibas = ibs0+ibs1
        igrl = igr0+igr1
        tab.append((sigma, ibas-igrl))
tab.sort()

tab2=[]
for i in range(501):
    sigma = N(2*k*i/500)
    v = sigma^2
    if i==0:
        tab2.append((0, 1-1/k))
    else:
        scal = v/k
        shap = k^2/v
        tab2.append((sigma, N((gamma(shap,1/scal)-gamma(shap-1,1/scal)/scal)/gamma(shap))))
tab2.sort()

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