質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

88.93%

sirsモデル コロナに対するシミュレーションについて

解決済

回答 1

投稿 編集

  • 評価
  • クリップ 1
  • VIEW 144

itta0602

score 1

こちらはsirモデルについて書いたものをsirsモデルに変えようとしています。
参考サイトを貼ります。
[リンク内容](https://rpubs.com/ktgrstsh/tokyor84   )
こちらを元に助けてもらいたいです。
sirは実際に動いたのでこれにこちらのサイトでいうrhoに関するものを付け加えていくと言った感じです。
初めて質問させてもらうので色々不備あると思いますが、よろしくお願い致します。
質問してください!

from matplotlib import rcParams​
import matplotlib.pyplot as plt​
from scipy. integrate import odeint​
​
​
#sirモデルの常微分方程式​def SIR_equation(v,t,β,γ,):​
    S, I, R =v​
    dSdt=-β*S*I
    dIdt=β*S*I-γ*I​
    dRdt=γ*I
    return [dSdt,dIdt,dRdt]​
for β in np.linspace(0.10,0.20,11):​
    for γ in np.linspace(0.10,0.20,11):​

        ​
#初期値設定​
          S0=.999999928​
          I0=.000000071​
          R0=0.0#係数決定​
​
​
#0 ≦t ≦100のあいだを1000分割​
        time_list=np.linspace(1,150,150)​
​
        var_list=odeint(​
           SIR_equation,​
           [S0,I0,R0],​
           time_list,​
           args=(β,γ)​
)​
#微分方程式をここまでで解いた​
#var_listにはS、I、Rのそれぞれの時間に対する結果が配列で入っているので、​
#それらを分けて配列(_1ist)に収納​
        S_list = var_list [:,0] ​
        I_list = var_list [:,1] ​
        R_list = var_list [:,2] ​
        tot_list = S_list + I_list + R_list​
​
#精題にプロットできるよう色々工夫してみる​
       ​
        plt.xlabel("time")​
        plt.xlim(0,100)​
        plt.ylim(0,1)​
        plt.ylabel("population ratio")​
        plt.fill_between(time_list, tot_list,label ="population",alpha = 0.5)​
        #plt.fill_between(time_list, S_list,label = "Susceptible",alpha = 0.5)​
        #plt.fill_between(time_list, R_list,label = "Recovered",alpha = 0.5)​
        plt.fill_between(time_list, I_list,label="Infected",alpha = 0.5)​
        plt.fill_between(time_list, I_list+R_list,label="Infected(cumulative)",alpha = 0.5)​
        ​
​
        plt.legend(loc ='center right')​
        plt.title(r"SIR model ($\beta=%.2f,\gamma=%.2f$)"%(β,γ))​
        plt.savefig(r"SIR_model_beta=%.2f_gamma=%.2f.png"%(β,γ))​
        plt.close()​

コード
  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • itta0602

    2020/07/14 18:48

    ありがとうございました

    キャンセル

  • Penpen7

    2020/07/14 19:00 編集

    dSdt=-β*S*I
    dIdt=β*S*I-γ*I​
    dRdt=γ*I
    にρR項を文字通り付け足せばいいだけだと思うのですが、違うんですか?

    キャンセル

  • Penpen7

    2020/07/14 19:07

    def SIR_equation(v,t,β,γ,):

    def SIR_equation(v,t,β,γ,ρ):
    としたかったのですかね?

    キャンセル

回答 1

checkベストアンサー

+1

「文字通り」ρを加えればいいだけの話でしょう。
(各定数は参考URLに準拠, for文は簡単化のため削除)

from matplotlib import rcParams
import matplotlib.pyplot as plt
from scipy. integrate import odeint
import numpy as np


def SIRS_equation(v, t, β, γ, ρ):
    S, I, R = v
    dSdt = -β*S*I+ρ*R
    dIdt = β*S*I-γ*I
    dRdt = γ*I-ρ*R
    return [dSdt, dIdt, dRdt]

ρ=0.01
β=1.0
γ=0.1
S0 = .9999999
I0 = 1-S0
R0 = 0.0
time_list = np.linspace(0, 150, 150)

var_list = odeint(
    SIRS_equation, [S0, I0, R0], time_list, args=(β, γ, ρ))
S_list = var_list[:, 0]
I_list = var_list[:, 1]
R_list = var_list[:, 2]
tot_list = S_list + I_list + R_list

# 精題にプロットできるよう色々工夫してみる

plt.xlabel("time")
plt.xlim(0, 100)
plt.ylim(0, 1)
plt.ylabel("population ratio")
plt.fill_between(time_list, tot_list, label="population", alpha=0.5)
plt.fill_between(time_list, I_list, label="Infected", alpha=0.5)
plt.fill_between(time_list, I_list+R_list, label="Infected(cumulative)", alpha=0.5)
plt.legend(loc ='center right')
plt.title(r"SIRS model ($\beta=%.2f,\gamma=%.2f, \rho=%.2f$)"%(β, γ, ρ))
plt.savefig(r"SIRS_model_beta=%.2f_gamma=%.2f_rho=%.2f.png" %(β, γ, ρ))
plt.close()


イメージ説明

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/14 20:49

    ありがとうございます!
    とても助かりました!
    for 文も同様にして入れたらいいですよね??

    キャンセル

  • 2020/07/15 01:15

    ρ=0.01
    β=1.0
    γ=0.1

    for ρ
    for β
    for γ
    に置き換えればいいでしょう

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 88.93%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る