[Introduction to SEIR model] Try fitting COVID-19 data ♬

I tried to fit the infection data of the new corona using the SEIR model. The method is the same as the previous [Introduction to minimize] Data analysis with SEIR model ♬. Infection data was obtained from the following reference sites. The following results were obtained in [Introduction to matplotlib] Reading the end time from COVID-19 data ♬ for Wuhan data. fig_Hubei_.png The initial value is selected and fitted so that such a curve can be obtained, and the result is obtained as follows. SEIR_Hubei_b0.00_ip0.00_gamma14.61_N112177_E0318_I0389_R00.02_.png So, this time, I would like to make similar fittings for each country and consider whether new knowledge can be obtained. The data is obtained from the following, and unless otherwise noted, it is ** data as of March 24, 2020 **, but the above is ** as of March 22, 2020 **. Also, as explained in the code explanation, please note that this time the fitting is performed for the infection data and the fitting using the healing data is executed by different programs. 【reference】 CSSEGISandData/COVID-19

What i did

・ The end of Japan ・ Code explanation ・ Fitting of national data

・ The end of Japan

I tried to write it gently as mentioned above, but since the Tokyo Metropolitan Government has said that it is a blockade of the city, I would like to mention whether the end of Japan can be seen first. From the conclusion, it seems that the end will surely be seen in another month. Actually, in the above data on March 22, there was a situation where I fell asleep for a while and this was the peak of infection, but since the number of infections increased linearly here, the saturation was slightly postponed. So, as a fitting, I got the following graph for the data up to March 24th. fig_Japan_.png

The fitting determines both the number of infections and the cure curve at the same time. The data for Japan is a little strange from the 53rd to the 59th, but it is generally moving naturally. So, I got a beautiful curve as follows. In other words, the result is that the peak appears in about 10 days, the healing curve peaks in about a month, and the rest gradually ends. SEIR_Japan_b1.59e-01_ip1.97e-23_gamma53_N1474_E01_I00_R04_.png

The obtained parameters are

beta_const=0.159
lp=1.97e-23
gamma=52.7
N=1474
R0 = 4.45

・ Code explanation

The entire code is below. The two reasons are that the provided data structure has changed and the current data needs to be read individually as in (2). In (2), the fitting function was changed from the normal maximum likelihood function to the normal square error so that the curve could be fitted more easily. ①collective_particles/fitting_COVID.pycollective_particles/fitting_COVID2.py

The simple code is explained below. Use the following Lib. It is the sum of the articles from the previous two articles.

#include package
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib.pyplot as plt
import pandas as pd

The following is the same as the infection data drawing.

#Read CSV data with pandas.
data = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv')
data_r = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Recovered.csv')
data_d = pd.read_csv('COVID-19/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Deaths.csv')
    
confirmed = [0] * (len(data.columns) - 4)
confirmed_r = [0] * (len(data_r.columns) - 4)
confirmed_d = [0] * (len(data_d.columns) - 4)
days_from_22_Jan_20 = np.arange(0, len(data.columns) - 4, 1)

#City
city = "Hubei" #Wuhan data

#Process the data
t_cases = 0

for i in range(0, len(data), 1):
    #if (data.iloc[i][1] == city): #for country/region
    if (data.iloc[i][0] == city):  #for province:/state  
        print(str(data.iloc[i][0]) + " of " + data.iloc[i][1])
        for day in range(4, len(data.columns), 1):
            confirmed[day - 4] += data.iloc[i][day] -  data_r.iloc[i][day] 
            confirmed_r[day - 4] += data_r.iloc[i][day]
            confirmed_d[day - 4] += data_d.iloc[i][day]

Use the SEIR model.

#define differencial equation of seir model
def seir_eq(v,t,beta,lp,ip,N0):
    N=N0 #int(26749*1) #58403*4
    a = -beta*v[0]*v[2]/N
    b = beta*v[0]*v[2]/N-(1/lp)*v[1] #N
    c = (1/lp)*v[1]-(1/ip)*v[2]
    d = (1/ip)*v[2]
    return [a,b,c,d]

Here, the picture obtained depends on the value of N0. It's good to fit this, but this time it's not done. When fitting Hubei, I referred to the following values. However, for ** N and S0, the actual number of infected people (values moved around) was adopted based on the idea that the residents are not all infected in the sense that they are participants in the game. 【reference】 ・ I tried to predict the behavior of the new coronavirus with the SEIR model. Initial value setting

