An attempt to estimate COVID-19 transmission rate

coronavirus.py

def clamp(x, x_min, x_max):
    return max(x_min, min(x_max, x))


def safe_segment(a, start=0, end=None):
    if end is None:
        end = len(a)
    start = clamp(start, 0, len(a))
    end = clamp(end, start, len(a))
    return a[start:end]


def safe_sum(a, start=0, end=None):
    return sum(safe_segment(a, start, end))


class State(object):
    report_period = 7
    transmission_rate = 2.5  # Just ignore it - the code below overwrites it anyway
    transmission_period = 7
    death_rate = 0.05
    death_period = 12
    cure_period = 28

    def __init__(self):
        self.infections = []
        self.deaths = []

    def __str__(self):
        ri = int(self.reported_infected_count())
        i = int(self.infected_count())
        c = int(self.cured_count())
        d = int(self.dead_count())
        return f'Day #{self.day_index()}: reported={ri} cured={c} died={d} (infected={i})'

    def day_index(self):
        return len(self.infections)

    def daily_transmission_rate(self):
        return pow(self.transmission_rate + 1, 1.0 / self.transmission_period) - 1

    def contagious_count(self):
        return safe_sum(self.infections, len(self.infections) - self.transmission_period)

    def reported_infected_count(self):
        return safe_sum(self.infections, 0, len(self.infections) - self.report_period)

    def infected_count(self):
        return sum(self.infections)

    def sick_count(self):
        infected_count_28 = safe_sum(self.infections, len(self.infections) - self.cure_period)
        dead_count_14 = safe_sum(self.infections, len(self.infections) - (self.cure_period - self.death_period))
        return infected_count_28 - dead_count_14

    def dead_count(self):
        return sum(self.deaths)

    def cured_count(self):
        max_day = len(self.infections) - self.cure_period - 1
        cured_count = safe_sum(self.infections, 0, max_day)
        dead_count = safe_sum(self.deaths, 0, max_day)
        return cured_count - dead_count

    def next_state(self, index=1):
        if index == 0:
            return self
        if index > 1:
            s = self.next_state(index - 1)
            return s.next_state()
        death_pool_day = self.day_index() - self.death_period
        death_pool_count = safe_sum(self.infections, death_pool_day, death_pool_day + 1)
        new_infections = self.contagious_count() * self.daily_transmission_rate()
        new_deaths = death_pool_count * self.death_rate

        r = State()
        r.infections = self.infections + [new_infections]
        r.deaths = self.deaths + [new_deaths]
        return r


def find_transmission_rate():
    min_rate = 1
    max_rate = 1000

    while (max_rate - min_rate) > 0.00001:
        next_rate = (min_rate + max_rate) / 2
        State.transmission_rate = next_rate

        s = State()
        s.infections = [1]
        s.deaths = [0]
        while s.reported_infected_count() < 68 and s.day_index() < 200:
            s = s.next_state()
        e = s.next_state(14)

        if e.reported_infected_count() < 3000:
            min_rate = next_rate
        else:
            max_rate = next_rate
    return s, e


def find_transmission_and_death_rate():
    min_rate = 0.001
    max_rate = 0.3

    while (max_rate - min_rate) > 0.001:
        next_rate = (min_rate + max_rate) / 2
        State.death_rate = next_rate

        s, e = find_transmission_rate()
        if e.dead_count() < 60:
            min_rate = next_rate
        else:
            max_rate = next_rate
    return s, e


def model_case(title, report_period, transmission_period, cure_period, death_period):

    print(f'"{title}" case:')
    print(f'  report_period:       {report_period}')
    print(f'  transmission_period: {transmission_period}')
    print(f'  cure_period:         {cure_period}')
    print(f'  death_period:        {death_period}')

    State.report_period = report_period
    State.transmission_period = transmission_period
    State.cure_period = cure_period
    State.death_period = death_period
    s, e = find_transmission_and_death_rate()

    print(f'Computed rates:')
    print(f'  transmission_rate: {State.transmission_rate}')
    print(f'  death_rate:        {State.death_rate}')
    print()

    print('Modelling Feb 29 ... Mar 14:')
    while s.reported_infected_count() <= 3100:
        print(s)
        s = s.next_state()
    print()


model_case('Pessimistic (=15)', 10, 10, 24, 12)
model_case('Realistic (=7.5)', 7, 7, 21, 12)
model_case('Ok-ish (=2.7)', 3, 3, 17, 12)
model_case('Official (=2)', 2, 2, 17, 12)

Output:

"Pessimistic (=15)" case:
  report_period:       10
  transmission_period: 10
  cure_period:         24
  death_period:        12
Computed rates:
  transmission_rate: 15.169240988790989
  death_rate:        0.033119140625

