간단한 SIQR 모델 표시

#! /usr/bin/python
#
#   siqr01.py
#
#                       May/14/2020
#
# ------------------------------------------------------------------
import matplotlib.pyplot as plt
import numpy as np
import  sys

# ------------------------------------------------------------------
def plot_graph(period,r0,p_kakuri,l_ss,l_ii,l_qq,l_rr):
    plt.rcParams["font.family"] = "TakaoExGothic"
    # plt.rcParams["font.family"] = "IPAGothic"

    plt.figure()
    str_out = "基本再生産数 : {}".format(format(r0,".3f"))
    str_out += "  "
    str_out += "隔離率 : {}".format(format(p_kakuri,".3f"))
    plt.title(str_out)
    plt.xlabel('時間')
    plt.ylabel('比率')

    x = np.linspace(0, 100, period+1)
    plt.plot(x, l_ss, label="未感染者")
    plt.plot(x, l_ii, label="感染 非隔離者")
    plt.plot(x, l_qq, label="感染 隔離者")
    plt.plot(x, l_rr, label="回復者")
    plt.legend()
    plt.show()
# ------------------------------------------------------------------
def next_state(beta,nu,p_kakuri,ss, ii, qq, rr):
    q_kakuri = 0.0
    beta_ss_ii = beta * ss * ii
    delta_ss = - beta_ss_ii
    delta_ii = (1 - qq) * beta_ss_ii - p_kakuri * ii - nu * ii
    delta_qq = q_kakuri * beta_ss_ii + p_kakuri * ii - nu * qq
    delta_rr = nu * ii + nu * qq
    ss = ss + delta_ss
    ii = ii + delta_ii
    qq = qq + delta_qq
    if ss < 0:
        ss = 0
    if 1.0 < ii:
        ii = 1.0
    rr = 1.0 - ss - ii - qq
#
    return ss, ii, qq, rr
#
# ------------------------------------------------------------------
def calc_proc(beta,nu,p_kakuri,period):
#
    ii = 0.01
    qq = 0
    rr = 0
    ss = 1 - ii - rr

    l_ss = [ss]
    l_ii = [ii]
    l_qq = [qq]
    l_rr = [rr]
    r0 = beta*ss/nu

    for t in range(period):
        ss, ii, qq, rr = next_state(beta,nu,p_kakuri,ss, ii, qq, rr)
        l_ss.append(ss)
        l_ii.append(ii)
        l_qq.append(qq)
        l_rr.append(rr)
#
    return r0,l_ss,l_ii,l_qq,l_rr
#
# ------------------------------------------------------------------
sys.stderr.write("*** start ***\n")
beta = 0.5
nu = 0.3
p_kakuri = 0.0
period = 100
#
sys.stderr.write("beta = %f\n" % beta)
sys.stderr.write("nu = %f\n" % nu)
sys.stderr.write("p_kakuri = %f\n" % p_kakuri)
#
r0,l_ss,l_ii,l_qq,l_rr = calc_proc(beta,nu,p_kakuri,period)
plot_graph(period,r0,p_kakuri,l_ss,l_ii,l_qq,l_rr)
#
p_kakuri = 0.1
sys.stderr.write("p_kakuri = %f\n" % p_kakuri)
r0,l_ss,l_ii,l_qq,l_rr = calc_proc(beta,nu,p_kakuri,period)
plot_graph(period,r0,p_kakuri,l_ss,l_ii,l_qq,l_rr)
#
sys.stderr.write("*** end ***\n")
# ------------------------------------------------------------------