#solve seir model
N0 = 70939 #Hubei
N,S0,E0,I0=int(N0*1.),N0,int(318),int(389) #N is total population, S0 initial value of fresh peaple
ini_state=[S0,E0,I0,0] #initial value of diff. eq.
beta,lp,ip=1, 2, 7.4 
t_max=60 #len(days_from_22_Jan_20)
dt=0.01
t=np.arange(0,t_max,dt)
obs_i = confirmed

A function that returns an evaluation value for an input parameter. The following is v [2]; returns the number of infections.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip,N):
    v=odeint(seir_eq,ini_state,t,args=(beta,lp,ip,N))
    est=v[0:int(t_max/dt):int(1/dt)]
    return est[:,2]

The objective function of minimization below (the following uses the maximum likelihood function).

#define logscale likelihood function
def y(params):
    est_i=estimate_i(ini_state,params[0],params[1],params[2],params[3])
    return np.sum(est_i-obs_i*np.log(np.abs(est_i)))

Find the minimum value with minimize @ scipy and calculate the basic reproduction number.

#optimize logscale likelihood function
mnmz=minimize(y,[beta,lp,ip,N],method="nelder-mead")
print(mnmz)
#R0
beta_const,lp,gamma_const,N = mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3] #Infection rate, infection waiting time, removal rate (recovery rate)
print(beta_const,lp,gamma_const,N)
R0 = N*beta_const*(1/gamma_const)
print(R0)

Graph display of calculation results.

#plot reult with observed data
fig, (ax1,ax2) = plt.subplots(2,1,figsize=(1.6180 * 4, 4*2))
lns1=ax1.plot(obs_i,"o", color="red",label = "data")
lns2=ax1.plot(confirmed_r,"*", color="green",label = "data_r")
lns3=ax1.plot(estimate_i(ini_state,mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3]), label = "estimation")
lns0=ax1.plot(t,odeint(seir_eq,ini_state,t,args=(mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3])))
lns_ax1 = lns1+lns2 + lns3 + lns0
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)
ax1.set_ylim(0,N0)

The following displays a long-span graph so that the end can be seen.

t_max=200 #len(days_from_22_Jan_20)
t=np.arange(0,t_max,dt)

lns4=ax2.plot(obs_i,"o", color="red",label = "data")
lns5=ax2.plot(confirmed_r,"*", color="green",label = "recovered")
lns6=ax2.plot(t,odeint(seir_eq,ini_state,t,args=(mnmz.x[0],mnmz.x[1],mnmz.x[2],mnmz.x[3])))
ax2.legend(['data','data_r','Susceptible','Exposed','Infected','Recovered'], loc=1)
ax2.set_title('SEIR_b{:.2f}_ip{:.2f}_gamma{:.2f}_N{:.2f}_E0{:d}_I0{:d}_R0{:.2f}'.format(beta_const,lp,gamma_const,N,E0,I0,R0))
ax1.grid()
ax2.grid()
plt.savefig('./fig/SEIR_{}_b{:.2f}_ip{:.2f}_gamma{:.2f}_N{:.2f}_E0{:d}_I0{:d}_R0{:.2f}_.png'.format(city,beta_const,lp,gamma_const,N,E0,I0,R0)) 
plt.show()
plt.close()

On the other hand, in the simultaneous fitting of the number of infections and the number of cures, the evaluation function and objective function in the middle were changed as follows.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip,N):
    v=odeint(seir_eq,ini_state,t,args=(beta,lp,ip,N))
    est=v[0:int(t_max/dt):int(1/dt)] 
    return est[:,2],est[:,3] #v0-S,v1-E,v2-I,v3-R
    
#define logscale likelihood function
def y(params):
    est_i_2,est_i_3=estimate_i(ini_state,params[0],params[1],params[2],params[3])
    return np.sum((est_i_2-obs_i_2)*(est_i_2-obs_i_2)+(est_i_3-obs_i_3)*(est_i_3-obs_i_3))

In this case, it is possible to fit only the number of infections or only the number of cures.

・ Fitting of national data

·Iran Iran, which has a relatively high cure rate, saturates faster than Japan and is likely to end, but the number of infections is still increasing and the cure rate is sluggish, so it may take longer. fig_Iran_.png SEIR_Iran_b1.60e-01_ip1.35e-23_gamma17_N38098_E01_I00_R0350_.png

