SIQR 모델 표시

#! /usr/bin/python
#
#   siqr02.py
#
#                       May/14/2020
# ------------------------------------------------------------------
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import  sys

# ------------------------------------------------------------------
def SIQR_EQ(v, t, beta, gamma, p_kakuri):
    q_kakuri = 0.0
    rvalue = [0,0,0,0]
    beta_ss_ii = beta * v[0] * v[1]
    rvalue[0] = - beta_ss_ii
    rvalue[1] = (1 - q_kakuri) * beta_ss_ii - p_kakuri * v[1] - gamma * v[1]
    rvalue[2] = q_kakuri * beta_ss_ii + p_kakuri * v[1] - gamma * v[2]
    rvalue[3] = gamma * v[1] + gamma * v[2]
#   rvalue[3] = -(rvalue[0] + rvalue[1] + rvalue[2])
#
    return rvalue
#
# ------------------------------------------------------------------
def calc_proc(times,beta_const,gamma_const,p_kakuri):
    S_0=999
    I_0=1
    Q_0=0
    R_0=0
    ini_state = [S_0,I_0,Q_0,R_0]

    args  =(beta_const, gamma_const, p_kakuri)

    N_total = S_0 + I_0 + Q_0 + R_0
    R0 = N_total*beta_const *(1/gamma_const)
    print(R0)

#Numerical Solution using scipy.integrate
#Solver SIR model
    result = odeint(SIQR_EQ, ini_state, times, args)
#
    return R0,result
# ------------------------------------------------------------------
def plot_graph_proc(R0,p_kakuri,times,result):
    plt.rcParams["font.family"] = "TakaoExGothic"
#    plt.rcParams["font.family"] = "IPAGothic"

    str_out = "基本再生産数 : {}".format(format(R0,".3f"))
    str_out += "  "
    str_out += "隔離率 : {}".format(format(p_kakuri,".3f"))
    plt.title(str_out)
    plt.xlabel('日数')
    plt.ylabel('人数')
    plt.plot(times,result)
    plt.legend(['未感染者','感染 非隔離者','感染 隔離者','回復者'])
    plt.show()
#
# ------------------------------------------------------------------
sys.stderr.write("*** start ***\n")
t_max = 160
dt = 0.01
#
beta_const = 0.2/1000
gamma_const = 0.1
p_kakuri = float(sys.argv[1])
sys.stderr.write("p_kakuri = %.3f\n" % p_kakuri)

times =np.arange(0,t_max, dt)
R0,result = calc_proc(times,beta_const,gamma_const,p_kakuri)
#
plot_graph_proc(R0,p_kakuri,times,result)
#
sys.stderr.write("*** end ***\n")
# ------------------------------------------------------------------

./siqr02.py 0.0

./siqr02.py 0.0