確率的な状態変化をコード上で実現したい。

実現したいこと

現在、SIRモデルと呼ばれるモデルの拡張したモデルを以下の通り計算しています。

以下のコードに関して、状態変化（SからE、EからP1など）を確率的に行うことを目指しています。

しかし、tやNの扱い方がうまくわからず苦戦しています。

どなたかお知恵をお書いいただけないでしょうか。

例えば1試行に関して、Eを確率σでdt秒あたりにP1に変化させる、といったことを実現したいです。

発生している問題・エラーメッセージ

N_やtの定義の文（37、38行目）でエラーが出て苦戦しております。

該当のソースコード

python

1import numpy as np 2import matplotlib.pyplot as plt 3 4#各パラメータの設定5beta=4/9 #曝露率6sigma=1/3 #感染率7rho1=1 #P1からP2になる確率8rho2=1 #P2からIaかIsになる確率9eta=0.54 #PからIsになる確率10gamma=1/7 #IsまたはIaからRになる確率11R_0=412N=100 #全人口13S0=9914E0=115P1_0=016P2_0=017Ia_0=018Is_0=019 20dt=0.001 #刻み幅21t=022S=S0 23E=E0 24P1=P1_0 25P2=P2_0 26Ia=Ia_0 27Is=Is_0 28SList=[S0]29EList=[E0]30P1List=[P1_0]31P2List=[P2_0]32IaList=[Ia_0]33IsList=[Is_0]34tList=[t]35 36for i in range(30000):37N_=S+E+P1+P2+Ia+Is 38t=dt*(i+1)39P1_ = np.random.binomial(P1, 1 - rho1*dt) # P1からP2に遷移する40P2_ = np.random.binomial(P2, 1 - rho2*dt) # P2からIaまたはIsに遷移する41P1_ += np.random.binomial(S, beta*dt*(P1+P2+Ia+Is)/N_) # SからP1に遷移する42P2_ += np.random.binomial(S, beta*dt*(P1+P2+Ia+Is)/N_) # SからP2に遷移する43E_ = np.random.binomial(E, 1 - sigma*dt) # EからIaまたはIsに遷移する44E_ += np.random.binomial(S, beta*dt*(P1+P2+Ia+Is)/N_) # SからEに遷移する45Ia_ = np.random.binomial(Ia, 1 - gamma*dt) # IaからRに遷移する46Is_ = np.random.binomial(Is, 1 - gamma*dt) # IsからRに遷移する47Ia_ += np.random.binomial(P2, (1-eta)*rho2*dt) # P2からIaに遷移する48Is_ += np.random.binomial(P2, eta*rho2*dt) # P2からIsに遷移する49S_ = N_ - E_ - P1_ - P2_ - Ia_ - Is_ # Sの更新50S=S_ 51E=E_ 52P1=P1_ 53P2=P2_ 54Ia=Ia_ 55Is=Is_ 56 tList.append(t)57 SList.append(S)58 EList.append(E)59 P1List.append(P1)60 P2List.append(P2)61 IaList.append(Ia)62 IsList.append(Is)63plt.plot(tList,SList,color='green',label='S')64plt.plot(tList,EList,color='red',label='E')65plt.plot(tList,P1List,color='orange',label='P1')66plt.plot(tList,P2List,color='blue',label='P2')67plt.plot(tList,IaList,color='skyblue',label='Ia')68plt.plot(tList,IsList,color='purple',label='Is')69plt.xlim(0,20)70plt.ylim(0,1)71plt.legend()72plt.show()