·Korea Somehow, the distortion became more noticeable compared to the previous data, and the credibility of the data was broken. Still, the cure rate is increasing overall, and the end seems to be relatively early. fig_Korea, South_.png SEIR_Korea, South_b1.98e-01_ip5.42e-32_gamma56_N16115_E01_I00_R057_.png

·Italy Although the mortality rate is extremely high and dangerous, the data are very clean and sophisticated. Therefore, it was easy to fit and I got on cleanly. The healing rate is still about 10%, and it seems that it will take some time (healing peak is 2 months), but if the peak is seen after this, the end is near. fig_Italy_.png SEIR_Italy_b1.61e-01_ip6.22e-26_gamma78_N116637_E01_I00_R0242_.png

The following is the initial expansion period, and it seems that the fitting accuracy is poor and it is not reliable, but I will post it. ·Spain It stands up more and more, but it's difficult to match. fig_Spain_.png SEIR_Spain_b1.50e-01_ip1.47e-30_gamma69_N104007_E01_I00_R0226_.png ·Germany This has not been completed yet, and the number of cures is too small in the first place, so it seems that it is in the initial expansion period. fig_Germany_.png SEIR_Germany_b1.22e-01_ip1.51e-37_gamma160_N29502_E01_I00_R023_.png

Summary

・ I tried fitting COVID-19 data ・ Japan, Iran, and South Korea are expected to see an end following China. ・ Other Europe and the United States are still in the expansion period, and it is difficult to predict the end.

・ I talked roughly this time, but I would like to make a more detailed analysis while watching the daily updates.

bonus

Simultaneous fitting of Wuhan is as follows Wuhan data does not seem to be explained by this model (> <) raw data fig_Hubei_.png Infection data only fitting SEIR_Hubei_b5.47e-01_ip4.97e+00_gamma15_N23774_E0318_I0389_R0866_.png Healing data only fitting SEIR_Hubei_b2.11e-01_ip1.15e+00_gamma6_N176229_E0318_I0389_R06739_.png Fitting infection and cure data SEIR_Hubei_b3.68e+00_ip2.88e+01_gamma23_N-1287313_E0318_I0389_R0-204877_.png

Recommended Posts

[Introduction to SEIR model] Try fitting COVID-19 data ♬
[Introduction to minimize] Data analysis with SEIR model ♬
[Introduction to SIR model] Predict the end time of each country with COVID-19 data fitting ♬
Fitting to ARMA, ARIMA model
[Introduction to matplotlib] Read the end time from COVID-19 data ♬
[Introduction to infectious disease model] I tried fitting and playing ♬
[Introduction to Tensorflow] Understand Tensorflow properly and try to make a model
Introduction to Statistical Modeling for Data Analysis GLM Model Selection
Try to put data in MongoDB
[Introduction to SIR model] Consider the fitting result of Diamond Princess ♬
[Introduction to Python3, Day 17] Chapter 8 Data Destinations (8.1-8.2.5)
[Introduction to Python3, Day 17] Chapter 8 Data Destinations (8.3-8.3.6.1)
[Introduction to Python3 Day 19] Chapter 8 Data Destinations (8.4-8.5)
[Introduction to Python3 Day 18] Chapter 8 Data Destinations (8.3.6.2 to 8.3.6.3)
Try converting to tidy data with pandas
[Introduction to Data Scientists] Basics of Python ♬
Try using django-import-export to add csv data to django
Try to aggregate doujin music data with pandas
An introduction to statistical modeling for data analysis
[Introduction to Python] How to handle JSON format data
Preparing to try "Data Science 100 Knock (Structured Data Processing)"
[Introduction to cx_Oracle] (5th) Handling of Japanese data
Introduction to MQTT (Introduction)
Introduction to Scrapy (1)
Introduction to Supervisor
Introduction to Tkinter 1: Introduction
Introduction to PyQt
Introduction to Scrapy (2)
[Linux] Introduction to Linux
Introduction to Scrapy (4)
Introduction to discord.py (2)
Introduction to discord.py
Reading Note: An Introduction to Data Analysis with Python
Try to divide twitter data into SPAM and HAM
[CovsirPhy] COVID-19 Python package for data analysis: SIR-F model
[CovsirPhy] COVID-19 Python Package for Data Analysis: SIR model
Try to decipher the login data stored in Firefox
Try writing JSON format data to object storage Cloudian/S3
[Technical book] Introduction to data analysis using Python -1 Chapter Introduction-