import pandas as pd
from datetime import datetime, timedelta
import statsmodels.api as sm
from IPython.core.display import display, HTML
import os
from src import *
pd.set_option('display.max_rows', None)
These parameters configure various things such as the relative location of data files, which date is being analysed and the assumed generation period for an infection.
ALL_VIC_CASES='archive/2021-07-12/all-vic-cases.csv'
QUARANTINE='archive/2021-07-16/quarantine.csv'
GENERATION_DAYS=5
ENABLE_MELBOURNE_ANIMATION=False
TODAY=datetime.today().strftime("%Y-%m-%d")
RUN_DATE=TODAY
#RUN_DATE='2022-02-20'
YESTERDAY=add_days(RUN_DATE,-1)
TOMORROW=add_days(RUN_DATE, 1)
LAST_WEEK=add_days(RUN_DATE, -7)
#RUN_DATE=YESTERDAY
LOAD_DATE=RUN_DATE
PREV_DAY=add_days(RUN_DATE, -1)
NEXT_WEEK=add_days(RUN_DATE, 7)
PREV_WEEK=add_days(RUN_DATE, -7)
PREV_FORTNIGHT=add_days(RUN_DATE, -14)
NSW_6_MONTHS=lambda date: f"archive/{date}/last-6-months-nsw.csv"
NSW_14_DAYS=lambda date: f"archive/{date}/last-14-days-nsw.csv"
os.makedirs(f"archive/{RUN_DATE}", exist_ok=True)
AMEND_TOTAL=None
date=RUN_DATE
index=load_index()
try:
parse_statistics_html(load_statistics(index[date], date))
except KeyError as e:
pass
if False:
try:
with open(f'archive/{RUN_DATE}/statistics.json') as f:
pass
except FileNotFoundError as e:
with open(f'archive/{YESTERDAY}/statistics.json') as f:
with open(f'archive/{RUN_DATE}/statistics.json', 'w') as wf:
wf.write(f.read())
Currently, we parse tne contents of https://www.health.nsw.gov.au/news/Pages/2022-nsw-health.aspx and then for each date with a statistics HTML we populate a per-date cache by downloading the page NSW health if required or using a previously downloaded page if it is exists. We then parse each of these to obtain the "total" and "cumulative_corrected" statistics for each date. We calculate "cumulative" as the cumulative sum of total and "correction" as the difference between "cumulative_corrected" and "cumulative".
#nsw_df = update_df(load_data(NSW_6_MONTHS(LOAD_DATE)),load_data(NSW_14_DAYS(LOAD_DATE)))
#nsw_df = load_data(NSW_6_MONTHS(LOAD_DATE))
nsw_df = load_nswhealth_stats(limit_date=RUN_DATE)
if AMEND_TOTAL:
if len(nsw_df[nsw_df["date"]==RUN_DATE])==0:
nsw_df = append_recent_date(nsw_df, RUN_DATE, AMEND_TOTAL)
vic_df = load_vic_data(ALL_VIC_CASES)
quarantine_df = load_quarantine(QUARANTINE)
merge_rats(nsw_df)
Next, we truncate and refindex the data frames to each outbreak.
#avalon = select_outbreak(nsw_df[(nsw_df['date'] >= '2020-12-17') & (nsw_df['date'] <= '2021-01-16')], generation_days=GENERATION_DAYS)
bondi = select_outbreak(nsw_df[(nsw_df['date'] >= '2021-06-17')], generation_days=GENERATION_DAYS)
vic_outbreak = select_outbreak(vic_df[(vic_df['date'] >= '2020-05-27') & (vic_df['date'] <= '2020-10-29')], generation_days=GENERATION_DAYS)
/usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs)
#merge_rats(bondi)
#merge_rats_alt(bondi)
bondi[
['date', 'pcr', 'rat', 'total', 'rat_added', 'rat_added_last_week']
][bondi.date>'2021-12-31'].to_csv(f'archive/{RUN_DATE}/rat-redistribution.csv', index=False)
slice=bondi.loc[bondi['date'] > '2021-12-31', ['date', 'rat','pcr', 'total']]
ax=slice.plot.bar(x='date', figsize=(10,10))
ax.set_title('(A) redistribution of NSW RAT tests since January 1, 2022')
#ax.set_yscale('log')
ax.figure.savefig(f'archive/{RUN_DATE}/rat-redistribution.png')
_=_
bondi[['rat', 'rat_added']].sum()
rat 435718 rat_added 421995 dtype: int64
if False:
ax=bondi[bondi.date > '2021-12-31'][['date','rat_ratio']].plot(x='date', figsize=(10,10))
_=ax.set_title("RAT:PCR ratio since Jan 1, 2022")
ax.figure.savefig(f'archive/{RUN_DATE}/rat-ratio.png')
if False:
bondi[['date', 'rat_ratio', 'rat_ratio_pre_last_week', 'rat_ratio_last_week']][bondi.rat_ratio_last_week.notna()]
display(HTML(summary(bondi)))
Date: 2022-03-09 (#265) New Cases Reported Today: 14775 ↓Δ -116.0 Cumulative (since 2021-06-16): 1432855 ↑Δ 14775.0 Corrections not available on this day. Projection (from yesterday): 10560 Projection Error: -4215 (-28.5%) Projection (for tomorrow): 12158 Projection (for next week): 10116 Cumulative Growth Rate: 0.85% per day ↑Δ 0.104 Linear Growth Rate : 1.03% per day ↓Δ -0.019 Reff: 1.11 ↑Δ 0.074 Doubling Period: 82.0 days ↓Δ -11.388 Growth Decay Rate: -1.18% per day ↑Δ 0.01 ln-ln Gradient: 7.344 ↑Δ 1.655
display(HTML(weekly_summary(bondi, 0.0)))
Delta values are with respect to values from 1 week prior.
Date: 2022-03-09 (#265) New Cases Reported Today: 14775 ↑Δ 3804.0 Cumulative (since 2021-06-16): 1432855 ↑Δ 78911.0 Corrections not available on this day. Projection (from 7 days ago): 6959 Projection Error: -7816 (-52.9%) Projection (for next week): 10116 Care Factor: 0.0 Cumulative Growth Rate: 0.85% per day ↑Δ 0.218 Linear Growth Rate : 1.03% per day ↑Δ 0.221 Reff: 1.11 ↑Δ 0.11 Doubling Period: 82.0 days ↓Δ -28.307 Growth Decay Rate: -1.18% per day ↑Δ 0.17 ln-ln Gradient: 7.344 ↑Δ 4.167
If the model was perfect, tbe following would be true:
$ \space\space\space\space\space\text{cumulative}_{t+1} = (1+\frac{\text{g}_t}{100})\times\text{cumulative}_t $
In practice, the growth rate, $g_t$, is estimated by fitting a linear regession through a set of points:
$ \space\space\space\space\space(t-k, \ln{(\text{cumulative}_{t-k})}) $
for k = 0..4. This yields parameters of a linear reqregssion:
$ \space\space\space\space\space\ln(\text{cumulative}_t) = m \times t + b $
Raisng e to each side, yields:
$ \space\space\space\space\space\text{cumulative}_t = (e^m)^t \times e^b = e^m \times (e^m)^{t-1} \times e^b = e^m \times \text{cumulative}_{t-1} $
The growth rate, $g_t$, is thus:
$ \space\space\space\space\space g_t = (e^m-1) * 100 $
The minimum growth rate $g_{t,min}$ is defined as:
$ \space\space\space\space\space g_{t,min} = min(\{k: 0 \le k < 5: g_{t-k}\}) $
The 7-day forward projection of the cumulative total, $projection_t$ used on the so-called "hedgehog plot" uses $g_{t,min}$ since this was observed to provide closest fit to the Melbourne 2020 outbreak, particularly in the later stages.
$ \space\space\space\space\space \text{projection}_t = (1+\frac{g_{t-7,min}}{100})^7 \times \text{cumulative}_{t-7} $
The intuitive justification is that because the growth rate eventually starts to decay, a fit to the last 5 days growth is going to tend overestimate the growth in the next 5 days so choosing the minimum observed growth rate estimate is more likely to be closer to true growth rate rather than the very latest estimate.
The replication factor, $R_{eff}$ is defined as:
$ \space\space\space\space\space R_{eff} = (1+\frac{g_t}{100})^5 $
The doubling period, in days, is defined as:
$ \space\space\space\space\space \text{doubling period}_t = \ln{(2)}/\ln{(1+\frac{g_t}{100})} $
The table below documents various calculated parameters.
display(bondi[[
"date",
"cumulative",
"total",
'rolling-average',
"ols-growth-rate",
"ols-growth-rate-min",
"ols-growth-rate-decay",
"linear-growth-rate",
"doubling-period",
"Reff",
"ltlc-gradient"
]].tail(56))
date | cumulative | total | rolling-average | ols-growth-rate | ols-growth-rate-min | ols-growth-rate-decay | linear-growth-rate | doubling-period | Reff | ltlc-gradient | |
---|---|---|---|---|---|---|---|---|---|---|---|
210 | 2022-01-13 | 710182 | 51578 | 48549.000000 | 7.707883 | 7.592022 | -13.522447 | 7.262645 | 9.334989 | 0.917548 | 0.646545 |
211 | 2022-01-14 | 748666 | 38484 | 46520.285714 | 7.689807 | 7.592022 | -8.475658 | 5.140343 | 9.356138 | 0.917992 | 0.523695 |
212 | 2022-01-15 | 788703 | 40037 | 43164.857143 | 6.912675 | 6.912675 | -4.718554 | 5.076309 | 10.369904 | 0.943774 | 0.460278 |
213 | 2022-01-16 | 819376 | 30673 | 41249.000000 | 5.566509 | 5.566509 | -7.035703 | 3.743458 | 12.795544 | 0.930177 | 0.222124 |
214 | 2022-01-17 | 845060 | 25684 | 40517.714286 | 4.477576 | 4.477576 | -13.146301 | 3.039311 | 15.824453 | 0.832494 | -0.082787 |
215 | 2022-01-18 | 868509 | 23449 | 38072.142857 | 3.727676 | 3.727676 | -17.159211 | 2.699914 | 18.939079 | 0.657738 | -0.473835 |
216 | 2022-01-19 | 897838 | 29329 | 34176.285714 | 3.225265 | 3.225265 | -17.516064 | 3.266625 | 21.835910 | 0.560179 | -0.662448 |
217 | 2022-01-20 | 924200 | 26362 | 30574.000000 | 3.059395 | 3.059395 | -14.145557 | 2.852413 | 23.001181 | 0.557547 | -0.939613 |
218 | 2022-01-21 | 947062 | 22862 | 28342.285714 | 2.943123 | 2.943123 | -9.848531 | 2.413992 | 23.896319 | 0.633511 | -1.128480 |
219 | 2022-01-22 | 964919 | 17857 | 25173.714286 | 2.674205 | 2.674205 | -7.279533 | 1.850622 | 26.264794 | 0.700396 | -1.185136 |
220 | 2022-01-23 | 986121 | 21202 | 23820.714286 | 2.333761 | 2.333761 | -7.518663 | 2.150040 | 30.046109 | 0.741540 | -1.369569 |
221 | 2022-01-24 | 999212 | 13091 | 22021.714286 | 1.984346 | 1.984346 | -10.397286 | 1.310132 | 35.276202 | 0.718044 | -2.277414 |
222 | 2022-01-25 | 1014558 | 15346 | 20864.142857 | 1.741087 | 1.741087 | -12.613433 | 1.512580 | 40.156758 | 0.650806 | -3.009679 |
223 | 2022-01-26 | 1036283 | 21725 | 19777.857143 | 1.726047 | 1.726047 | -11.029095 | 2.096435 | 40.503646 | 0.652460 | -2.852606 |
224 | 2022-01-27 | 1053428 | 17145 | 18461.142857 | 1.699075 | 1.699075 | -7.450349 | 1.627544 | 41.141159 | 0.706002 | -2.649100 |
225 | 2022-01-28 | 1064649 | 11221 | 16798.142857 | 1.658234 | 1.658234 | -3.762377 | 1.053962 | 42.145954 | 0.766670 | -3.063376 |
226 | 2022-01-29 | 1079211 | 14562 | 16327.428571 | 1.516983 | 1.516983 | -3.107226 | 1.349319 | 46.038172 | 0.830609 | -2.839699 |
227 | 2022-01-30 | 1093411 | 14200 | 15327.142857 | 1.323729 | 1.323729 | -6.238169 | 1.298688 | 52.709048 | 0.840233 | -2.841065 |
228 | 2022-01-31 | 1107576 | 14165 | 15480.571429 | 1.277141 | 1.277141 | -7.653465 | 1.278919 | 54.619172 | 0.807154 | -3.011739 |
229 | 2022-02-01 | 1117188 | 9612 | 14661.428571 | 1.230336 | 1.230336 | -7.402037 | 0.860374 | 56.683893 | 0.766487 | -3.865704 |
230 | 2022-02-02 | 1128495 | 11307 | 13173.142857 | 1.114383 | 1.114383 | -6.667428 | 1.001954 | 62.546001 | 0.717652 | -3.705237 |
231 | 2022-02-03 | 1140275 | 11780 | 12406.714286 | 1.031742 | 1.031742 | -6.150055 | 1.033084 | 67.528172 | 0.732842 | -3.339800 |
232 | 2022-02-04 | 1150389 | 10114 | 12248.571429 | 0.967725 | 0.967725 | -7.048299 | 0.879181 | 71.972480 | 0.782130 | -3.272867 |
233 | 2022-02-05 | 1157423 | 7034 | 11173.142857 | 0.903836 | 0.903836 | -7.298896 | 0.607729 | 77.035530 | 0.738546 | -4.352358 |
234 | 2022-02-06 | 1165711 | 8288 | 10328.571429 | 0.801386 | 0.801386 | -7.612499 | 0.710982 | 86.839710 | 0.686659 | -4.347215 |
235 | 2022-02-07 | 1174093 | 8382 | 9502.428571 | 0.719415 | 0.719415 | -8.695822 | 0.713913 | 96.694822 | 0.689098 | -5.689141 |
236 | 2022-02-08 | 1180830 | 6737 | 9091.714286 | 0.667567 | 0.667567 | -9.251901 | 0.570531 | 104.178106 | 0.680755 | -7.124020 |
237 | 2022-02-09 | 1191221 | 10391 | 8960.857143 | 0.707008 | 0.667567 | -6.517077 | 0.872298 | 98.385625 | 0.748876 | -5.584670 |
238 | 2022-02-10 | 1201756 | 10535 | 8783.000000 | 0.756730 | 0.667567 | -1.311990 | 0.876634 | 91.943802 | 0.864583 | -4.198438 |
239 | 2022-02-11 | 1210810 | 9054 | 8631.571429 | 0.794676 | 0.667567 | 3.296757 | 0.747764 | 87.570015 | 1.034100 | -4.711161 |
240 | 2022-02-12 | 1218768 | 7958 | 8763.571429 | 0.798739 | 0.667567 | 4.871769 | 0.652954 | 87.126308 | 1.157719 | -4.289469 |
241 | 2022-02-13 | 1224549 | 5781 | 8405.428571 | 0.694845 | 0.694845 | 0.193400 | 0.472092 | 100.101760 | 1.074826 | -4.643874 |
242 | 2022-02-14 | 1233897 | 9348 | 8543.428571 | 0.642760 | 0.642760 | -4.502633 | 0.757600 | 108.185375 | 0.947782 | -2.718633 |
243 | 2022-02-15 | 1237222 | 3325 | 8056.000000 | 0.556492 | 0.556492 | -8.879064 | 0.268747 | 124.902829 | 0.768608 | -6.067327 |
244 | 2022-02-16 | 1245685 | 8463 | 7780.571429 | 0.541321 | 0.541321 | -9.517001 | 0.679385 | 128.393722 | 0.703106 | -4.722641 |
245 | 2022-02-17 | 1253790 | 8105 | 7433.428571 | 0.568659 | 0.541321 | -5.564931 | 0.646440 | 122.237771 | 0.724162 | -2.949518 |
246 | 2022-02-18 | 1263209 | 9419 | 7485.571429 | 0.604400 | 0.541321 | -1.009278 | 0.745641 | 115.029652 | 0.900207 | -0.952552 |
247 | 2022-02-19 | 1268004 | 4795 | 7033.714286 | 0.633203 | 0.541321 | 3.753763 | 0.378153 | 109.813039 | 1.012524 | -3.529996 |
248 | 2022-02-20 | 1274061 | 6057 | 7073.142857 | 0.564797 | 0.541321 | 1.942813 | 0.475409 | 123.071345 | 1.042425 | -4.190637 |
249 | 2022-02-21 | 1279788 | 5727 | 6555.857143 | 0.497243 | 0.497243 | -3.305880 | 0.447496 | 139.744377 | 0.950549 | -4.887105 |
250 | 2022-02-22 | 1289261 | 9473 | 7434.142857 | 0.502038 | 0.497243 | -5.944347 | 0.734762 | 138.413005 | 0.917458 | -4.393923 |
251 | 2022-02-23 | 1298010 | 8749 | 7475.000000 | 0.588086 | 0.497243 | -2.621282 | 0.674032 | 118.211128 | 1.001223 | -1.557331 |
252 | 2022-02-24 | 1306125 | 8115 | 7476.428571 | 0.640528 | 0.497243 | 4.283687 | 0.621303 | 108.561255 | 1.045336 | 1.184808 |
253 | 2022-02-25 | 1313521 | 7396 | 7187.428571 | 0.652412 | 0.497243 | 8.185936 | 0.563067 | 106.589965 | 1.199527 | 2.686031 |
254 | 2022-02-26 | 1320327 | 6806 | 7474.714286 | 0.596769 | 0.502038 | 4.597665 | 0.515478 | 116.496323 | 1.216977 | 3.149438 |
255 | 2022-02-27 | 1326248 | 5921 | 7455.285714 | 0.540032 | 0.540032 | -2.383680 | 0.446447 | 128.699364 | 1.096003 | 1.300156 |
256 | 2022-02-28 | 1332591 | 6343 | 7543.285714 | 0.498875 | 0.498875 | -6.657282 | 0.475990 | 139.288429 | 0.878252 | 2.440375 |
257 | 2022-03-01 | 1342973 | 10382 | 7673.142857 | 0.537386 | 0.498875 | -5.513137 | 0.773061 | 129.331232 | 0.862460 | 0.308489 |
258 | 2022-03-02 | 1353944 | 10971 | 7990.571429 | 0.630144 | 0.498875 | 1.044688 | 0.810299 | 110.344459 | 0.995263 | 3.176799 |
259 | 2022-03-03 | 1365808 | 11864 | 8526.142857 | 0.749608 | 0.498875 | 9.301931 | 0.868643 | 92.814107 | 1.329463 | 5.819960 |
260 | 2022-03-04 | 1375172 | 9364 | 8807.285714 | 0.800867 | 0.498875 | 13.649762 | 0.680933 | 86.895748 | 1.571783 | 7.321348 |
261 | 2022-03-05 | 1385422 | 10250 | 9299.285714 | 0.780983 | 0.537386 | 10.378156 | 0.739847 | 89.099256 | 1.577327 | 5.947243 |
262 | 2022-03-06 | 1393949 | 8527 | 9671.571429 | 0.727598 | 0.630144 | 3.340645 | 0.611715 | 95.611236 | 1.359148 | 4.483346 |
263 | 2022-03-07 | 1403189 | 9240 | 10085.428571 | 0.677933 | 0.677933 | -2.925835 | 0.658500 | 102.590396 | 1.105233 | 2.979471 |
264 | 2022-03-08 | 1418080 | 14891 | 10729.571429 | 0.744687 | 0.677933 | -2.828884 | 1.050082 | 93.425157 | 1.031535 | 5.688811 |
265 | 2022-03-09 | 1432855 | 14775 | 11273.000000 | 0.848494 | 0.677933 | 1.908311 | 1.031158 | 82.037543 | 1.105632 | 7.343633 |
This section presents the 7-day old, 7-day forward projection, e.g. the prediction for the current date and a comparison with what actually happened. This is about measuring the effectiveness for a 7-day model for a known period. In particular, we are not trying to predict what will happen in the next 7 days/
display(bondi[[
"date",
"cumulative",
"correction",
"total",
"7-day-delta",
"7-day-projection",
"7-day-projection-error",
"7-day-projection-relative-error",
"7-day-forward-projection-cumulative",
"7-day-forward-projection-total"
]].tail(28))
date | cumulative | correction | total | 7-day-delta | 7-day-projection | 7-day-projection-error | 7-day-projection-relative-error | 7-day-forward-projection-cumulative | 7-day-forward-projection-total | |
---|---|---|---|---|---|---|---|---|---|---|
238 | 2022-02-10 | 1201756 | NaN | 10535 | 61481.0 | 1225221.0 | 23465.0 | 38.166263 | 1259051 | 8349 |
239 | 2022-02-11 | 1210810 | NaN | 9054 | 60421.0 | 1230616.0 | 19806.0 | 32.779994 | 1268537 | 8412 |
240 | 2022-02-12 | 1218768 | NaN | 7958 | 61345.0 | 1232667.0 | 13899.0 | 22.657103 | 1276874 | 8467 |
241 | 2022-02-13 | 1224549 | NaN | 5781 | 58838.0 | 1232697.0 | 8148.0 | 13.848193 | 1285366 | 8870 |
242 | 2022-02-14 | 1233897 | NaN | 9348 | 59804.0 | 1234511.0 | 614.0 | 1.026687 | 1290496 | 8242 |
243 | 2022-02-15 | 1237222 | NaN | 3325 | 56392.0 | 1237127.0 | -95.0 | -0.168464 | 1286229 | 7118 |
244 | 2022-02-16 | 1245685 | NaN | 8463 | 54464.0 | 1248014.0 | 2329.0 | 4.276219 | 1293661 | 6965 |
245 | 2022-02-17 | 1253790 | NaN | 8105 | 52034.0 | 1259051.0 | 5261.0 | 10.110697 | 1302078 | 7010 |
246 | 2022-02-18 | 1263209 | NaN | 9419 | 52399.0 | 1268537.0 | 5328.0 | 10.168133 | 1311859 | 7063 |
247 | 2022-02-19 | 1268004 | NaN | 4795 | 49236.0 | 1276874.0 | 8870.0 | 18.015273 | 1316839 | 7090 |
248 | 2022-02-20 | 1274061 | NaN | 6057 | 49512.0 | 1285366.0 | 11305.0 | 22.832849 | 1323129 | 7124 |
249 | 2022-02-21 | 1279788 | NaN | 5727 | 45891.0 | 1290496.0 | 10708.0 | 23.333551 | 1325004 | 6556 |
250 | 2022-02-22 | 1289261 | NaN | 9473 | 52039.0 | 1286229.0 | -3032.0 | -5.826399 | 1334811 | 6604 |
251 | 2022-02-23 | 1298010 | NaN | 8749 | 52325.0 | 1293661.0 | -4349.0 | -8.311515 | 1343869 | 6649 |
252 | 2022-02-24 | 1306125 | NaN | 8115 | 52335.0 | 1302078.0 | -4047.0 | -7.732875 | 1352271 | 6691 |
253 | 2022-02-25 | 1313521 | NaN | 7396 | 50312.0 | 1311859.0 | -1662.0 | -3.303387 | 1359928 | 6729 |
254 | 2022-02-26 | 1320327 | NaN | 6806 | 52323.0 | 1316839.0 | -3488.0 | -6.666284 | 1367431 | 6831 |
255 | 2022-02-27 | 1326248 | NaN | 5921 | 52187.0 | 1323129.0 | -3119.0 | -5.976584 | 1377203 | 7397 |
256 | 2022-02-28 | 1332591 | NaN | 6343 | 52803.0 | 1325004.0 | -7587.0 | -14.368502 | 1379829 | 6849 |
257 | 2022-03-01 | 1342973 | NaN | 10382 | 53712.0 | 1334811.0 | -8162.0 | -15.195859 | 1390579 | 6903 |
258 | 2022-03-02 | 1353944 | NaN | 10971 | 55934.0 | 1343869.0 | -10075.0 | -18.012300 | 1401939 | 6959 |
259 | 2022-03-03 | 1365808 | NaN | 11864 | 59683.0 | 1352271.0 | -13537.0 | -22.681501 | 1414223 | 7020 |
260 | 2022-03-04 | 1375172 | NaN | 9364 | 61651.0 | 1359928.0 | -15244.0 | -24.726282 | 1423919 | 7068 |
261 | 2022-03-05 | 1385422 | NaN | 10250 | 65095.0 | 1367431.0 | -17991.0 | -27.638067 | 1438385 | 7688 |
262 | 2022-03-06 | 1393949 | NaN | 8527 | 67701.0 | 1377203.0 | -16746.0 | -24.735233 | 1456611 | 9121 |
263 | 2022-03-07 | 1403189 | NaN | 9240 | 70598.0 | 1379829.0 | -23360.0 | -33.088756 | 1471147 | 9906 |
264 | 2022-03-08 | 1418080 | NaN | 14891 | 75107.0 | 1390579.0 | -27501.0 | -36.615762 | 1486760 | 10011 |
265 | 2022-03-09 | 1432855 | NaN | 14775 | 78911.0 | 1401939.0 | -30916.0 | -39.178315 | 1502250 | 10116 |
display(bondi[[
"date",
"cumulative",
"total",
"one-day-projection-cumulative",
"one-day-projection-total",
"one-day-error",
"one-day-relative-error",
]].tail(28))
date | cumulative | total | one-day-projection-cumulative | one-day-projection-total | one-day-error | one-day-relative-error | |
---|---|---|---|---|---|---|---|
238 | 2022-02-10 | 1201756 | 10535 | 1210850.0 | 9094.0 | -2113.0 | -20.056953 |
239 | 2022-02-11 | 1210810 | 9054 | 1220432.0 | 9622.0 | 40.0 | 0.441794 |
240 | 2022-02-12 | 1218768 | 7958 | 1228503.0 | 9735.0 | 1664.0 | 20.909776 |
241 | 2022-02-13 | 1224549 | 5781 | 1233058.0 | 8509.0 | 3954.0 | 68.396471 |
242 | 2022-02-14 | 1233897 | 9348 | 1241828.0 | 7931.0 | -839.0 | -8.975182 |
243 | 2022-02-15 | 1237222 | 3325 | 1244107.0 | 6885.0 | 4606.0 | 138.526316 |
244 | 2022-02-16 | 1245685 | 8463 | 1252428.0 | 6743.0 | -1578.0 | -18.645870 |
245 | 2022-02-17 | 1253790 | 8105 | 1260920.0 | 7130.0 | -1362.0 | -16.804442 |
246 | 2022-02-18 | 1263209 | 9419 | 1270844.0 | 7635.0 | -2289.0 | -24.301943 |
247 | 2022-02-19 | 1268004 | 4795 | 1276033.0 | 8029.0 | 2840.0 | 59.228363 |
248 | 2022-02-20 | 1274061 | 6057 | 1281257.0 | 7196.0 | 1972.0 | 32.557372 |
249 | 2022-02-21 | 1279788 | 5727 | 1286152.0 | 6364.0 | 1469.0 | 25.650428 |
250 | 2022-02-22 | 1289261 | 9473 | 1295734.0 | 6473.0 | -3109.0 | -32.819593 |
251 | 2022-02-23 | 1298010 | 8749 | 1305643.0 | 7633.0 | -2276.0 | -26.014402 |
252 | 2022-02-24 | 1306125 | 8115 | 1314491.0 | 8366.0 | -482.0 | -5.939618 |
253 | 2022-02-25 | 1313521 | 7396 | 1322091.0 | 8570.0 | 970.0 | 13.115197 |
254 | 2022-02-26 | 1320327 | 6806 | 1328206.0 | 7879.0 | 1764.0 | 25.918307 |
255 | 2022-02-27 | 1326248 | 5921 | 1333410.0 | 7162.0 | 1958.0 | 33.068738 |
256 | 2022-02-28 | 1332591 | 6343 | 1339239.0 | 6648.0 | 819.0 | 12.911871 |
257 | 2022-03-01 | 1342973 | 10382 | 1350190.0 | 7217.0 | -3734.0 | -35.966095 |
258 | 2022-03-02 | 1353944 | 10971 | 1362476.0 | 8532.0 | -3754.0 | -34.217482 |
259 | 2022-03-03 | 1365808 | 11864 | 1376046.0 | 10238.0 | -3332.0 | -28.084963 |
260 | 2022-03-04 | 1375172 | 9364 | 1386185.0 | 11013.0 | 874.0 | 9.333618 |
261 | 2022-03-05 | 1385422 | 10250 | 1396242.0 | 10820.0 | 763.0 | 7.443902 |
262 | 2022-03-06 | 1393949 | 8527 | 1404091.0 | 10142.0 | 2293.0 | 26.891052 |
263 | 2022-03-07 | 1403189 | 9240 | 1412702.0 | 9513.0 | 902.0 | 9.761905 |
264 | 2022-03-08 | 1418080 | 14891 | 1428640.0 | 10560.0 | -5378.0 | -36.115775 |
265 | 2022-03-09 | 1432855 | 14775 | 1445013.0 | 12158.0 | -4215.0 | -28.527919 |
output=pd.DataFrame(columns=['cumulative', 'min', 'vic'])
output[['min', 'cumulative']] = bondi[['min','cumulative']]
output=output.reindex([r for r in range(0, len(bondi)+14)])
output['vic'] = vic_outbreak['cumulative']
x=11
output['vic-offset'] = vic_outbreak['cumulative'].shift(-x)
ax=output.plot(figsize=(10,10))
ax.set_yscale('log')
ax.grid()
ax.set_title(f'Cumulative Case Growth Projections ({RUN_DATE})')
ax.legend([
'Sydney (2021) - actual',
'Sydney (2021) - model',
'Melbourne (2020)',
f'Melbourne (2020) +{x} days'
])
ax.figure.savefig(f'archive/{RUN_DATE}/cumulative-partial.png')
_=_
df=vic_outbreak
ax=df[['min', 'cumulative']].plot(figsize=(10,10))
#ax.set_yscale('log')
ax.grid()
ax.plot(bondi['cumulative'])
ax.plot(bondi['min'])
ax.legend(['model (Melbourne 2020) ', 'cumulative (Melbourne 2020)','cumulative (Sydney 2021)', 'model (Sydney 2021) '])
ax.set_title("7 Day Projection vs Actual (Melbourne 2020, Sydney 2021)")
ax.figure.savefig(f'archive/{RUN_DATE}/cumulative-full.png')
_=_
VIC_EXTRA_DAYS=0
vic_growth_params=derive_growth_params(vic_outbreak[(vic_outbreak.index >= 70) & (vic_outbreak.index < 120)], generation_days=GENERATION_DAYS)
bondi_growth_params=derive_growth_params(bondi.tail(8), generation_days=GENERATION_DAYS)
N=1
bondi_growth_params_3=derive_growth_params(bondi.tail(8+N).head(8), generation_days=GENERATION_DAYS)
gp=derive_growth_params(bondi[(bondi.index>15)])
bondi_projection_1=select_outbreak(project_ols_growth_rate_min(bondi, 84, vic_growth_params[1]))
bondi_projection_2=select_outbreak(project_ols_growth_rate_min(bondi, 200, gp[1]))
bondi_projection_3=select_outbreak(project_ols_growth_rate_min(bondi.head(len(bondi)-N), 84+N, gp[1]))
vic_partial=vic_outbreak.head(len(bondi)+VIC_EXTRA_DAYS)
vic_partial_growth_params=derive_growth_params(vic_partial)
vic_projection=select_outbreak(project_ols_growth_rate_min(vic_partial, len(vic_outbreak)-len(bondi)-VIC_EXTRA_DAYS, vic_partial_growth_params[1]), generation_days=GENERATION_DAYS)
/usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs)
This plot shows the estimated and projected daily cumulative growth rate as a percentage. By definition, the cumulative growth rate is always an (eventually small) non-negative value. The values plotted for each outbreak are $g_t$ and $g'_t$ as calculated above.
gp=GrowthPlot(RUN_DATE)
gp.add(bondi, offset=0, legend="Sydney 2021")
#gp.add(avalon, offset=0, legend="Avalon 2020")
gp.add(vic_outbreak, offset=0, legend="Melbourne 2020")
gp.ax.plot(bondi['ols-growth-rate-min'],color="C3")
gp.ax.plot(vic_outbreak['ols-growth-rate-min'], color="C4")
gp.ax.plot(bondi_projection_1['ols-growth-rate-min'], linestyle='dotted', color='C4')
gp.ax.plot(bondi_projection_2['ols-growth-rate-min'], linestyle='dashed', color='C3')
#gp.ax.plot(bondi_projection_3['ols-growth-rate-min'], linestyle='dotted', color='C5')
gp.ax.set_yscale('log')
gp.legend = gp.legend+[
'Sydney 2021 (retrospective model)',
'Melbourne 2020 (retrospective model)',
'Sydney 2021 (projection - Melbourne 2020 decay)',
f'Sydney 2021 (projection - Sydney 2021 decay ({RUN_DATE})',
# f'Sydney 2021 (projection - Sydney 2021 decay ({PREV_DAY})'
]
gp.ax.legend(gp.legend)
gp.ax.figure.savefig(f'archive/{RUN_DATE}/cumulative-growth.png')
#gp.add(vic_outbreak.shift(-11), offset=0, legend="Melbourne 2020 (shifted)")
_=_
This plot display the relative error between the projected 7-day cumulative growth and the growth that occurred in practice.
The relative error is defined as:
$ \space\space\space\text{error}_t = \text{projection}_t - \text{cumulative}_t $
$ \space\space\space\text{relative error}_t = \frac{\text{projection}_t - \text{cumulative}_t}{\text{cumulative}_t - \text{cumulative}_{t-7}} \times 100 $
output=pd.DataFrame()
#output['vic'] = modeling_errors(vic_outbreak)
output["bondi"] = modeling_errors(bondi)
#output["avalon"] = modeling_errors(avalon)
ax=output.loc[output.index >= 15, ['bondi', ]].plot(figsize=(10,10))
ax.grid()
ax.set_title("Growth Projection Relative Error % vs Day Of Outbreak")
ax.set_xlabel("Day Of Outbreak")
ax.set_ylabel("% overshoot of projection vs actual")
ax.legend(['Sydney 2011',])# 'Melbourne 2020'])
ax.figure.savefig(f'archive/{RUN_DATE}/modellng-error.png')
_=_
df=bondi
ax=df["Reff"].plot(figsize=(10,10))
ax.set_title("Reff vs Day Of Outbreak (Sydney 2021)")
ax.set_yticks([r/2 for r in range(0,16)])
ax.hlines(1,0,df.index.max(), linestyle="dashed", color="C1")
ax.vlines(df[df['date']=="2021-10-11"].index, 0, 7, linestyle="dashed", color="C2")
ax.legend(["Reff", "Reff == 1", "Restrictions Relaxed (2021-10-11)"])
ax.figure.savefig(f"archive/{RUN_DATE}/reff.png")
plot_vic=False
GRAPH_DATE=bondi.tail(1).date.values[0]
pp = PhasePlot(f"New Cases vs Error Modeling % - Sydney 2021 vs Melbourne 2020 - ({GRAPH_DATE})")
bondi_idx=pp.add(bondi, offset=20, legend="Sydney 2021", color="C0") # 20
if plot_vic:
vic_idx=pp.add(vic_outbreak, offset=15, legend="Melbourne 2020", color="C1") # 15
pp.add_horizon(horizon(bondi, 7), legend=f"7-Day Projection for {NEXT_WEEK}", color="blue")
pp.add_horizon(horizon(bondi.head(len(bondi)-7), 7), legend=f"7-Day Projection for {RUN_DATE}", color="red")
# for i in range(0,7):
# pp.add_horizon(horizon(bondi.head(len(bondi)-7-i),7), legend=f"today - as projected {7+i} days ago", color=f"C{i}")
if plot_vic:
first_case=pp.frames[vic_idx].head(1)['date']
last_case=pp.frames[vic_idx].tail(1)[['date', 'total', '7-day-projection-relative-error']]
pp.add_label(vic_idx, "2020-08-04", "peak daily new cases")
pp.add_label(vic_idx, "2020-07-16", "last model undershoot, prior to recovery")
pp.add_label(vic_idx, "2020-08-02", "stage 4 restrictions announced")
pp.add_label(vic_idx, "2020-07-15", "similar state (VIC)")
pp.add_label(vic_idx, "2020-07-09", "stage 3 restrictions announced")
pp.add_label(vic_idx, last_case.values[0][0], f"last day plotted (VIC)")
pp.add_label(vic_idx, first_case.values[0], f"first day plotted (VIC)")
# pp.add_label(vic_idx, "2020-08-09", "1 week after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-08-16", "2 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-08-23", "3 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-08-30", "4 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-09-06", "5 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-09-13", "6 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-09-20", "7 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-09-27", "8 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-10-03", "9 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-10-10", "10 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-10-17", "11 weeks after stage 4 restrictions announced")
# pp.add_label(vic_idx, "2020-10-24", "12 weeks after stage 4 restrictions announced")
first_case_nsw=pp.frames[bondi_idx].head(1)['date']
pp.add_label(bondi_idx, GRAPH_DATE, "current state (NSW)")
pp.add_label(bondi_idx, "2021-10-11", "relaxation (NSW)")
#pp.add_label(bondi_idx, PREV_WEEK, "a week ago (NSW)")
#pp.add_label(bondi_idx, PREV_FORTNIGHT, "two weeks ago (NSW)")
#pp.add_label(bondi_idx, first_case_nsw.values[0], f"first day plotted (NSW)")
pp.ax.grid()
pp.ax.figure.savefig(f'archive/{RUN_DATE}/hedgehog.png')
_=_
The decay rate estimates are calculated by takimg the Nth root of the ratio between two growth rate estimates taken on dates N days apart.
ax=plot_decay_rate_estimates(bondi, "Sydney 2021")
ax.figure.savefig(f"archive/{RUN_DATE}/decay-rate-sydney.png")
ax=plot_decay_rate_estimates(vic_outbreak, "Melbourne 2020")
ax.figure.savefig(f"archive/{RUN_DATE}/decay-rate-melbourne.png")
df=bondi
df['total-1'] = df['total'].rolling(window=5).mean()
shifted=df.shift(5)
df['total-2'] = shifted['Reff']*shifted["total"]
df['total-3'] = shifted['Reff']*shifted["total-1"]
df['total-4'] = df['total-2'].rolling(window=5).mean()
df['total-5'] = df['total-3'].rolling(window=5).mean()
cols=['total', 'total-1', 'total-2', 'total-3', 'total-4', 'total-5']
ax=df[cols].plot(figsize=(10,10))
_=ax.legend(['total', 'total (rolling average)', 'total (forecast#1)', 'total (forecast#2, rolling average)', 'total (rolling average forecast #1)', 'total (rolling average forecast #2)'])
cols=['total', 'total-2']
ax=df[cols].plot(figsize=(10,10))
_=ax.legend(['total (actual)', 'total (forecast)'])
cols=['total-1', 'total-4']
ax=df[cols].plot(figsize=(10,10))
ax.set_title("actual and forecast total (rolling averages)")
_=ax.legend(['total (actual)', 'total (forecast)'])
caution: These projections need to be taken with a large grain of salt.
We calculate a set of curves that estimate the growth decay rate.
We do this by taking the ratio of growth rate estimates for days separated by a period of n days for n=7..28, and then, assuming the change is due to exponential decay, take the nth root of that ratio to produce a daily growth decay factor estimate which is then converted into a rate estmate - the growth decay rate estimate.
We then take a weighted average of all such curves to produce a single curve which represents the weighted average estimate of the growth decay rate as a function of time. The weight for each curve is period n, so the curve is naturally weighted towards the longer period estimates. This seems reasonable, because the weighted average estimate falls within the bounds of the individual period based estimates.
Then, we calculate 4 statistics from that curve: min, max, mean and last.
Finally, we project the cumulative growth into the future using these 4 growth decay rate estimates to progressively update the growth rate.
It is expected that the last and mean statistics will converge over time and, if Melbourne 2020 is any guide, 'last' should start to drive beneath 'mean'. Until that happens, it seems safe to assume that the eventual result will be somewhere between the 'last' curve and the 'mean' curves since it seems improbable, at this point, that the growth decay rate will weaken further in the other directon.
ax=plot_new_cases_projection(bondi)
ax.figure.savefig(f"archive/{RUN_DATE}/new-cases-projections.png")
/Users/jon/shares/staff/personal/projects/bondi-outbreak/src/functions.py:640: RuntimeWarning: invalid value encountered in sqrt "mid": -np.sqrt(S.mean()*S.tail(1).values[0]) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs) /usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs)
The decay rate is calculated as an exponential fit through growth rates after July 2 (day 15).
gp=derive_growth_params(bondi[(bondi.date >= '2021-07-02')])
ax=bondi["ols-growth-rate"].plot(figsize=(10,10))
ax.plot((factor(gp[1])**bondi.index)*gp[0])
ax.set_title(f"Daily Cumulative Growth Rate % (Sydney 2021) ({RUN_DATE})")
ax.set_xlabel("Day Of Outbreak (t)")
ax.set_ylabel("Daily Cumulative Growth Rate % (gt)")
ax.legend(["observed", f"trend: gt=g0 * (1+r/100)^t"])
ax.text(34, 62, f"g0={round(gp[0],3)}, r={round(gp[1],3)}")
ax.grid()
ax.figure.savefig(f'archive/{RUN_DATE}/growth-rate-trend.png')
This plots the actual cumulative blue) and daily totals (orange) for the Sydney 2021 outbreak together with the implied daily growth rate (green) and the rate of decay of the growth rate (red</red>).
The dotted extensions of each observed plot are projections into the future assuming a constant (negative) decay rate. This decay rate is based on expoential fit through recent growth rate estimates (see plot above). This assumption is not sound and is likely to underestimate the true decay rate particularly in the later stages of the outbreak. As such peak estimates and timing are indicative only and are likely to be worst-case by a factor of 2 or more.
df=select_outbreak(project_ols_growth_rate_min(bondi, 365-len(bondi), gp[1]))
ax=plot_derivatives(df, len(bondi), "Sydney 2021")
ax.figure.savefig(f'archive/{RUN_DATE}/derivatives-sydney-partial.png')
factor(gp[1])
/usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs)
0.9882297396504495
df=vic_partial
gp=derive_growth_params(df[df.index>30])
ax=df["ols-growth-rate"].plot(figsize=(10,10))
ax.plot((factor(gp[1])**bondi.index)*gp[0])
ax.set_title(f"Daily Cumulative Growth Rate % (Melbourne 2021) ({RUN_DATE})")
ax.set_xlabel("Day Of Outbreak")
ax.set_ylabel("Daily Cumulative Growth Rate %")
ax.legend(["observed", f"trend: y={round(gp[0],3)} * ({round(gp[1],3)}^x)"])
#ax.figure.savefig(f'archive/{RUN_DATE}/growth-rate-trend.png')
<matplotlib.legend.Legend at 0x16b4bd1f0>
df=select_outbreak(project_ols_growth_rate_min(vic_partial, 110, gp[1]))
ax=plot_derivatives(df, len(vic_partial), dataset="Melbourne 2020")
#ax.figure.savefig(f'archive/{RUN_DATE}/derivatives-melbourne-partial.png')
/usr/local/lib/python3.9/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: invalid value encountered in log result = getattr(ufunc, method)(*inputs, **kwargs)
ax=plot_derivatives(vic_outbreak, None, dataset="Melbourne 2020")
ax.figure.savefig(f'archive/{RUN_DATE}/derivatives-melbourne-full.png')
The concurrent hospitalised and ICU statistics are the number of beds occupied by COVID-19 cases at the time of the cumulative case observation - it is not the total number of people admitted to hospital or ICU since the beginning of the outbreak. So if, over 7 days, there is a net increase of 1000 cumulative cases and 2 people leave the ICU and 5 people enter the ICU, then there will be a increase in occupancy of 3 people in the ICU - a rate of 3/1000 = 0.3%.
ax=plot_health_outcomes(bondi)
ax.figure.savefig(f'archive/{RUN_DATE}/health-outcomes.png')
ax.plot(bondi['pcr_cumulative'], bondi['hospitalised'], color="black", linestyle="dashed")
[<matplotlib.lines.Line2D at 0x16b1e3550>]
ax=bondi.shift(-186).plot(y="hospitalised", figsize=(10,10))
_=ax.plot(bondi[bondi.index<=97].shift(-46)['hospitalised'])
dec_x=15
june_x=50
level=1266
ax.vlines(dec_x, 200, level, color="C0", linestyle="dotted")
ax.vlines(june_x, 200, level, color="C1", linestyle="dotted")
ax.legend([
'december from 2021-12-20',
'june from 2021-08-02',
f'day={dec_x} (december)',
f'day={june_x} (june)',
])
ax.hlines(level, 0, 50, color='black', linestyle="dotted")
#ax.set_yscale('log')
_=ax.set_title("Comparison of hospitalisations - NSW Waves 2021")
ax.figure.savefig(f'archive/{RUN_DATE}/hospitalisation-velocity.png')
dec_x=10
june_x=19
june_peak_x=46
june_peak_y=244
crop_june=100
december = bondi.shift(-191)
december=december[december['icu'].notna()]
level=december.icu.max()
june=bondi[bondi.index<=crop_june].shift(-49)
june=june[june['icu'].notna()]
def get_parms(df, peak_y):
max=df['icu'].max()
min=df['icu'].min()
period=(df[df['icu']==max].index.min()-df[df['icu']==min].index.min())
gf=np.exp(np.log(max/min)/period)
peak_x=int(np.ceil(np.log(peak_y/min)/np.log(gf)))
min_date=df['date'].min()
return {
"max": max,
"min": min,
"period": period,
"gf": gf,
"peak_x": peak_x,
"min_date": min_date
}
p=get_parms(december, june_peak_y)
dec_min=p["min"]
dec_max=p["max"]
dec_period=p["period"]
dec_gf=p["gf"]
dec_peak_x = p["peak_x"]
dec_min_date = p["min_date"]
dec_proj=pd.DataFrame(index=range(0, dec_peak_x+1))
p=get_parms(june, june_peak_y)
june_min=p["min"]
june_max=p["max"]
june_period=p["period"]
june_gf=p["gf"]
june_peak_x = p["peak_x"]
june_min_date = p["min_date"]
june_proj=pd.DataFrame(index=range(0, june_peak_x+1))
df=june.copy().reindex()
df=df[df.index<=8]
p=get_parms(june, june_peak_y)
june_trunc_gf=p["gf"]
ax=december.rolling(1)['icu'].mean().plot(figsize=(10,10))
ax.plot(np.exp(np.log(dec_gf)*dec_proj.index)*dec_min, linestyle="dashed", color="C0")
ax.vlines(dec_x, 48, level, color='C0', linestyle="dotted")
ax.vlines(dec_peak_x, 48, 244, color="C0", linestyle="dotted")
ax.plot(june.rolling(1)['icu'].mean())
ax.plot(np.exp(np.log(june_trunc_gf)*june_proj.index)*june_min, linestyle="dashed", color="C1")
ax.vlines(june_x, 48, level, color='C1', linestyle="dotted")
ax.vlines(june_peak_x, 48, 244, color='C1', linestyle="dotted")
ax.legend([
'December Wave ICU from 2021-12-25',
"December (projection)",
'June Wave ICU from 2021-08-06',
"June (projection)",
f"{dec_x} days",
f"{dec_peak_x} days - {add_days(dec_min_date, dec_peak_x)}",
f"{june_x} days",
f"{june_peak_x} days"
])
ax.hlines(level, 0, 25, color='black', linestyle="dotted")
ax.hlines(244, 0, 55, color='black', linestyle="dotted")
#ax.set_yscale('log')
ax.set_title("Comparison of ICU - NSW Waves 2021")
ax.figure.savefig(f'archive/{RUN_DATE}/icu-velocity.png')
ax.text(36, 70, "June projection taken from day 8")
ax.set_yscale('log')
_=_
dec_gf
0.965579094854883
start_date='2021-08-20'
end_date=add_days(start_date, 9)
start_date='2021-12-25'
end_date=add_days(RUN_DATE, -1)
start=bondi[bondi['date']==start_date][['date', 'total', 'rolling-average', 'cumulative', 'hospitalised', 'icu']]
end=bondi[bondi['date']==end_date][['date', 'total', 'rolling-average', 'cumulative', 'hospitalised', 'icu']]
display(start,end)
end.index,\
end['total'].min()/start['total'].min(),\
end['rolling-average'].min()/start['rolling-average'].min(),\
end['cumulative'].min()/start['cumulative'].min(),\
end['hospitalised'].min()/start['hospitalised'].min(),\
end['icu'].min()/start['icu'].min()
date | total | rolling-average | cumulative | hospitalised | icu | |
---|---|---|---|---|---|---|
191 | 2021-12-25 | 6288 | 4214.571429 | 120004 | 388.0 | 52.0 |
date | total | rolling-average | cumulative | hospitalised | icu | |
---|---|---|---|---|---|---|
265 | 2022-03-09 | 14775 | 11273.0 | 1432855 | 1038.0 | 39.0 |
(Int64Index([265], dtype='int64'), 2.3497137404580153, 2.674767812351705, 11.94006033132229, 2.6752577319587627, 0.75)
bondi.tail(168)[["date", "cumulative", "rolling-average", "hospitalised-average", "total", "hospitalised", "deaths", "icu"]]
date | cumulative | rolling-average | hospitalised-average | total | hospitalised | deaths | icu | |
---|---|---|---|---|---|---|---|---|
98 | 2021-09-23 | 50857 | 1107.571429 | 1235.857143 | 1063 | 1244.0 | 266.0 | 233.0 |
99 | 2021-09-24 | 51900 | 1073.142857 | 1227.428571 | 1043 | 1186.0 | 277.0 | 232.0 |
100 | 2021-09-25 | 52907 | 1026.857143 | 1222.857143 | 1007 | 1187.0 | 288.0 | 229.0 |
101 | 2021-09-26 | 53868 | 1009.428571 | 1209.714286 | 961 | 1146.0 | 297.0 | 222.0 |
102 | 2021-09-27 | 54655 | 988.285714 | 1202.285714 | 787 | 1155.0 | 309.0 | 214.0 |
103 | 2021-09-28 | 55518 | 965.571429 | 1186.428571 | 863 | 1155.0 | 316.0 | 213.0 |
104 | 2021-09-29 | 56381 | 941.000000 | 1165.000000 | 863 | 1082.0 | 331.0 | 212.0 |
105 | 2021-09-30 | 57322 | 923.571429 | 1143.000000 | 941 | 1090.0 | 337.0 | 213.0 |
106 | 2021-10-01 | 58186 | 898.000000 | 1124.285714 | 864 | 1055.0 | 352.0 | 210.0 |
107 | 2021-10-02 | 58999 | 870.285714 | 1098.285714 | 813 | 1005.0 | 362.0 | 202.0 |
108 | 2021-10-03 | 59666 | 828.285714 | 1074.714286 | 667 | 981.0 | 372.0 | 195.0 |
109 | 2021-10-04 | 60289 | 804.857143 | 1046.714286 | 623 | 959.0 | 378.0 | 193.0 |
110 | 2021-10-05 | 60897 | 768.428571 | 1021.428571 | 608 | 978.0 | 385.0 | 190.0 |
111 | 2021-10-06 | 61491 | 730.000000 | 1003.857143 | 594 | 959.0 | 395.0 | 188.0 |
112 | 2021-10-07 | 62078 | 679.428571 | 978.285714 | 587 | 911.0 | 403.0 | 181.0 |
113 | 2021-10-08 | 62724 | 648.285714 | 949.857143 | 646 | 856.0 | 414.0 | 170.0 |
114 | 2021-10-09 | 63304 | 615.000000 | 922.285714 | 580 | 812.0 | 425.0 | 163.0 |
115 | 2021-10-10 | 63781 | 587.857143 | 895.571429 | 477 | 794.0 | 431.0 | 159.0 |
116 | 2021-10-11 | 64277 | 569.714286 | 868.428571 | 496 | 769.0 | 439.0 | 153.0 |
117 | 2021-10-12 | 64637 | 534.285714 | 838.142857 | 360 | 766.0 | 444.0 | 155.0 |
118 | 2021-10-13 | 65081 | 512.857143 | 803.428571 | 444 | 716.0 | 448.0 | 150.0 |
119 | 2021-10-14 | 65487 | 487.000000 | 774.857143 | 406 | 711.0 | 454.0 | 143.0 |
120 | 2021-10-15 | 65886 | 451.714286 | 749.285714 | 399 | 677.0 | 458.0 | 145.0 |
121 | 2021-10-16 | 66205 | 414.428571 | 726.428571 | 319 | 652.0 | 460.0 | 138.0 |
122 | 2021-10-17 | 66506 | 389.285714 | 701.428571 | 301 | 619.0 | 470.0 | 137.0 |
123 | 2021-10-18 | 66771 | 356.285714 | 678.142857 | 265 | 606.0 | 475.0 | 132.0 |
124 | 2021-10-19 | 67044 | 343.857143 | 652.857143 | 273 | 589.0 | 479.0 | 128.0 |
125 | 2021-10-20 | 67327 | 320.857143 | 629.428571 | 283 | 552.0 | 486.0 | 124.0 |
126 | 2021-10-21 | 67699 | 316.000000 | 602.571429 | 372 | 523.0 | 487.0 | 124.0 |
127 | 2021-10-22 | 68044 | 308.285714 | 574.714286 | 345 | 482.0 | 492.0 | 125.0 |
128 | 2021-10-23 | 68376 | 310.142857 | 548.571429 | 332 | 469.0 | 494.0 | 123.0 |
129 | 2021-10-24 | 68672 | 309.428571 | 528.714286 | 296 | 480.0 | 498.0 | 119.0 |
130 | 2021-10-25 | 68966 | 313.571429 | 509.857143 | 294 | 474.0 | 502.0 | 116.0 |
131 | 2021-10-26 | 69248 | 314.857143 | 491.000000 | 282 | 457.0 | 503.0 | 109.0 |
132 | 2021-10-27 | 69552 | 317.857143 | 471.857143 | 304 | 418.0 | 506.0 | 97.0 |
133 | 2021-10-28 | 69845 | 306.571429 | 451.571429 | 293 | 381.0 | 508.0 | 82.0 |
134 | 2021-10-29 | 70113 | 295.571429 | 434.571429 | 268 | 363.0 | 510.0 | 80.0 |
135 | 2021-10-30 | 70349 | 281.857143 | 416.571429 | 236 | 343.0 | 513.0 | 81.0 |
136 | 2021-10-31 | 70526 | 264.857143 | 396.571429 | 177 | 340.0 | 514.0 | 78.0 |
137 | 2021-11-01 | 70661 | 242.142857 | 378.714286 | 135 | 349.0 | 518.0 | 77.0 |
138 | 2021-11-02 | 70834 | 226.571429 | 361.000000 | 173 | 333.0 | 522.0 | 72.0 |
139 | 2021-11-03 | 71024 | 210.285714 | 345.428571 | 190 | 309.0 | 526.0 | 68.0 |
140 | 2021-11-04 | 71332 | 212.428571 | 334.142857 | 308 | 302.0 | 530.0 | 64.0 |
141 | 2021-11-05 | 71581 | 209.714286 | 323.000000 | 249 | 285.0 | 533.0 | 61.0 |
142 | 2021-11-06 | 71851 | 214.571429 | 312.571429 | 270 | 270.0 | 536.0 | 55.0 |
143 | 2021-11-07 | 72095 | 224.142857 | 302.428571 | 244 | 269.0 | 537.0 | 52.0 |
144 | 2021-11-08 | 72282 | 231.571429 | 290.285714 | 187 | 264.0 | 544.0 | 48.0 |
145 | 2021-11-09 | 72504 | 238.571429 | 279.000000 | 222 | 254.0 | 548.0 | 42.0 |
146 | 2021-11-10 | 72720 | 242.285714 | 268.428571 | 216 | 235.0 | 551.0 | 41.0 |
147 | 2021-11-11 | 72981 | 235.571429 | 257.857143 | 261 | 228.0 | 552.0 | 40.0 |
148 | 2021-11-12 | 73267 | 240.857143 | 250.857143 | 286 | 236.0 | 554.0 | 34.0 |
149 | 2021-11-13 | 73517 | 238.000000 | 244.857143 | 250 | 228.0 | 554.0 | 32.0 |
150 | 2021-11-14 | 73712 | 231.000000 | 238.428571 | 195 | 224.0 | 555.0 | 32.0 |
151 | 2021-11-15 | 73877 | 227.857143 | 231.571429 | 165 | 216.0 | 556.0 | 32.0 |
152 | 2021-11-16 | 74089 | 226.428571 | 225.285714 | 212 | 210.0 | 558.0 | 32.0 |
153 | 2021-11-17 | 74320 | 228.571429 | 221.285714 | 231 | 207.0 | 558.0 | 33.0 |
154 | 2021-11-18 | 74582 | 228.714286 | 217.142857 | 262 | 199.0 | 561.0 | 29.0 |
155 | 2021-11-19 | 74798 | 218.714286 | 211.428571 | 216 | 196.0 | 564.0 | 28.0 |
156 | 2021-11-20 | 74980 | 209.000000 | 206.142857 | 182 | 191.0 | 564.0 | 28.0 |
157 | 2021-11-21 | 75156 | 206.285714 | 201.571429 | 176 | 192.0 | 566.0 | 32.0 |
158 | 2021-11-22 | 75336 | 208.428571 | 199.571429 | 180 | 202.0 | 567.0 | 30.0 |
159 | 2021-11-23 | 75509 | 202.857143 | 197.571429 | 173 | 196.0 | 569.0 | 34.0 |
160 | 2021-11-24 | 75757 | 205.285714 | 195.857143 | 248 | 195.0 | 571.0 | 35.0 |
161 | 2021-11-25 | 76033 | 207.285714 | 194.714286 | 276 | 191.0 | 571.0 | 31.0 |
162 | 2021-11-26 | 76294 | 213.714286 | 193.142857 | 261 | 185.0 | NaN | 27.0 |
163 | 2021-11-27 | 76529 | 221.285714 | 190.714286 | 235 | 174.0 | 571.0 | 26.0 |
164 | 2021-11-28 | 76714 | 222.571429 | 186.857143 | 185 | 165.0 | 571.0 | 24.0 |
165 | 2021-11-29 | 76864 | 218.285714 | 182.285714 | 150 | 170.0 | 571.0 | 25.0 |
166 | 2021-11-30 | 77043 | 219.142857 | 177.142857 | 179 | 160.0 | NaN | 26.0 |
167 | 2021-12-01 | 77294 | 219.571429 | 171.285714 | 251 | 154.0 | 574.0 | 25.0 |
168 | 2021-12-02 | 77565 | 218.857143 | 164.571429 | 271 | 144.0 | 574.0 | 24.0 |
169 | 2021-12-03 | 77902 | 229.714286 | 158.142857 | 337 | 140.0 | 574.0 | 25.0 |
170 | 2021-12-04 | 78227 | 242.571429 | 153.142857 | 325 | 139.0 | 575.0 | 25.0 |
171 | 2021-12-05 | 78513 | 257.000000 | 150.714286 | 286 | 148.0 | 576.0 | 26.0 |
172 | 2021-12-06 | 78721 | 265.285714 | 148.142857 | 208 | 152.0 | 576.0 | 24.0 |
173 | 2021-12-07 | 78981 | 276.857143 | 147.428571 | 260 | 155.0 | 578.0 | 28.0 |
174 | 2021-12-08 | 79384 | 298.571429 | 147.000000 | 403 | 151.0 | 579.0 | 25.0 |
175 | 2021-12-09 | 79804 | 319.857143 | 148.000000 | 420 | 151.0 | 580.0 | 25.0 |
176 | 2021-12-10 | 80320 | 345.428571 | 150.571429 | 516 | 158.0 | 580.0 | 24.0 |
177 | 2021-12-11 | 80880 | 379.000000 | 152.142857 | 560 | 150.0 | 583.0 | 25.0 |
178 | 2021-12-12 | 81365 | 407.428571 | 153.285714 | 485 | 156.0 | 585.0 | 23.0 |
179 | 2021-12-13 | 81901 | 454.285714 | 156.000000 | 536 | 171.0 | 585.0 | 24.0 |
180 | 2021-12-14 | 82705 | 532.000000 | 157.857143 | 804 | 168.0 | 586.0 | 21.0 |
181 | 2021-12-15 | 84065 | 668.714286 | 160.000000 | 1360 | 166.0 | 587.0 | 24.0 |
182 | 2021-12-16 | 85807 | 857.571429 | 165.857143 | 1742 | 192.0 | 587.0 | 26.0 |
183 | 2021-12-17 | 88020 | 1100.000000 | 174.000000 | 2213 | 215.0 | 588.0 | 24.0 |
184 | 2021-12-18 | 90502 | 1374.571429 | 182.000000 | 2482 | 206.0 | 589.0 | 26.0 |
185 | 2021-12-19 | 93068 | 1671.857143 | 192.142857 | 2566 | 227.0 | 589.0 | 28.0 |
186 | 2021-12-20 | 95569 | 1952.571429 | 205.000000 | 2501 | 261.0 | 589.0 | 33.0 |
187 | 2021-12-21 | 98626 | 2274.428571 | 221.571429 | 3057 | 284.0 | 591.0 | 39.0 |
188 | 2021-12-22 | 102389 | 2617.714286 | 241.000000 | 3763 | 302.0 | 593.0 | 40.0 |
189 | 2021-12-23 | 108104 | 3185.285714 | 263.142857 | 5715 | 347.0 | 594.0 | 45.0 |
190 | 2021-12-24 | 113716 | 3670.857143 | 287.000000 | 5612 | 382.0 | 595.0 | 53.0 |
191 | 2021-12-25 | 120004 | 4214.571429 | 313.000000 | 6288 | 388.0 | 595.0 | 52.0 |
192 | 2021-12-26 | 126398 | 4761.428571 | 346.000000 | 6394 | 458.0 | 595.0 | 52.0 |
193 | 2021-12-27 | 132722 | 5307.571429 | 383.000000 | 6324 | 520.0 | 598.0 | 55.0 |
194 | 2021-12-28 | 138784 | 5736.857143 | 422.000000 | 6062 | 557.0 | 599.0 | 60.0 |
195 | 2021-12-29 | 149985 | 6799.428571 | 468.142857 | 11201 | 625.0 | 602.0 | 61.0 |
196 | 2021-12-30 | 162211 | 7729.571429 | 525.142857 | 12226 | 746.0 | 603.0 | 63.0 |
197 | 2021-12-31 | 183362 | 9949.428571 | 589.428571 | 21151 | 832.0 | 609.0 | 69.0 |
198 | 2022-01-01 | 210627 | 12946.142857 | 662.714286 | 27265 | 901.0 | 614.0 | 79.0 |
199 | 2022-01-02 | 232700 | 15186.000000 | 749.571429 | 22073 | 1066.0 | 615.0 | 83.0 |
200 | 2022-01-03 | 257812 | 17870.000000 | 847.285714 | 25112 | 1204.0 | 619.0 | 95.0 |
201 | 2022-01-04 | 285746 | 20994.571429 | 959.714286 | 27934 | 1344.0 | 621.0 | 105.0 |
202 | 2022-01-05 | 328079 | 25442.000000 | 1083.428571 | 42333 | 1491.0 | 629.0 | 119.0 |
203 | 2022-01-06 | 370339 | 29732.571429 | 1206.714286 | 42260 | 1609.0 | 635.0 | 131.0 |
204 | 2022-01-07 | 423024 | 34237.428571 | 1336.142857 | 52685 | 1738.0 | 646.0 | 134.0 |
205 | 2022-01-08 | 486549 | 39417.428571 | 1463.857143 | 63525 | 1795.0 | 655.0 | 145.0 |
206 | 2022-01-09 | 530633 | 42561.857143 | 1586.857143 | 44084 | 1927.0 | 671.0 | 151.0 |
207 | 2022-01-10 | 561436 | 43374.857143 | 1704.857143 | 30803 | 2030.0 | 689.0 | 159.0 |
208 | 2022-01-11 | 602004 | 45179.714286 | 1825.142857 | 40568 | 2186.0 | 700.0 | 170.0 |
209 | 2022-01-12 | 658604 | 47217.857143 | 1932.428571 | 56600 | 2242.0 | 721.0 | 175.0 |
210 | 2022-01-13 | 710182 | 48549.000000 | 2043.000000 | 51578 | 2383.0 | 743.0 | 182.0 |
211 | 2022-01-14 | 748666 | 46520.285714 | 2155.428571 | 38484 | 2525.0 | 772.0 | 184.0 |
212 | 2022-01-15 | 788703 | 43164.857143 | 2267.000000 | 40037 | 2576.0 | 792.0 | 193.0 |
213 | 2022-01-16 | 819376 | 41249.000000 | 2370.285714 | 30673 | 2650.0 | 812.0 | 191.0 |
214 | 2022-01-17 | 845060 | 40517.714286 | 2476.857143 | 25684 | 2776.0 | 829.0 | 203.0 |
215 | 2022-01-18 | 868509 | 38072.142857 | 2571.714286 | 23449 | 2850.0 | 863.0 | 209.0 |
216 | 2022-01-19 | 897838 | 34176.285714 | 2660.428571 | 29329 | 2863.0 | 897.0 | 217.0 |
217 | 2022-01-20 | 924200 | 30574.000000 | 2717.285714 | 26362 | 2781.0 | 922.0 | 212.0 |
218 | 2022-01-21 | 947062 | 28342.285714 | 2748.428571 | 22862 | 2743.0 | 968.0 | 209.0 |
219 | 2022-01-22 | 964919 | 25173.714286 | 2775.000000 | 17857 | 2762.0 | 998.0 | 204.0 |
220 | 2022-01-23 | 986121 | 23820.714286 | 2783.857143 | 21202 | 2712.0 | 1032.0 | 189.0 |
221 | 2022-01-24 | 999212 | 22021.714286 | 2789.571429 | 13091 | 2816.0 | 1056.0 | 196.0 |
222 | 2022-01-25 | 1014558 | 20864.142857 | 2802.857143 | 15346 | 2943.0 | 1085.0 | 183.0 |
223 | 2022-01-26 | 1036283 | 19777.857143 | 2793.000000 | 21725 | 2794.0 | 1113.0 | 175.0 |
224 | 2022-01-27 | 1053428 | 18461.142857 | 2784.571429 | 17145 | 2722.0 | 1142.0 | 181.0 |
225 | 2022-01-28 | 1064649 | 16798.142857 | 2783.714286 | 11221 | 2737.0 | 1212.0 | 189.0 |
226 | 2022-01-29 | 1079211 | 16327.428571 | 2773.857143 | 14562 | 2693.0 | 1260.0 | 186.0 |
227 | 2022-01-30 | 1093411 | 15327.142857 | 2766.857143 | 14200 | 2663.0 | 1312.0 | 182.0 |
228 | 2022-01-31 | 1107576 | 15480.571429 | 2761.571429 | 14165 | 2779.0 | 1338.0 | 185.0 |
229 | 2022-02-01 | 1117188 | 14661.428571 | 2733.857143 | 9612 | 2749.0 | 1368.0 | 186.0 |
230 | 2022-02-02 | 1128495 | 13173.142857 | 2709.285714 | 11307 | 2622.0 | 1395.0 | 170.0 |
231 | 2022-02-03 | 1140275 | 12406.714286 | 2688.714286 | 11780 | 2578.0 | 1433.0 | 160.0 |
232 | 2022-02-04 | 1150389 | 12248.571429 | 2654.000000 | 10114 | 2494.0 | 1464.0 | 160.0 |
233 | 2022-02-05 | 1157423 | 11173.142857 | 2603.142857 | 7034 | 2337.0 | 1482.0 | 152.0 |
234 | 2022-02-06 | 1165711 | 10328.571429 | 2554.285714 | 8288 | 2321.0 | 1510.0 | 147.0 |
235 | 2022-02-07 | 1174093 | 9502.428571 | 2457.142857 | 8382 | 2099.0 | 1524.0 | 137.0 |
236 | 2022-02-08 | 1180830 | 9091.714286 | 2359.857143 | 6737 | 2068.0 | 1542.0 | 132.0 |
237 | 2022-02-09 | 1191221 | 8960.857143 | 2257.571429 | 10391 | 1906.0 | 1562.0 | 132.0 |
238 | 2022-02-10 | 1201756 | 8783.000000 | 2145.714286 | 10535 | 1795.0 | 1586.0 | 121.0 |
239 | 2022-02-11 | 1210810 | 8631.571429 | 2034.571429 | 9054 | 1716.0 | 1605.0 | 108.0 |
240 | 2022-02-12 | 1218768 | 8763.571429 | 1936.428571 | 7958 | 1650.0 | 1637.0 | 104.0 |
241 | 2022-02-13 | 1224549 | 8405.428571 | 1850.000000 | 5781 | 1716.0 | 1659.0 | 93.0 |
242 | 2022-02-14 | 1233897 | 8543.428571 | 1785.714286 | 9348 | 1649.0 | 1673.0 | 100.0 |
243 | 2022-02-15 | 1237222 | 8056.000000 | 1716.428571 | 3325 | 1583.0 | 1689.0 | 96.0 |
244 | 2022-02-16 | 1245685 | 7780.571429 | 1655.285714 | 8463 | 1478.0 | 1716.0 | 92.0 |
245 | 2022-02-17 | 1253790 | 7433.428571 | 1605.571429 | 8105 | 1447.0 | 1730.0 | 92.0 |
246 | 2022-02-18 | 1263209 | 7485.571429 | 1557.714286 | 9419 | 1381.0 | 1745.0 | 92.0 |
247 | 2022-02-19 | 1268004 | 7033.714286 | 1507.285714 | 4795 | 1297.0 | 1757.0 | 81.0 |
248 | 2022-02-20 | 1274061 | 7073.142857 | 1445.000000 | 6057 | 1280.0 | 1778.0 | 77.0 |
249 | 2022-02-21 | 1279788 | 6555.857143 | 1393.428571 | 5727 | 1288.0 | 1785.0 | 74.0 |
250 | 2022-02-22 | 1289261 | 7434.142857 | 1352.000000 | 9473 | 1293.0 | 1799.0 | 71.0 |
251 | 2022-02-23 | 1298010 | 7475.000000 | 1318.857143 | 8749 | 1246.0 | 1805.0 | 69.0 |
252 | 2022-02-24 | 1306125 | 7476.428571 | 1285.142857 | 8115 | 1211.0 | 1817.0 | 59.0 |
253 | 2022-02-25 | 1313521 | 7187.428571 | 1251.285714 | 7396 | 1144.0 | 1823.0 | 64.0 |
254 | 2022-02-26 | 1320327 | 7474.714286 | 1227.428571 | 6806 | 1130.0 | 1834.0 | 59.0 |
255 | 2022-02-27 | 1326248 | 7455.285714 | 1208.285714 | 5921 | 1146.0 | 1841.0 | 58.0 |
256 | 2022-02-28 | 1332591 | 7543.285714 | 1186.571429 | 6343 | 1136.0 | 1847.0 | 55.0 |
257 | 2022-03-01 | 1342973 | 7673.142857 | 1158.714286 | 10382 | 1098.0 | 1856.0 | 49.0 |
258 | 2022-03-02 | 1353944 | 7990.571429 | 1133.857143 | 10971 | 1072.0 | 1861.0 | 45.0 |
259 | 2022-03-03 | 1365808 | 8526.142857 | 1108.714286 | 11864 | 1035.0 | 1870.0 | 43.0 |
260 | 2022-03-04 | 1375172 | 8807.285714 | 1088.142857 | 9364 | 1000.0 | 1872.0 | 42.0 |
261 | 2022-03-05 | 1385422 | 9299.285714 | 1068.857143 | 10250 | 995.0 | 1882.0 | 45.0 |
262 | 2022-03-06 | 1393949 | 9671.571429 | 1049.285714 | 8527 | 1009.0 | 1887.0 | 43.0 |
263 | 2022-03-07 | 1403189 | 10085.428571 | 1039.285714 | 9240 | 1066.0 | 1892.0 | 43.0 |
264 | 2022-03-08 | 1418080 | 10729.571429 | 1035.285714 | 14891 | 1070.0 | 1897.0 | 43.0 |
265 | 2022-03-09 | 1432855 | 11273.000000 | 1030.428571 | 14775 | 1038.0 | 1906.0 | 39.0 |
ax=bondi.plot(x='rolling-average', y='hospitalised', figsize=(10,10))
ax.plot
ax.set_title("NSW Currently Hospitalised vs Case Total (rolling 7 day average)")
ax.set_xlabel("Daily Cases (7-day rolling average)")
ax.set_ylabel("Currently Hospitalised")
#ax.scatter(x=bondi['rolling-average'], y=bondi['hospitalised'])
legends=["hospitalised"]
for (start,end) in [
('2021-06-16', '2021-09-11'),
('2021-09-11', '2021-09-25'),
('2021-09-18', '2021-10-18'),
('2021-10-11', '2021-12-05'),
('2021-12-05', RUN_DATE),
]:
out,m,b=hospitalised_vs_cases(bondi, start,end)
ax.plot(out['rolling-average'], out['hospitalised'])
legends.append(f"{start} - {end} m={round(m,3)}, b={round(b,2)}")
out=pd.DataFrame(index=range(0,100,1))
min=bondi['rolling-average'].min()
max=bondi['rolling-average'].max()
#max=25000
out['rolling-average']=out.index*(max-min)/100+min
out['hospitalised']=out['rolling-average']*m+b
ax.plot(out['rolling-average'], out['hospitalised'], linestyle="dashed")
_=ax.legend(legends)
ax.plot(bondi['pcr'].rolling(window=7).mean(), bondi['hospitalised'], linestyle="dashed")
ax.figure.savefig(f'archive/{RUN_DATE}/hospitalised-phase-plot.png')
ax=bondi.plot(x='rolling-average', y='hospitalised-average', figsize=(10,10))
ax.plot
ax.set_title("NSW Currently Hospitalised vs Case Total (both rolling 7 day average)")
ax.set_xlabel("Daily Cases (7-day rolling average)")
ax.set_ylabel("Currently Hospitalised (7-day rolling average)")
#ax.scatter(x=bondi['rolling-average'], y=bondi['hospitalised'])
legends=["hospitalised-average"]
for (start,end) in [
('2021-06-16', '2021-09-11'),
('2021-09-11', '2021-09-25'),
('2021-09-18', '2021-10-18'),
('2021-10-11', '2021-12-05'),
('2021-12-05', RUN_DATE),
]:
out,m,b=hospitalised_vs_cases(bondi, start,end,hospitalised='hospitalised-average')
ax.plot(out['rolling-average'], out['hospitalised-average'])
legends.append(f"{start} - {end} m={round(m,3)}, b={round(b,2)}")
out=pd.DataFrame(index=range(0,100,1))
min=bondi['rolling-average'].min()
max=bondi['rolling-average'].max()
#max=25000
out['rolling-average']=out.index*(max-min)/100+min
out['hospitalised-average']=out['rolling-average']*m+b
ax.plot(out['rolling-average'], out['hospitalised-average'], linestyle="dashed")
_=ax.legend(legends)
ax.figure.savefig(f'archive/{RUN_DATE}/hospitalised-average-phase-plot.png')
ax=bondi.plot(x='rolling-average', y='icu', figsize=(10,10))
ax.plot
ax.set_title("NSW Currently In ICU vs Case Total (rolling 7 day average)")
ax.set_xlabel("Daily Cases (7-day rolling average)")
ax.set_ylabel("Currently In ICU")
#ax.scatter(x=bondi['rolling-average'], y=bondi['hospitalised'])
legends=["icu"]
for (start,end) in [
('2021-06-16', '2021-09-11'),
('2021-09-11', '2021-09-25'),
('2021-09-18', '2021-10-18'),
('2021-10-11', '2021-12-05'),
('2021-12-05', RUN_DATE),
]:
out,m,b=icu_vs_cases(bondi, start,end)
ax.plot(out['rolling-average'], out['icu'])
legends.append(f"{start} - {end} m={round(m,3)}, b={round(b,2)}")
out=pd.DataFrame(index=range(0,100,1))
min=bondi['rolling-average'].min()
max=bondi['rolling-average'].max()
#max=25000
out['rolling-average']=out.index*(max-min)/100+min
out['icu']=out['rolling-average']*m+b
ax.plot(out['rolling-average'], out['icu'], linestyle="dashed")
_=ax.legend(legends)
ax.figure.savefig(f'archive/{RUN_DATE}/icu-phase-plot.png')
#ax.set_yscale('log')
#ax.set_xscale('log')
#slice[['log-cumulative', 'deaths', 'log-deaths', 'log-deaths-fit']]
bondi["icu-average"] = bondi["icu"].rolling(7).mean()
ax=bondi.plot(x='rolling-average', y='icu-average', figsize=(10,10))
ax.plot
ax.set_title("NSW Currently In ICU vs Case Total (both rolling 7 day average)")
ax.set_xlabel("Daily Cases (7-day rolling average)")
ax.set_ylabel("Currently In ICU (7-day rolling average)")
#ax.scatter(x=bondi['rolling-average'], y=bondi['hospitalised'])
legends=["icu-average"]
for (start,end) in [
('2021-06-16', '2021-09-11'),
('2021-09-11', '2021-09-25'),
('2021-09-18', '2021-10-18'),
('2021-10-11', '2021-12-05'),
('2021-12-05', RUN_DATE),
]:
out,m,b=icu_vs_cases(bondi, start,end, icu='icu-average')
ax.plot(out['rolling-average'], out['icu-average'])
legends.append(f"{start} - {end} m={round(m,3)}, b={round(b,2)}")
out=pd.DataFrame(index=range(0,100,1))
min=bondi['rolling-average'].min()
max=bondi['rolling-average'].max()
#max=25000
out['rolling-average']=out.index*(max-min)/100+min
out['icu-average']=out['rolling-average']*m+b
ax.plot(out['rolling-average'], out['icu-average'], linestyle="dashed")
_=ax.legend(legends)
ax.figure.savefig(f'archive/{RUN_DATE}/icu-average-phase-plot.png')
#slice[['log-cumulative', 'deaths', 'log-deaths', 'log-deaths-fit']]
legends=['icu']
ax=bondi.plot(x='hospitalised', y='icu', figsize=(10,10))
ax.set_title("ICU beds occupied vs Hospital Beds occupied (NSW, since 2021-06-16)")
for (start,end) in [
('2021-06-16', '2021-09-11'),
('2021-09-11', '2021-09-25'),
('2021-09-18', '2021-10-18'),
('2021-10-11', '2021-12-05'),
('2021-12-05', RUN_DATE),
]:
out,m,b=icu_vs_hospitalised(bondi, start,end)
ax.plot(out['hospitalised'], out['icu'])
legends.append(f"{start} - {end} m={round(m,3)}, b={round(b,2)}")
ax.legend(legends)
_=ax.set_ylabel('icu')
ax.figure.savefig(f'archive/{RUN_DATE}/icu-hospitalisation-phase-plot.png')
%%capture phaseplot
animate_phaseplot(
df=bondi,
outbreak="Sydney 2021",
fn=f'archive/{RUN_DATE}/animated-hedgehog-sydney2021.gif',
offset=15
)
if ENABLE_MELBOURNE_ANIMATION:
animate_phaseplot(
df=vic_outbreak,
offset=20,
outbreak="Melbourne 2020",
fn=f'archive/{RUN_DATE}/animated-hedgehog-melbourne2020.gif'
)
display(HTML(f"<img src='animated-hedgehog-sydney2021.gif'>"))
display(HTML(f"<img src='../latest/animated-hedgehog-melbourne2020.gif'>"))
%%capture derivatives
animate_derivatives(
bondi,
"Sydney 2021",
f'archive/{RUN_DATE}/animated-derivatives-sydney.gif',
derive_growth_params(bondi[(bondi.index>15)])[1]
)
if ENABLE_MELBOURNE_ANIMATION:
animate_derivatives(
vic_outbreak,
"Melbourne 2020",
f'archive/{RUN_DATE}/animated-derivatives-melbourne.gif',
derive_growth_params(vic_outbreak[(vic_outbreak.index>30)])[1]
)
display(HTML(f"<img src='animated-derivatives-sydney.gif'>"))
display(HTML(f"<img src='../latest/animated-derivatives-melbourne.gif'>"))
The minutephysics YouTube channel popularised the technique of using log-log plots to visualise when exponential growth of daily case totals slows and then reverses.
This section creates static and animated log-log plots. In addition to plotting the log-log curve, we fit trend lines to the last 14 days to better help visualise changes in growth rates. This is most evident in the animated plots.
ax=plot_log_log_day(bondi.index.max(), vic_outbreak, bondi)
ax.figure.savefig(f"archive/{RUN_DATE}/log-log.png")
Also, we plot the calculated gradients against time to see how they change with time.
Note: we calculate and plot the gradient of the natural log-log curves not the log base 10 curves.
df=bondi
ax=df["ltlc-gradient"].plot(figsize=(10,10))
max_x=df.tail(1).index.values[0]+15
ax.plot(vic_outbreak["ltlc-gradient"])
ax.hlines(0,xmin=0,xmax=max_x, color="black", linestyle="dotted")
ax.set_ylim(top=20, bottom=-20)
ax.set_xlim(left=0, right=max_x)
ax.set_ylabel("ln-ln gradient")
ax.set_xlabel("day of outbreak")
ax.set_title("gradient of ln daily cases vs ln cumulative cases plot vs day of outbreak")
ax.grid()
# df=vic_outbreak[vic_outbreak.index < 60]
# df=bondi
# endog=df['ltlc-gradient']
# exog=sm.add_constant(df.index)
# window=31
# model = RollingOLS(endog=endog, exog=exog, window=window, min_nobs=window)
# params=model.fit().params.tail(1)
# ax.plot((vic_outbreak.index*params["x1"].values[0]+params["const"].values[0]))
import statsmodels.api as sm
slice=bondi[bondi.index>65]
model=sm.OLS(endog=slice['ltlc-gradient'], exog=sm.add_constant(slice.index))
b,m=model.fit().params.const,model.fit().params.x1
projection=np.array([r for r in range(65, max_x)])
x0=-b/m
ax.plot(projection, m*projection+b, linestyle="dashed", color="C0")
ax.set_xticks([r for r in range(0, max_x, 10)])
ax.vlines(df[df['date']=="2021-10-11"].index, -20, 20, linestyle="dashed", color="C2")
ax.figure.savefig(f"archive/{RUN_DATE}/ltlc-gradients.png")
#print(x0)
%%capture log_log
animate_log_log_plot(vic_outbreak, bondi, f"archive/{RUN_DATE}/animated-log-log.gif")
display(HTML(f"<img src='animated-log-log.gif'>"))
%%capture new_cases
_=animate_new_cases_plot(bondi, 21, f"archive/{RUN_DATE}/animated-new-cases.gif")
display(HTML(f"<img src='animated-new-cases.gif'>"))
ax=plot_linear_growth(bondi)
ax=plot_linear_growth_error(bondi)
ax=plot_linear_growth(vic_outbreak)
ax=plot_linear_growth_error(vic_outbreak)
slice=bondi[bondi.index>=171]
december=slice
june=bondi[(bondi.index>=25) & (bondi.index<75)]
june_icu=bondi[(bondi.index>=60) & (bondi.index<90)]
dec_hosp=slice[slice.index>180]
dec_icu=slice[slice.index>185]
early=bondi[(bondi['date'] >= '2021-06-16') & (bondi['date'] <= '2021-09-11')]
late=bondi[(bondi['date'] >= '2021-09-18') & (bondi['date'] <= '2021-10-18')]
#early
def fit(df, col):
model=sm.OLS(endog=np.log(df[col]), exog=sm.add_constant(df.index))
b,m=model.fit().params.const,model.fit().params.x1
avg,v=slice.head()[['rolling-average', col]].min()
intercept=v-(avg*m)
return m,b
ax=bondi[['rolling-average', 'hospitalised', 'icu']].plot(figsize=(10,10))
ax.set_yscale('log')
m_cases, b_cases=fit(slice, 'rolling-average')
ax.plot(slice.index, np.exp(slice.index*m_cases+b_cases), color="C0")
m_cases_june, b_cases_june=fit(june, 'rolling-average')
ax.plot(june.index, np.exp(june.index*m_cases_june+b_cases_june), color="C0")
m_hosp, b_hosp=fit(dec_hosp, 'hospitalised')
ax.plot(dec_hosp.index, np.exp(dec_hosp.index*m_hosp+b_hosp), color="C1")
m_hosp_june, b_hosp_june=fit(june, 'hospitalised')
ax.plot(june.index, np.exp(june.index*m_hosp_june+b_hosp_june), color="C1")
m_icu, b_icu=fit(dec_icu, 'icu')
ax.plot(dec_icu.index, np.exp(dec_icu.index*m_icu+b_icu), color="C2")
m_icu_june, b_icu_june=fit(june_icu, 'icu')
ax.plot(june_icu.index, np.exp(june_icu.index*m_icu_june+b_icu_june), color="C2")
x=december.min()['rolling-average']
y=december.min().icu
m=0.1275
b=y-m*x
x,y,m,b,m*x+b,(y-b)/m
ax.plot(december.index, (december.icu-b)/m-125, linewidth=1, linestyle="dashed")
ax.plot(december.index, december['rolling-average']-(december.icu-b)/m)
ax.plot(early.index, (early.icu-b)/m-125, color="C3", linewidth=1, linestyle="dashed")
ax.plot(late.index, (late.icu-b)/m-500, color="C3", linewidth=1, linestyle="dashed")
ax.legend([
f"dec: daily cases (7-day average) g={round(np.exp(m_cases),3)}, A={round(np.exp(b_cases),12)}",
f"dec: hospitalised g={round(np.exp(m_hosp),3)}, A={round(np.exp(b_hosp),7)}",
f"dec: icu g={round(np.exp(m_icu),3)}, A={round(np.exp(b_icu),5)}",
f"june: daily cases (7-day average) g={round(np.exp(m_cases_june),3)}, A={round(np.exp(b_cases_june),12)}",
f"june: hospitalised g={round(np.exp(m_hosp_june),3)}, A={round(np.exp(b_hosp_june),5)}",
f"june: icu g={round(np.exp(m_icu_june),3)}, A={round(np.exp(b_icu_june),5)}",
f"hypothesis: Delta daily cases (7-day average) corrected assuming delta explains all ICU",
f"hypothesis: Omicron daily cases (7-day average) corrected assuming delta explains all ICU",
"???"
])
<matplotlib.legend.Legend at 0x17be46d00>
add_days('2021-06-16', 117)
'2021-10-11'
slice=bondi[bondi.date >= '2021-12-15']
col='hospitalised'
model=sm.OLS(endog=np.log(slice[col]), exog=sm.add_constant(slice.index))
b,m=model.fit().params.const,model.fit().params.x1
b,m
(4.119902783148169, 0.01359990219258975)
ax=slice[col].plot()
ax.set_yscale("log")
ax.set_ylabel(col)
ax.set_xlabel("day")
ax.set_title(f"{col} vs day")
ax.text(195, 200, f"day {slice.index.values[0]} = {slice.date.head(1).values[0]}")
ax.plot(slice.index, np.exp(slice.index*m+b))
ax.legend(["hospitalised (actual)", f"hospitalised = e^({round(m,3)}*t+{round(b,3)})"])
<matplotlib.legend.Legend at 0x17bf8d430>
df=bondi
endog=df["hospitalised"]#.rolling(window=7).mean()
window=7
x='rolling-average'
#x='total'
exog=sm.add_constant(df[x])
model = RollingOLS(endog=endog, exog=exog, window=window, min_nobs=window)
params=model.fit().params
params['date']=df['date']
#params[params.date >= '2021-10-11'][x].plot()
#np.log(params[x]).plot()
params[x].plot()
<AxesSubplot:>
params[params.date > '2021-12-12'][x].plot()
<AxesSubplot:>
bondi[['pcr', 'rat', 'rat_added', 'rat_added_last_week']][bondi['rat']>0].plot(figsize=(10,10))
<AxesSubplot:>