Rate of change tracking

This whole thing we’re facing has had me focusing my attention on our immediate area. To give me some kind of daily indicator it meant coming up with a simple program.

Workup
data example
# the data
p = '''060420
March  6    1
March  7    1
March  8    1
March  9    7
March 10   16
March 11   24
March 12   27
March 13   35
March 14   55
March 15   64
March 16   90
March 17  103
March 18  135
March 19  168
March 20  212
March 21  250
March 22  289
March 23  338
March 24  411
March 25  478
March 26  530
March 27  626
March 28  679
March 29  708
March 30  758
March 31  890
April  1 1003
April  2 1136
April  3 1210
April  4 1273
April  5 1334
April  6 1348
'''
# the program
def case_delta(u):
    t = tuple(map(lambda x: int(x[2]), map(lambda x: x.split(), u.split('\n')[1:-1])))
    return [f"{(b / a):.3f}" for (a, b) in zip(t[:-1], t[1:])]
# the output
print (case_delta(p))
['1.000', '1.000', '7.000', '2.286', '1.500', '1.125', '1.296', '1.571', '1.164', '1.406', '1.144', '1.311', '1.244', '1.262', '1.179', '1.156', '1.170', '1.216', '1.163', '1.109', '1.181', '1.085', '1.043', '1.071', '1.174', '1.127', '1.133', '1.065', '1.052', '1.048', '1.010']
testing

Now to test the program. What does it represent? A ratio. What better a test than a sequence that represents a near constant ratio than Fibonacci? Enter the next problem. Where do we get a Fibonacci Sequence from?

We write the code that generates it. After all, we have a general term that describes the sequence.

Tn = Tn-1 + Tn-2

n is a natural number, the subscript of T. In order to start the sequence we need the first two terms. A natural sequence will start with natural terms, so 1 and 1 represent the base of the sequence. We’re not counting, remember.

In programming terms that means,

terms.append(terms[-2] + terms[-1])

In Pythonic terms it means,

a, b = b, a + b

if we are generating a new term. All we need to do is iterate enough times to reach the desired term.

def fib_term(t, a=1, b=1):
    for x in range(t - 1):
        a, b = b, a + b
    return b

That gives the terms but not the sequence.

def fib_seq(t):
    return [fib_term(x + 1) for x in range(t)]

gives us the sequence. It will be a list. object in this instance.

Now we test the presumptive ratio in the original code.

def fib_ratio(u):
    return [f"{(b / a):.3f}" for (a, b) in zip(u[:-1], u[1:])]

print (fib_ratio(fib_seq(10)))

Output

['2.000', '1.500', '1.667', '1.600', '1.625', '1.615', '1.619', '1.618', '1.618']

Any more, and due to rounding there would be no change. This, to my thinking represents a done deal. The code works.


Back to the main program. I call it my optimism meter. As long as the numbers are less than 1.2 I’m optimistic. Any trend greater than 1.2 is not optimistic.

Rule of 72

This axiom states that the doubling rate is equivalent to the percentage divided into 72. Twenty into seventy-two means the doubling rate is less than four days. Ergo, the urgency to keep this number smaller than 1.2. 1.02 even better. 1.002 even better still. Ideally, 1.0 on a flat line is what we are counting on.


Today's report from our region...
data

Source: https://covid19stats.alberta.ca/

data = '''310520
March  6    1    0
March  7    1    0
March  8    1    0
March  9    7    0
March 10   16    0
March 11   24    0
March 12   27    0
March 13   35    0
March 14   53    0
March 15   62    0
March 16   90    0
March 17  102    0
March 18  133    0
March 19  167    0
March 20  212    0
March 21  248    0
March 22  289   12
March 23  343   17
March 24  411   21
March 25  473   27
March 26  524   34
March 27  613   54
March 28  665   75
March 29  695   96
March 30  737  121
March 31  872  144
April  1  982  184
April  2 1099  211
April  3 1160  247
April  4 1215  292
April  5 1260  351
April  6 1294  417
April  7 1344  482
April  8 1375  552
April  9 1417  635
April 10 1470  708
April 11 1524  771
April 12 1594  825
April 13 1662  871
April 14 1790  920
April 15 1921  993
April 16 2066 1047
April 17 2282 1101
April 18 2469 1146
April 19 2666 1174
April 20 2855 1210
April 21 3129 1245
April 22 3424 1280
April 23 3772 1328
April 24 4012 1370
April 25 4242 1438
April 26 4449 1511
April 27 4631 1622
April 28 4892 1758
April 29 5129 1908
April 30 5358 2114
May    1 5497 2311
May    2 5616 2485
May    3 5685 2664
May    4 5761 2897
May    5 5824 3172
May    6 5889 3507
May    7 5973 3758
May    8 6051 3967
May    9 6137 4154
May   10 6210 4337
May   11 6258 4609
May   12 6329 4819
May   13 6392 5034
May   14 6453 5168
May   15 6524 5285
May   16 6578 5345
May   17 6624 5421
May   18 6663 5487
May   19 6708 5554
May   20 6736 5609
May   21 6766 5687
May   22 6793 5780
May   23 6826 5848
May   24 6856 5902
May   25 6868 5956 138
May   26 6887 6024 139
May   27 6914 6083 141
May   28 6945 6145 143
May   29 6977 6205 143
May   30 7001 6247 143
May   31 7010 6283 143
'''

If you can access data from your region, compose it like above and run the code.


