Per determinare il tasso di crescita λ viene innanzitutto linearizzata la relazione esponenziale che lo coinvolge nel modello sviluppato. Su tale base si calcola il suo valore giornaliero e per rirdurre i fattori di disturbo, si esegue una media mobile di 14 giorni.
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
.
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))
######### 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 richiede inoltre l'unità temporale da porre in ascissa.
print()
print("Si intendono elaborare i nuovi positivi nazionali (i) o quelli di una regione (r)? ")
scelta = input(" :> ")
print("Inserire la data nel formato (YYYMMDD) ")
dataISO = input(" :> ")
print("Si vuole l'asse temporale espresso in mesi (m) o in giorni trascorsi (g) dall'inizio pandemia? ")
unitaAssex = input(" :> ")
Vengono composti i nomi dei file e, una volta letti i dati contenuti, questi vengono caricano 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)
Per ridurre i fattori giornalieri o settimanali di distrurbo, si esegue una media mobile su un intervallo di 14 giorni dei valori giornalieri di λ.
media14Giorni = np.mean(partiziona(valoriLambdaLin, intervalloGiorni), axis=1)
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 comune
codificaItaliana = dataISO[6:] + '-' + dataISO[4:6] + '-' + dataISO[:4]
plt.rcParams['figure.figsize'] = [12, 6]
fig, ax = plt.subplots()
ax.grid(which='both', color='.85', linestyle='-', linewidth=1)
ax.set_ylabel('tasso di crescita λ (1/giorni)')
ax.set_title(regioneScelta + ": tasso di crescita giornaliero λ (regressione lineare)\ne media mobile su 14 giorni")
# parte variabile
if unitaAssex == "m":
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, valoriLambdaLin, s=4, label='λ giornaliero', zorder=3)
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, media14Giorni, color='tomato', linewidth=2, label='media mobile λ su 14 giorni', zorder=4)
ax.set_xlabel('date dal 24 febbraio 2020')
ax.text(datetime.date(2020, 5, 13), 0.16, 'aggiornato il\n' + codificaItaliana)
else:
xValoriMediGiornalieri = np.arange(intervalloGiorni-1, numGiorni)
xMediaMobile = np.arange(intervalloGiorni+5, numGiorni-intervalloGiorni/2)
ax.scatter(xValoriMediGiornalieri, valoriLambdaLin, s=4, label='λ giornaliero', zorder=3)
ax.plot(xMediaMobile, media14Giorni, color='tomato', linewidth=2, label='media mobile λ su 14 giorni', zorder=4)
ax.text(85, 0.16, 'aggiornato il\n' + codificaItaliana)
ax.set_xlabel('giorni dal 24 febbraio 2020')
x_major_ticks = np.arange(0, numGiorni+20, 20)
ax.set_xticks(x_major_ticks)
plt.legend()
plt.show()