確率的なシミュレーションにおいて、設定値と異なる初期値がグラフに現れる

実現したいこと

感染症モデルであるSIRモデルを確率的にシミュレーションすることを目指しています。

具体的には、S状態の人1人について、Δt=0.01の間にβIΔtの確率でIに変化させる、また、I状態の人1人について、γΔtの確率でRに変化させる操作をS状態の人、I状態の人の数だけ確率的にシミュレーションし、これをt=1まで行う操作をまず10回繰り返してそのグラフの平均を取り、得られたt=1での値をt=1からt=2までのシミュレーションの初期条件とする、という操作をt=13まで繰り返し、その時間発展のグラフを得るということを目標にしています。

しかしながら、実際に描画した場合、初期条件をS=99,I=1,R=0で設定しているにも関わらず、t＝0でR=99になってしまうという不具合が発生しています。

Pythonに精通されている方、どなたかこの解決策についてご教授いただけないでしょうか。

また、S状態の人一人一人に関して確率的な状態変化を行うことを、19行目のfor文で行うという発想は正しいでしょうか。併せてご教授いただければと思います。

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

グラフを描画した際、設定した初期条件と異なる初期条件が現れる。

該当のソースコード

Python

1import numpy as np 2import matplotlib.pyplot as plt 3 4def sir_simulation_stochastic(S0, I0, R0, beta, gamma, t_max, dt, num_simulations):5 t_values = np.arange(0, t_max + dt, dt)6 S_values = []7 I_values = []8 R_values = []9 10 for _ in range(num_simulations):11 S = S0 12 I = I0 13 R = R0 14 S_temp = [S0]15 I_temp = [I0]16 R_temp = [R0]17 18 for t in t_values:19 for _ in range(S):20 prob_infection_s = beta * I * dt 21 if np.random.random() < prob_infection_s:22 S -= 123 I += 124 25 for _ in range(I):26 prob_recovery_i = gamma * dt 27 if np.random.random() < prob_recovery_i:28 I -= 129 R += 130 31 S_temp.append(S)32 I_temp.append(I)33 R_temp.append(R)34 35 S_values.append(S_temp[:len(t_values)])36 I_values.append(I_temp[:len(t_values)])37 R_values.append(R_temp[:len(t_values)])38 S0, I0, R0 = S_temp[-1], I_temp[-1], R_temp[-1]39 40 S_avg = np.mean(S_values, axis=0)41 I_avg = np.mean(I_values, axis=0)42 R_avg = np.mean(R_values, axis=0)43 44 return t_values, S_avg, I_avg, R_avg 45 46# 初期条件とパラメータの設定47S0 = 9948I0 = 149R0 = 050beta = 4/951gamma = 1/752t_max = 1353dt = 0.0154num_simulations = 1055 56# シミュレーションの繰り返し57t_values, S_avg, I_avg, R_avg = sir_simulation_stochastic(S0, I0, R0, beta, gamma, t_max, dt, num_simulations)58 59# グラフの描画60plt.figure(figsize=(10, 6))61plt.plot(t_values, S_avg, label='S(t)')62plt.plot(t_values, I_avg, label='I(t)')63plt.plot(t_values, R_avg, label='R(t)')64 65plt.xlabel('Days')66plt.ylabel('Population')67plt.title('Stochastic SIR Model Simulation')68plt.legend()69plt.grid(True)70plt.show()71