[PYTHON] [Introduction to minimize] Data analysis with SEIR model ♬

In preparation for data analysis of COVID-19, we will explain how to output the following graph. This figure is a figure in which the SEIR model of infectious disease is fitted to the given data using minimize @ scipy. SEIR_b6.88_ip1.22_gamma2.01_N762_E01_I00_R03.41_.png The method is almost as shown below. Both were helpful in many ways. 【reference】 (1) Beginning of mathematical model of infectious disease: Outline of SEIR model by Python and introduction to parameter estimationMathematical prediction model for infectious diseases (SIR model): Case study (1)I tried to predict the behavior of the new coronavirus with the SEIR model.

What i did

・ Code explanation ・ Try fitting with SIR model and SEIR model ・ Other models

・ Code explanation

The entire code is placed below. -Collive_particles / minimize_params.py I will leave the explanation to the above reference, and here I will explain the code used this time. The data was used from Reference (2). Also, for the minimize code, I referred to Reference ①. The Lib to be used is as follows. This time, it is implemented in the keras-gpu environment of the conda environment of Windows 10, but it is previously installed including scipy.

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

The differential equations of the SEIR model are defined below. Here, what to do with N will be decided as appropriate, but in this proposition, N = 763 people from the paper published in Reference (2), which was fixed. The differential equations of the SEIR model are as follows Here, in Reference (1), N is pushed into the variable beta, which is not explicitly expressed, but the following formula is easier to understand in terms of Uwan. The following is quoted from Reference ③.

\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dE}{dt} &=  \beta \frac{SI}{N}  -\epsilon E \\
\frac{dI}{dt} &=  \epsilon E -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

If you look at the following while looking at the above, you can understand it because it is just right.

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

First, a solution is obtained based on the initial value as an ordinary differential equation.

#solve seir model
beta,lp,ip=6.87636378, 1.21965986, 2.01373496 #2.493913  , 0.95107715, 1.55007883
plt.plot(t,odeint(seir_eq,ini_state,t,args=(beta,lp,ip))) #0.0001,1,3

The result is as follows. If the parameters are estimated and changed here, most fittings can be performed, but complicated ones are difficult. SEIR_calc_.png As shown in Reference (2), substitute the data and output it to the graph as shown below.

#show observed i
data_day = [1,2,3,4,5,6,7,8,9,10,11,12,13,14]
obs_i = data_influ

plt.plot(obs_i,"o", color="red",label = "data")

The result is as follows, and since this is an output for confirmation on the way, headings etc. are omitted. SEIR_data_.png Now, the next is the estimation part. This evaluation function calculates with a differential equation using the value of the argument and returns the evaluation result.

#function which estimate i from seir model func 
def estimate_i(ini_state,beta,lp,ip):
    return est[:,2]

The log-likelihood function is defined as follows.

#define logscale likelihood function
def y(params):
    return np.sum(est_i-obs_i*np.log(np.abs(est_i)))

Use the following Scipy minimize function to find the maximum value.

#optimize logscale likelihood function

The output is as follows

>python minimize_params.py
 final_simplex: (array([[6.87640764, 1.21966435, 2.01373196],
       [6.87636378, 1.21965986, 2.01373496],
       [6.87638203, 1.2196629 , 2.01372646],
       [6.87631916, 1.21964456, 2.0137297 ]]), array([-6429.40676091, -6429.40676091, -6429.40676091, -6429.4067609 ]))
           fun: -6429.406760912483
       message: 'Optimization terminated successfully.'
          nfev: 91
           nit: 49
        status: 0
       success: True
             x: array([6.87640764, 1.21966435, 2.01373196])

Take out beta_const, lp, gamma_const from the calculation result as follows, and calculate R0 below. This is the basic reproduction number, which is a guideline for the spread of infection. If R0 <1, it does not spread, but if R0> 1, it spreads and is an index showing the degree.

beta_const,lp,gamma_const = mnmz.x[0],mnmz.x[1],mnmz.x[2] #Infection rate, infection waiting time, removal rate (recovery rate)
R0 = beta_const*(1/gamma_const)

This time, it is as follows, and since R0 = 3.41, the infection has spread.

6.876407637532918 1.2196643491443309 2.0137319643699927

The obtained result is displayed on the graph.

#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(estimate_i(ini_state,mnmz.x[0],mnmz.x[1],mnmz.x[2]), label = "estimation")
lns_ax1 = lns1+lns2
labs_ax1 = [l.get_label() for l in lns_ax1]
ax1.legend(lns_ax1, labs_ax1, loc=0)

lns3=ax2.plot(obs_i,"o", color="red",label = "data")
ax2.legend(['data','Susceptible','Exposed','Infected','Recovered'], loc=0)

The result is as above.

・ Other models