Modelling Feb 29 ... Mar 14:
Day #27: reported=78 cured=1 died=1 (infected=1059)
Day #28: reported=101 cured=1 died=1 (infected=1374)
Day #29: reported=131 cured=2 died=2 (infected=1782)
Day #30: reported=171 cured=3 died=3 (infected=2312)
Day #31: reported=222 cured=4 died=4 (infected=3000)
Day #32: reported=288 cured=5 died=5 (infected=3891)
Day #33: reported=374 cured=7 died=7 (infected=5047)
Day #34: reported=485 cured=9 died=9 (infected=6547)
Day #35: reported=629 cured=12 died=12 (infected=8492)
Day #36: reported=816 cured=16 died=16 (infected=11015)
Day #37: reported=1059 cured=21 died=20 (infected=14288)
Day #38: reported=1374 cured=27 died=27 (infected=18534)
Day #39: reported=1782 cured=35 died=35 (infected=24040)
Day #40: reported=2312 cured=46 died=45 (infected=31183)
Day #41: reported=3000 cured=60 died=59 (infected=40447)

"Realistic (=7.5)" case:
  report_period:       7
  transmission_period: 7
  cure_period:         21
  death_period:        12
Computed rates:
  transmission_rate: 7.524774305522442
  death_rate:        0.075166015625

Modelling Feb 29 ... Mar 14:
Day #23: reported=74 cured=1 died=1 (infected=473)
Day #24: reported=96 cured=1 died=1 (infected=616)
Day #25: reported=126 cured=1 died=2 (infected=802)
Day #26: reported=164 cured=2 died=3 (infected=1044)
Day #27: reported=214 cured=3 died=4 (infected=1360)
Day #28: reported=279 cured=4 died=5 (infected=1770)
Day #29: reported=363 cured=6 died=7 (infected=2304)
Day #30: reported=473 cured=8 died=9 (infected=3000)
Day #31: reported=616 cured=11 died=12 (infected=3904)
Day #32: reported=802 cured=14 died=16 (infected=5082)
Day #33: reported=1044 cured=19 died=20 (infected=6615)
Day #34: reported=1360 cured=25 died=27 (infected=8611)
Day #35: reported=1770 cured=33 died=35 (infected=11208)
Day #36: reported=2304 cured=43 died=46 (infected=14588)
Day #37: reported=3000 cured=56 died=60 (infected=18988)

"Ok-ish (=2.7)" case:
  report_period:       3
  transmission_period: 3
  cure_period:         17
  death_period:        12
Computed rates:
  transmission_rate: 2.746389038860798
  death_rate:        0.217658203125

Modelling Feb 29 ... Mar 14:
Day #17: reported=72 cured=0 died=1 (infected=161)
Day #18: reported=94 cured=0 died=1 (infected=210)
Day #19: reported=123 cured=1 died=2 (infected=275)
Day #20: reported=161 cured=1 died=2 (infected=359)
Day #21: reported=210 cured=2 died=3 (infected=468)
Day #22: reported=275 cured=3 died=5 (infected=611)
Day #23: reported=359 cured=5 died=6 (infected=796)
Day #24: reported=468 cured=7 died=9 (infected=1039)
Day #25: reported=611 cured=10 died=11 (infected=1354)
Day #26: reported=796 cured=13 died=15 (infected=1765)
Day #27: reported=1039 cured=18 died=20 (infected=2301)
Day #28: reported=1354 cured=24 died=26 (infected=2999)
Day #29: reported=1765 cured=31 died=35 (infected=3910)
Day #30: reported=2301 cured=41 died=45 (infected=5096)
Day #31: reported=2999 cured=54 died=59 (infected=6641)

"Official (=2)" case:
  report_period:       2
  transmission_period: 2
  cure_period:         17
  death_period:        12
Computed rates:
  transmission_rate: 1.981666974723339
  death_rate:        0.257369140625

Modelling Feb 29 ... Mar 14:
Day #16: reported=82 cured=0 died=1 (infected=139)
Day #17: reported=107 cured=0 died=1 (infected=180)
Day #18: reported=139 cured=0 died=2 (infected=232)
Day #19: reported=180 cured=1 died=3 (infected=301)
Day #20: reported=232 cured=1 died=4 (infected=389)
Day #21: reported=301 cured=2 died=5 (infected=502)
Day #22: reported=389 cured=4 died=7 (infected=649)
Day #23: reported=502 cured=6 died=9 (infected=837)
Day #24: reported=649 cured=8 died=12 (infected=1081)
Day #25: reported=837 cured=12 died=16 (infected=1396)
Day #26: reported=1081 cured=16 died=21 (infected=1801)
Day #27: reported=1396 cured=21 died=27 (infected=2324)
Day #28: reported=1801 cured=28 died=35 (infected=2999)
Day #29: reported=2324 cured=37 died=46 (infected=3870)
Day #30: reported=2999 cured=48 died=59 (infected=4994)