output...
['1', '1', '7', '2.286', '1.5', '1.125', '1.296', '1.514', '1.170', '1.452', '1.133', '1.304', '1.256', '1.269', '1.170', '1.165', '1.187', '1.198', '1.151', '1.108', '1.170', '1.085', '1.045', '1.060', '1.183', '1.126', '1.119', '1.056', '1.047', '1.037', '1.027', '1.039', '1.023', '1.031', '1.037', '1.037', '1.046', '1.043', '1.077', '1.073', '1.075', '1.105', '1.082', '1.080', '1.071', '1.096', '1.094', '1.102', '1.064', '1.057', '1.049', '1.041', '1.056', '1.048', '1.045', '1.026', '1.022', '1.012', '1.013', '1.011', '1.011', '1.014', '1.013', '1.014', '1.012', '1.008', '1.011', '1.010', '1.010', '1.011', '1.008', '1.007', '1.006', '1.007', '1.004', '1.004', '1.004', '1.005', '1.004', '1.002', '1.003', '1.004', '1.004', '1.005', '1.003', '1.001']
(1, 1, 1, 7, 16, 24, 27, 35, 53, 62, 90, 102, 133, 167, 212, 248, 289, 343, 411, 473, 524, 613, 665, 695, 737, 872, 982, 1099, 1160, 1215, 1260, 1294, 1344, 1375, 1417, 1470, 1524, 1594, 1662, 1790, 1921, 2066, 2282, 2469, 2666, 2855, 3129, 3424, 3772, 4012, 4242, 4449, 4631, 4892, 5129, 5358, 5497, 5616, 5685, 5761, 5824, 5889, 5973, 6051, 6137, 6210, 6258, 6329, 6392, 6453, 6524, 6578, 6624, 6663, 6708, 6736, 6766, 6793, 6826, 6856, 6868, 6887, 6914, 6945, 6977, 7001, 7010)
(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 17, 21, 27, 34, 54, 75, 96, 121, 144, 184, 211, 247, 292, 351, 417, 482, 552, 635, 708, 771, 825, 871, 920, 993, 1047, 1101, 1146, 1174, 1210, 1245, 1280, 1328, 1370, 1438, 1511, 1622, 1758, 1908, 2114, 2311, 2485, 2664, 2897, 3172, 3507, 3758, 3967, 4154, 4337, 4609, 4819, 5034, 5168, 5285, 5345, 5421, 5487, 5554, 5609, 5687, 5780, 5848, 5902, 5956, 6024, 6083, 6145, 6205, 6247, 6283)

The last dozen of so data points will change tomorrow.


graphs

I’ve lobbed off the first four data points to remove the starting spike that messes with the graph. Since this is a daily change then spikes are a norm. Taking out the big one at the start makes the graph easier to interpret going forward.

Day-on-Day Rate of Change

Cases


code

Latest update, May 20, 2020 (added Recoveries graph)
May 23, 2020 (implemented Decimal class; added input parameter to main())

import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from decimal import Decimal, Context, setcontext
from data import data

def case_delta(u):
    '''
    ~ Extract serial data objects required by the graphs ~
    The data is imported as a string object 
    with columnal data, Month, Day, Cases, Recoveries, respectively. 
    > u: the string object `data` in 'data.py' companion file
    > s: row data minus first and last line of string
    > r: column 3 data - recoveries
    > t: column 2 data - cases to date
    >>>: list of day[n] / day[n-1], t, and r
    '''
    context = Context(prec=4)
    setcontext(context)
    s = tuple(map(lambda x: x.split(), u.split('\n')[1:-1]))
    r = tuple(map(lambda x: int(x[3]), s))
    t = tuple(map(lambda x: int(x[2]), s))
    return [f"{(Decimal(b) / Decimal(a))}" for (a, b) in zip(t[:-1], t[1:])], t, r

def roc_plot(sample):
    y = [*map(lambda x: float(x), sample[4:])]
    m = datetime.now()
    n = len(y)
    o = m - timedelta(days=n)
    x = range(n)
    t = [1.2] * len(x)
    u = [1.1] * len(x)
    v = [1.0] * len(x)
    plt.title(f'''
Alberta Reported Cases Day-on-Day Rate of Change
{str(o)[:10]} - {str(m)[:10]}
    ''')
    plt.plot(x, t, 'r')
    plt.plot(x, u, 'y')
    plt.plot(x, v, 'g')
    plt.plot(x, y)
    plt.legend(['20% Change', '10% Change', ' 0% Change', 'Plots'])
    plt.xlabel('Days Running')
    plt.ylabel('Day to Day Ratios')
    plt.show()

def cases_plot(cases, recoveries):
    y, z = cases[5:], recoveries[5:]
    n = len(y)
    x = range(n)
    m = datetime.now()
    o = m - timedelta(days=n)
    plt.title(f'''
Alberta Reported Cases
{str(o)[:10]} - {str(m)[:10]}
    ''')
    plt.plot(x, y)
    plt.plot(x, z, 'g')
    plt.legend([
        f'{y[-1]} Cases', 
        f'{z[-1]} Recovered'
    ])
    plt.xlabel('Days Running')
    plt.ylabel('Cases to Date')
    plt.show()

def main(n):
    ratios, cases, recoveries = case_delta(n)
    print (ratios)
    print (cases)
    print (recoveries)
    roc_plot(ratios)
    cases_plot(cases, recoveries)

main(data)