When I was doing this calculation, I found a slightly simpler SIR model and a slightly more complicated model, so I will summarize them. The SIR model is as follows because it is a model that ignores Exposed from S-E-I-R.

\frac{dS}{dt} &= -\beta \frac{SI}{N} \\
\frac{dI}{dt} &=  \beta \frac{SI}{N} -\gamma I \\
\frac{dR}{dt} &=  \gamma I \\

You can understand the meaning of $ R_0 $ by rewriting as follows.

\frac{dS}{dt} &= -R_0 \frac{S}{N} \frac{dR}{dt} \\
R_0 &= \frac{\beta}{\gamma}\\ 
\frac{dI}{dt} &= 0 \Infected person peak at\\
R_0 &=  \frac{\beta}{\gamma}=  \frac{N}{S_0} \\

On the other hand, in a complicated situation such as a transition of death and death without notification, it will be as follows. It is assumed that deaths do not contribute to infection.


-Data fitting can now be done automatically by minimize @ scipy. ・ I was able to fit the actual data with the SEIR model.

・ Next time, I would like to analyze the data of each country of COVID-19.

Recommended Posts

[Introduction to minimize] Data analysis with SEIR model ♬
[Introduction to SEIR model] Try fitting COVID-19 data ♬
Reading Note: An Introduction to Data Analysis with Python
20200329_Introduction to Data Analysis with Python Second Edition Personal Summary
Introduction to Statistical Modeling for Data Analysis GLM Model Selection
Introduction to Data Analysis with Python P32-P43 [ch02 3.US Baby Names 1880-2010]
Introduction to Data Analysis with Python P17-P26 [ch02 1.usa.gov data from bit.ly]
Data analysis with python 2
I tried fMRI data analysis with python (Introduction to brain information decoding)
Data analysis with Python
[Technical book] Introduction to data analysis using Python -1 Chapter Introduction-
Introduction to RDB with sqlalchemy Ⅰ
Introduction to Machine Learning with scikit-learn-From data acquisition to parameter optimization
[Introduction to Python] How to get data with the listdir function
How to deal with imbalanced data
Introduction to RDB with sqlalchemy II
How to Data Augmentation with PyTorch
Data analysis starting with python (data visualization 1)
Introduction to image analysis opencv python
Data analysis starting with python (data visualization 2)
[Stock price analysis] Learn pandas with Nikkei 225 (004: Change read data to Nikkei 225)
Introduction to Structural Equation Modeling (SEM), Covariance Structure Analysis with Python
Links to people who are just starting data analysis with python
Introduction to Statistical Modeling for Data Analysis Generalized Linear Models (GLM)
From the introduction of JUMAN ++ to morphological analysis of Japanese with Python
"Introduction to data analysis by Bayesian statistical modeling starting with R and Stan" implemented in Python
[Introduction to WordCloud] Let's play with scraping ♬
Introduction to Statistical Modeling for Data Analysis GLM Likelihood-Ratio Test and Test Asymmetry
[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-
Introduction to Python Image Inflating Image inflating with ImageDataGenerator
I tried factor analysis with Titanic data!
Convert Excel data to JSON with python
[Introduction to Python] Let's use foreach with Python
[Introduction to Python3 Day 19] Chapter 8 Data Destinations (8.4-8.5)
[Python] Introduction to CNN with Pytorch MNIST
Send data to DRF API with Vue.js
Convert FX 1-minute data to 5-minute data with Python
Try converting to tidy data with pandas
[Introduction to Pytorch] I played with sinGAN ♬
Data analysis starting with python (data preprocessing-machine learning)
How to read problem data with paiza
I tried to predict the behavior of the new coronavirus with the SEIR model.
[Introduction to Data Scientists] Basics of Python ♬
Introduction to Statistical Modeling for Data Analysis Expanding the range of applications of GLM
How to replace with Pandas DataFrame, which is useful for data analysis (easy)
[In-Database Python Analysis Tutorial with SQL Server 2017] Step 2: Import data to SQL Server using PowerShell
Introduction to Time Series Analysis ~ Seasonal Adjustment Model ~ Implemented in R and Python
An introduction to data analysis using Python-To increase the number of video views-
[Introduction to Python] How to get the index of data with a for statement
How to create sample CSV data with hypothesis
Markov Chain Chatbot with Python + Janome (1) Introduction to Janome
Markov Chain Chatbot with Python + Janome (2) Introduction to Markov Chain
Try to aggregate doujin music data with pandas
Convert data with shape (number of data, 1) to (number of data,) with numpy.
I tried to save the data with discord
[Introduction to StyleGAN2] Independent learning with 10 anime faces ♬
Introduction to Tornado (1): Python web framework started with Tornado
I tried principal component analysis with Titanic data!
Save data to flash with STM32 Nucleo Board
I tried to get CloudWatch data with Python