Sulla base del tasso di crescita λ (notebook tassoCrescitaLambda.ipynb
) e su una stima del tempo di generazione si calcola l'indice Rt. Data l'incertezza del tempo di generazione, si rappresentano le fasce di confidenza al 68%.
Carica i moduli necessari
import json
import matplotlib.pyplot as plt
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
).
Si richiedono inoltre l'unità temporale da porre in ascissa e l'eventuale visualizzazione delle fasce di confidenza.
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(" :> ")
print("Si vuole l'asse temporale espresso in mesi (m) o in giorni trascorsi (g) dall'inizio pandemia? ")
unitaAssex = input(" :> ")
print("Si vogliono visualizzare le fasce di errore (s/n)? ")
fasce = input(" :> ")
Vengono composti i nomi dei file e, una volta letti i dati contenuti, questi vengono caricati nell'array nuoviPositivi
.
if scelta == "i":
regioneScelta = "ITALIA"
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)
if fasce == 's':
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)
Predisposizioni per l'asse x associato ai giorni trascorsi dal 24/02/2020 e valore massimo per l'asse y associato all'indice Rt. Codifica italiana della data inserita inizialmente.
xValoriMediGiornalieri = np.arange(intervalloGiorni-1, numGiorni)
xMediaMobile = np.arange(intervalloGiorni+5, numGiorni-intervalloGiorni/2)
codificaItaliana = dataISO[6:] + '-' + dataISO[4:6] + '-' + dataISO[:4]
rtMassimo = 3
Parte grafica comune e parte dipendente dalle scelte iniziali. Le variabili temporali date1
, date2
come le precedenti xValoriMediGiornalieri
e xMediaMobile
, tutte associate all'asse x, controllano la traslazione temporale dei valori in ordinata e sono scelte in modo da ottenere un accordo visivo soddisfacente con i dati giornalieri e in coerenza con il sito CovidStat.
# parte grafica comune
plt.rcParams['figure.figsize'] = [12, 6]
fig, ax = plt.subplots()
plt.ylim([0.4, rtMassimo])
y_major_ticks = np.arange(0.4, rtMassimo, 0.1)
ax.set_yticks(y_major_ticks)
ax.set_ylabel('indice di riproduzione Rt')
ax.grid(which='both', color='.85', linestyle='-', linewidth=1)
# parte grafica variabile
if fasce == "s": # assieme all'Rt si visualizzano le fasce di confidenza
ax.set_title(regioneScelta + ": Indice di riproduzione Rt e intervalli di confidenza 68%")
if unitaAssex == "m":
ax.set_xlabel('date dal 24 febbraio 2020')
date2 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)
for i in range(numGiorni-intervalloGiorni-12)]
ax.fill_between(date2, media14Giorni, media14GiorniSup, interpolate=True, facecolor='papayawhip')
ax.fill_between(date2, media14Giorni, media14GiorniInf, interpolate=True, facecolor='papayawhip')
ax.plot(date2, np.full(xMediaMobile.shape, 1), color='b')
ax.plot(date2, media14GiorniSup, color='lightskyblue', linewidth=2, label='media mobile Rt ($g+\Delta g$)')
ax.plot(date2, media14GiorniInf, color='greenyellow', linewidth=2, label='media mobile Rt ($g-\Delta g$)')
ax.plot(date2, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)
date1 = [datetime.date(2020, 3, 1) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)
for i in range(numGiorni-intervalloGiorni+1)]
ax.scatter(date1, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)
ax.text(datetime.date(2020, 11, 1), 0.8*rtMassimo, 'aggiornato il\n' + codificaItaliana)
else:
ax.set_xlabel('giorni dal 24 febbraio 2020')
ax.fill_between(xMediaMobile, media14Giorni, media14GiorniSup, interpolate=True, facecolor='papayawhip')
ax.fill_between(xMediaMobile, media14Giorni, media14GiorniInf, interpolate=True, facecolor='papayawhip')
ax.plot(xMediaMobile, np.full(xMediaMobile.shape, 1), color='b')
ax.plot(xMediaMobile, media14GiorniSup, color='lightskyblue', linewidth=2, label='media mobile Rt ($g+\Delta g$)')
ax.plot(xMediaMobile, media14GiorniInf, color='greenyellow', linewidth=2, label='media mobile Rt ($g-\Delta g$)')
ax.plot(xMediaMobile, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)
ax.scatter(xValoriMediGiornalieri, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)
ax.text(140, 0.8*rtMassimo, 'aggiornato il\n' + codificaItaliana)
x_major_ticks = np.arange(0, numGiorni+20, 20)
ax.set_xticks(x_major_ticks)
else: # visualizza solo l'andamento di Rt
ax.set_title(regioneScelta + ": Indice di riproduzione Rt")
if unitaAssex == "m":
ax.set_xlabel('date dal 24 febbraio 2020')
date2 = [datetime.date(2020, 3, 8) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)
for i in range(numGiorni-intervalloGiorni-12)]
ax.plot(date2, np.full(xMediaMobile.shape, 1), color='b')
ax.plot(date2, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)
date1 = [datetime.date(2020, 3, 1) + datetime.timedelta(days=int(intervalloGiorni/2)) + datetime.timedelta(days=i)
for i in range(numGiorni-intervalloGiorni+1)]
ax.scatter(date1, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)
ax.text(datetime.date(2020, 11, 1), 0.8*rtMassimo, 'aggiornato il\n' + codificaItaliana)
else:
ax.set_xlabel('giorni dal 24 febbraio 2020')
ax.plot(xMediaMobile, np.full(xMediaMobile.shape, 1), color='b')
ax.plot(xMediaMobile, media14Giorni, color='r', linewidth=2, label='media mobile Rt su 14 giorni', zorder=4)
ax.scatter(xValoriMediGiornalieri, valoriMediDayRtLin, s=4, label='Rt giornaliero', zorder=3)
ax.text(140, 0.8*rtMassimo, 'aggiornato il\n' + codificaItaliana)
x_major_ticks = np.arange(0, numGiorni+20, 20)
ax.set_xticks(x_major_ticks)
plt.legend()
plt.show()
[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)