Sulla base del tasso di crescita λ e su una stima del tempo di generazione g si calcolano l'indice Rt e le fasce di confidenza al 68% degli ultimi 30 giorni. Gli esiti vengono riportati a video nella forma
------------------------------
data Rt 68% [ min, max ]
------------------------------
Carica i moduli necessari
import json
import numpy as np
from scipy.stats import linregress
import datetime
La funzione partiziona()
partiziona xDati
in un array di array. Ciascun elemento ha lunghezza intervalloGiorni
e il loro numero è pari a n-intervalloGiorni+1
. Segue la funzione per il calcolo di Rt.
def partiziona(xdati, intervalloGiorni):
numElem = len(xdati)
elementi = []
for i in range(intervalloGiorni, numElem + 1):
elementi = np.append(elementi, xdati[i-intervalloGiorni:i])
return np.reshape(elementi, (numElem - intervalloGiorni + 1, intervalloGiorni))
def rT(g, l):
return np.exp(g*l)
######### Predisposizioni ###############################################################################################
Elementi da inserire in input per avviare l'elaborazione dei dati nazionali o regionali. Questi devono essere stati preventivamente memorizzati nella medesima cartella e devono corrispondere alla data immessa (avviare eventualmente il notebook prelevaRinomina.ipynb
).
print()
print("Si intendono elaborare i nuovi positivi nazionali (i) o quelli di una regione (r)? ")
scelta = input(" :> ")
print("Inserire la data nel formato (YYYYMMDD) ")
dataISO = input(" :> ")
Vengono composti i nomi dei file e, una volta letti i dati contenuti, questi vengono caricati nell'array nuoviPositivi
.
if scelta == "i":
nomeFile = 'datiNazionali' + dataISO
# lettura del file csv: la colonna 8 corrisponde al campo nuovi_positivi
nuoviPositivi = np.genfromtxt(nomeFile + '.csv', delimiter=',', skip_header=1, usecols=8, dtype=int)
else:
print("Inserire la regione ")
regioneScelta = input(" :> ")
regioneScelta = regioneScelta.lower().capitalize()
nomeFile = 'datiRegionali' + dataISO
# lettura del file json e riportato l'array ad un array di Numpy
with open(nomeFile + '.json') as f:
datiGrezzi = json.load(f)
nuoviPositivi = []
for record in datiGrezzi:
if record['denominazione_regione'] == regioneScelta:
nuoviPositivi.append(record['nuovi_positivi'])
nuoviPositivi = np.array(nuoviPositivi)
Dovendo calcolare il logaritmo dei nuoviPositivi
nel caso vi sia un valore giornaliero nullo (per cui il suo logaritmo non avrebbe significato), lo si pone pari ad 1. L'array giorniDelFit
rappresenta i giorni nei quali si eseguirà il fit lineare.
nuoviPositivi[nuoviPositivi == 0] = 1
logPositivi = np.log(nuoviPositivi)
numGiorni = len(nuoviPositivi)
intervalloGiorni = 14
# i giorni da 0 a 13 nel quale viene eseguito il fit
giorniDelFit = np.arange(0, intervalloGiorni, 1)
Per ognuno degli intervalli in cui è stato partizionato l'insieme dei logPositivi
, viene eseguita una regressione lineare in un intervallo temporale di 14 giorni. Della retta ottimale si determina la sola pendenza in quanto questo è il termine associato, nel modello teorico, al parametro λ (o tasso di crescita).
valoriLambdaLin = []
for i in range(intervalloGiorni, numGiorni+1):
# fit lineare
esito = linregress(giorniDelFit, logPositivi[i-intervalloGiorni:i])
valoriLambdaLin = np.append(valoriLambdaLin, esito.slope)
generationTime = 5.8
valoriMediDayRtLin = rT(generationTime, valoriLambdaLin)
media14Giorni = np.mean(partiziona(valoriMediDayRtLin, intervalloGiorni), axis=1)
generationTimeError = 1.88
valoriMediDayRtLinSup = rT(generationTime+generationTimeError, valoriLambdaLin)
valoriMediDayRtLinInf = rT(generationTime-generationTimeError, valoriLambdaLin)
media14GiorniSup = np.mean(partiziona(valoriMediDayRtLinSup, intervalloGiorni), axis=1)
media14GiorniInf = np.mean(partiziona(valoriMediDayRtLinInf, intervalloGiorni), axis=1)
La variabile date1
definisce la data cui associare i valori giornalieri di Rt. Qui, convenzionalmente, associamo il valore dell'indice alla data associata ai file iniziali in modo da reneder più immediato il confronto con i valori presenti nel foglio nuoviPositivi-IT.xlsx
. Per un confronto con CovidStat è sufficiente sottrarre 4 giorni.
date1 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)
for i in range(numGiorni-intervalloGiorni-6)]
Si riportano a video i risultati degli ultimi 30 giorni con un arrotondamento a due cifre decimali sia per Rt che per i corrispondenti estremi degli intervalli di confidenza al 68%.
print()
print(" data Rt int.confidenza")
print("-----------------------------------------")
for i in range(30):
if media14Giorni[-30:][i] >= 1:
print(date1[-30:][i].strftime("%d/%m/%y"), ' ', '{:.2f}'.format(media14Giorni[-30:][i]), " 68% [",
'{:.2f}'.format(media14GiorniInf[-30:][i]), ",", '{:.2f}'.format(media14GiorniSup[-30:][i]), "]")
else:
print(date1[-30:][i].strftime("%d/%m/%y"), ' ', '{:.2f}'.format(media14Giorni[-30:][i]), " 68% [",
'{:.2f}'.format(media14GiorniSup[-30:][i]), ",", '{:.2f}'.format(media14GiorniInf[-30:][i]), "]")
print()
[1] G. Bonifazi, L. Lista et al., A simplified estimate of the Effective Reproduction Rt Number using its relation with the doubling time and application to Italian COVID-19 data, The European Physical Journal Plus (2021)
[2] D. Cereda et al., The early phase of the COVID-19 outbreak in Lombardy, Italy, arXiv:2003.09320 (2020)