SIS model with no demographics¶
Moving on from the SI model to something slighlty more complex, we will add a state transition from infected back to susceptible, with no period of immunity - this is the SIS model. It turns out that the equations governing this will be analogous to the SI model with demography - rather than births & mortality providing routes to add susceptibles and subtract infectives, this transition from infective -> susceptible will provide essentially the same mechanism.
$$ \dot{S} = -\frac{\beta*S*I}{N} + \gamma I $$
$$ \dot{I} = \frac{\beta*S*I}{N} - \gamma I $$
As before, subbing $S = N-I$ into the second equation gives us
$$ \dot{I} = \beta I ( 1-\frac{\gamma}{\beta}-\frac{I}{N})$$
And comparing against the SI model with births, it is clear this equation is of the same form, with solution $$ \frac{Nx}{1+(\frac{Nx}{I_0}-1)e^{-\beta x t}} $$
$$ x = (1-\frac{\gamma}{\beta}) $$
This notebook tests the implementation and behavior of the model as follows:
Construct the model¶
In the first few cells, we do all the necessary imports. Then we construct a single-patch LASER model with three components: Susceptibility, Transmission, and Infection_SIS - this component will require a new agent property itimer, and upon expiration of itimer agents will return to the susceptible state. Finally, we initialize with a single infection and run. The Susceptibility and Transmission components are previously described.
Sanity check¶
The first test ensures certain basic constraints are being obeyed by the model. We confirm that at each timestep, $S_t=N_t-I_t$.
Scientific test¶
Finally, we come to the scientific test. As before, we first test on a single instance of the model and show that the expected output is recovered. Then, we select a few values of $\beta$ and $\gamma$, run the model, fit the outputs to the logistic equation, and compare the fitted value of $\beta$ and $\gamma$ to the known values; all of the considerations noted in the SI with births model, in terms of how to approach this fit, are echoed again here. Of particular concern is the approximation of an exponential transition from infected back to susceptible - as we are doing a first-order finite timestep integration here, that approximation will probably produce an error linear in $\gamma \Delta t$ between the analytic result and the modeled result. In fact, in a lot of real disease models, we have compartment dwell times in the exposed and infective states that are only a handful of $\Delta t$ long, but when doing real epi modeling and calibrating model parameters to uncertain data, this is generally not likely to be a dominant source of bias, uncertainty, etc. But when comparing specifically against an analytic result, it can become significant.
In [1]:
Copied!
import matplotlib.pyplot as plt
import numba as nb
import numpy as np
import pandas as pd
from laser.core.propertyset import PropertySet
from scipy.optimize import curve_fit
from scipy.optimize import minimize
import laser.core
import laser.generic
print(f"{np.__version__=}")
print(f"{laser.core.__version__=}")
print(f"{laser.generic.__version__=}")
import matplotlib.pyplot as plt
import numba as nb
import numpy as np
import pandas as pd
from laser.core.propertyset import PropertySet
from scipy.optimize import curve_fit
from scipy.optimize import minimize
import laser.core
import laser.generic
print(f"{np.__version__=}")
print(f"{laser.core.__version__=}")
print(f"{laser.generic.__version__=}")
np.__version__='2.3.5' laser.core.__version__='0.8.0' laser.generic.__version__='0.0.0'
In [2]:
Copied!
import laser.core.distributions as dists
import laser.generic.SIS as SIS
from laser.generic import Model
from laser.core.utils import grid
import laser.core.distributions as dists
import laser.generic.SIS as SIS
from laser.generic import Model
from laser.core.utils import grid
To make sure we don't accumulate lots of finite time-step error, make inf mean quite long in units of timestep
In [3]:
Copied!
pop = 3e5
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 1
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
parameters = PropertySet({"prng_seed": 42, "nticks": 3000, "beta": 0.1, "inf_mean": 32})
infdurdist = dists.exponential(scale=parameters.inf_mean)
pop = 3e5
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 1
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
parameters = PropertySet({"prng_seed": 42, "nticks": 3000, "beta": 0.1, "inf_mean": 32})
infdurdist = dists.exponential(scale=parameters.inf_mean)
In [4]:
Copied!
# Run until we get an outbreak
outbreak = False
while not outbreak:
parameters.prng_seed += 1
model = Model(scenario, parameters)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run()
outbreak = np.any(model.nodes.I[200] > 0)
# Run until we get an outbreak
outbreak = False
while not outbreak:
parameters.prng_seed += 1
model = Model(scenario, parameters)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run()
outbreak = np.any(model.nodes.I[200] > 0)
300,000 agents in 1 node(s): 100%|██████████| 3000/3000 [00:01<00:00, 1976.46it/s] 300,000 agents in 1 node(s): 100%|██████████| 3000/3000 [00:01<00:00, 1950.55it/s]
Sanity checks¶
Check that the relationships between susceptible, infected, and total population hold.
In [5]:
Copied!
fig, axs = plt.subplots(1, 2, figsize=(18, 5))
# Panel 2: Susceptible over time
axs[0].plot(model.nodes.S[1:], lw=4)
axs[0].plot(model.nodes.S[:-1] + model.nodes.newly_recovered[:-1] - model.nodes.newly_infected[:-1])
axs[0].set_yscale("log")
axs[0].set_title("Susceptible over time")
axs[0].set_xlabel("Time")
axs[0].set_ylabel("Susceptible")
# Panel 3: Population minus cumulative infections (incidence)
axs[1].plot(model.nodes.I[1:], lw=4)
axs[1].plot(model.nodes.I[:-1] - model.nodes.newly_recovered[:-1] + model.nodes.newly_infected[:-1])
axs[1].set_yscale("log")
axs[1].set_title("Population minus cumulative infections")
axs[1].set_xlabel("Time")
axs[1].set_ylabel("Incidence")
plt.tight_layout()
print("S[t] = S[t-1] + recovered[t-1] - incidence[t-1]: ", np.allclose(model.nodes.S[1:], model.nodes.S[:-1] + model.nodes.newly_recovered[:-1] - model.nodes.newly_infected[:-1]))
print("I[t] = I[t-1] - recovered[t-1] + incidence[t-1]: ", np.allclose(model.nodes.I[1:], model.nodes.I[:-1] - model.nodes.newly_recovered[:-1] + model.nodes.newly_infected[:-1]))
fig, axs = plt.subplots(1, 2, figsize=(18, 5))
# Panel 2: Susceptible over time
axs[0].plot(model.nodes.S[1:], lw=4)
axs[0].plot(model.nodes.S[:-1] + model.nodes.newly_recovered[:-1] - model.nodes.newly_infected[:-1])
axs[0].set_yscale("log")
axs[0].set_title("Susceptible over time")
axs[0].set_xlabel("Time")
axs[0].set_ylabel("Susceptible")
# Panel 3: Population minus cumulative infections (incidence)
axs[1].plot(model.nodes.I[1:], lw=4)
axs[1].plot(model.nodes.I[:-1] - model.nodes.newly_recovered[:-1] + model.nodes.newly_infected[:-1])
axs[1].set_yscale("log")
axs[1].set_title("Population minus cumulative infections")
axs[1].set_xlabel("Time")
axs[1].set_ylabel("Incidence")
plt.tight_layout()
print("S[t] = S[t-1] + recovered[t-1] - incidence[t-1]: ", np.allclose(model.nodes.S[1:], model.nodes.S[:-1] + model.nodes.newly_recovered[:-1] - model.nodes.newly_infected[:-1]))
print("I[t] = I[t-1] - recovered[t-1] + incidence[t-1]: ", np.allclose(model.nodes.I[1:], model.nodes.I[:-1] - model.nodes.newly_recovered[:-1] + model.nodes.newly_infected[:-1]))
S[t] = S[t-1] + recovered[t-1] - incidence[t-1]: True I[t] = I[t-1] - recovered[t-1] + incidence[t-1]: True
Single-simulation check¶
As before, starting with a single infection induces some stochasticity in terms of when the outbreak really starts to take off, and so we fit the expected behavior with a free offset parameter below.
In [6]:
Copied!
def SIS_logistic(t, beta, popsize, gamma, t0):
x = 1 - gamma / beta
return popsize * x / (1 + (popsize * x - 1) * np.exp(-beta * x * (t - t0)))
t = np.arange(model.params.nticks+1)
def objective(t0):
return np.sum(
(1 - SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, t0) / np.squeeze(model.nodes.I)) ** 2
)
result = minimize(objective, x0=10)
t0_opt = result.x[0]
plt.plot(model.nodes.I, lw=4)
plt.plot(SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, 0), lw=3)
plt.plot(SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, t0_opt), "r:", lw=3)
plt.yscale("log")
plt.legend(["Model output", "Logistic growth with known inputs, t0=0", f"Logistic growth with known inputs, best-fit t0 = {t0_opt:.1f}"])
def SIS_logistic(t, beta, popsize, gamma, t0):
x = 1 - gamma / beta
return popsize * x / (1 + (popsize * x - 1) * np.exp(-beta * x * (t - t0)))
t = np.arange(model.params.nticks+1)
def objective(t0):
return np.sum(
(1 - SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, t0) / np.squeeze(model.nodes.I)) ** 2
)
result = minimize(objective, x0=10)
t0_opt = result.x[0]
plt.plot(model.nodes.I, lw=4)
plt.plot(SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, 0), lw=3)
plt.plot(SIS_logistic(t, model.params.beta, pop, 1 / model.params.inf_mean, t0_opt), "r:", lw=3)
plt.yscale("log")
plt.legend(["Model output", "Logistic growth with known inputs, t0=0", f"Logistic growth with known inputs, best-fit t0 = {t0_opt:.1f}"])
Out[6]:
<matplotlib.legend.Legend at 0x124a66360>
Scientific testing¶
Finally, we run the model for a range of $\beta$ & $\gamma$ parameters, we freely fit the model output to the logistic equation, and we compare the known input parameters against the parameters fitted from output.
We will use only relatively large values of $\gamma$ for this procedure. The reason why will become clear in a second test, where we demonstrate that there is an error between the expected final size and the modeled final size, and that this error shrinks with $\gamma \Delta t$, as we would expect a first-order approximation error to.
To make this a pass-fail test, we will raise a flag if the fitted parameters are more than 5% different than the known ones.
In [7]:
Copied!
NTESTS = 10
nticks = 730
t = np.arange(nticks)
betarange = [0.02, 0.1]
gammarange = [1 / 300, 1 / 100]
seeds = list(range(NTESTS))
pop = 1e5
betas = np.random.uniform(betarange[0], betarange[1], NTESTS)
gammas = np.random.uniform(gammarange[0], gammarange[1], NTESTS)
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 3
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
output = []
for i, (seed, beta, gamma) in enumerate(zip(seeds, betas, gammas)):
parameters = PropertySet({"prng_seed": seed, "nticks": nticks, "verbose": True, "beta": beta, "inf_mean": 1 / gamma})
model = Model(scenario, parameters)
infdurdist = dists.exponential(scale=parameters.inf_mean)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run(label=f"SIS {i+1:2} of {NTESTS}, seed={seed}, beta={beta:.3f}, gamma={gamma:.5f}")
cases = model.nodes.I[1:,0]
popt, pcov = curve_fit(
SIS_logistic,
t,
cases,
p0=[np.mean(betarange), pop, np.mean(gammarange), 1],
bounds=([betarange[0] / 2, pop - 1, gammarange[0] / 2, -300], [betarange[1] * 2, pop + 1, gammarange[1] * 2, 300]),
)
output.append(
{
"seed": seed,
"beta": beta,
"gamma": gamma,
"cases": [np.array(cases)],
"fitted_beta": popt[0],
"fitted_gamma": popt[2],
"fitted_t0": popt[3],
"recovered": [np.array(model.nodes.newly_recovered[1:,0])]
}
)
output = pd.DataFrame(output)
NTESTS = 10
nticks = 730
t = np.arange(nticks)
betarange = [0.02, 0.1]
gammarange = [1 / 300, 1 / 100]
seeds = list(range(NTESTS))
pop = 1e5
betas = np.random.uniform(betarange[0], betarange[1], NTESTS)
gammas = np.random.uniform(gammarange[0], gammarange[1], NTESTS)
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 3
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
output = []
for i, (seed, beta, gamma) in enumerate(zip(seeds, betas, gammas)):
parameters = PropertySet({"prng_seed": seed, "nticks": nticks, "verbose": True, "beta": beta, "inf_mean": 1 / gamma})
model = Model(scenario, parameters)
infdurdist = dists.exponential(scale=parameters.inf_mean)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run(label=f"SIS {i+1:2} of {NTESTS}, seed={seed}, beta={beta:.3f}, gamma={gamma:.5f}")
cases = model.nodes.I[1:,0]
popt, pcov = curve_fit(
SIS_logistic,
t,
cases,
p0=[np.mean(betarange), pop, np.mean(gammarange), 1],
bounds=([betarange[0] / 2, pop - 1, gammarange[0] / 2, -300], [betarange[1] * 2, pop + 1, gammarange[1] * 2, 300]),
)
output.append(
{
"seed": seed,
"beta": beta,
"gamma": gamma,
"cases": [np.array(cases)],
"fitted_beta": popt[0],
"fitted_gamma": popt[2],
"fitted_t0": popt[3],
"recovered": [np.array(model.nodes.newly_recovered[1:,0])]
}
)
output = pd.DataFrame(output)
SIS 1 of 10, seed=0, beta=0.030, gamma=0.00478: 100%|██████████| 730/730 [00:00<00:00, 1600.42it/s] SIS 2 of 10, seed=1, beta=0.066, gamma=0.00870: 100%|██████████| 730/730 [00:00<00:00, 1832.14it/s] SIS 3 of 10, seed=2, beta=0.089, gamma=0.00488: 100%|██████████| 730/730 [00:00<00:00, 1870.13it/s] SIS 4 of 10, seed=3, beta=0.068, gamma=0.00471: 100%|██████████| 730/730 [00:00<00:00, 1859.55it/s] SIS 5 of 10, seed=4, beta=0.064, gamma=0.00768: 100%|██████████| 730/730 [00:00<00:00, 1555.22it/s] SIS 6 of 10, seed=5, beta=0.050, gamma=0.00816: 100%|██████████| 730/730 [00:00<00:00, 1819.16it/s] SIS 7 of 10, seed=6, beta=0.100, gamma=0.00864: 100%|██████████| 730/730 [00:00<00:00, 1840.22it/s] SIS 8 of 10, seed=7, beta=0.081, gamma=0.00601: 100%|██████████| 730/730 [00:00<00:00, 1833.60it/s] SIS 9 of 10, seed=8, beta=0.086, gamma=0.00476: 100%|██████████| 730/730 [00:00<00:00, 1865.35it/s] SIS 10 of 10, seed=9, beta=0.080, gamma=0.00903: 100%|██████████| 730/730 [00:00<00:00, 1862.75it/s]
In [8]:
Copied!
plt.figure()
plt.plot(output["beta"], output["fitted_beta"], "o")
plt.xlim(betarange[0], betarange[1])
plt.ylim(betarange[0], betarange[1])
plt.figure()
plt.plot(output["beta"], 1 - output["beta"] / output["fitted_beta"], "o")
plt.xlim(betarange[0], betarange[1])
plt.ylim(-0.25, 0.25)
plt.figure()
plt.plot(output["gamma"], output["fitted_gamma"], "o")
plt.xlim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.ylim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.figure()
plt.plot(output["gamma"], 1 - output["gamma"] / output["fitted_gamma"], "o")
plt.xlim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.ylim(-0.25, 0.25)
plt.figure()
plt.plot(output["beta"], output["fitted_beta"], "o")
plt.xlim(betarange[0], betarange[1])
plt.ylim(betarange[0], betarange[1])
plt.figure()
plt.plot(output["beta"], 1 - output["beta"] / output["fitted_beta"], "o")
plt.xlim(betarange[0], betarange[1])
plt.ylim(-0.25, 0.25)
plt.figure()
plt.plot(output["gamma"], output["fitted_gamma"], "o")
plt.xlim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.ylim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.figure()
plt.plot(output["gamma"], 1 - output["gamma"] / output["fitted_gamma"], "o")
plt.xlim(gammarange[0]*0.9, gammarange[1]*1.1)
plt.ylim(-0.25, 0.25)
Out[8]:
(-0.25, 0.25)
In [9]:
Copied!
print(
"All fitted beta are within 10% of known beta: " + str(np.all(np.abs((output["beta"] - output["fitted_beta"]) / output["beta"]) < 0.10))
)
print(
"All fitted gamma are within 20% of known gamma: "
+ str(np.all(np.abs((output["gamma"] - output["fitted_gamma"]) / output["gamma"]) < 0.20))
)
print(
"All fitted beta are within 10% of known beta: " + str(np.all(np.abs((output["beta"] - output["fitted_beta"]) / output["beta"]) < 0.10))
)
print(
"All fitted gamma are within 20% of known gamma: "
+ str(np.all(np.abs((output["gamma"] - output["fitted_gamma"]) / output["gamma"]) < 0.20))
)
All fitted beta are within 10% of known beta: True All fitted gamma are within 20% of known gamma: True
In [10]:
Copied!
output
output
Out[10]:
| seed | beta | gamma | cases | fitted_beta | fitted_gamma | fitted_t0 | recovered | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0.029660 | 0.004783 | [[3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5,... | 0.029353 | 0.004690 | -28.519006 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 1 | 1 | 0.066200 | 0.008704 | [[4, 4, 4, 4, 5, 6, 7, 8, 8, 9, 10, 12, 12, 12... | 0.064577 | 0.008226 | -30.855025 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 2 | 2 | 0.088750 | 0.004881 | [[4, 4, 4, 6, 7, 9, 10, 10, 10, 12, 16, 19, 19... | 0.086452 | 0.004554 | -24.615804 | [[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 3 | 3 | 0.067663 | 0.004713 | [[4, 4, 4, 4, 5, 6, 7, 8, 8, 9, 10, 12, 12, 12... | 0.065749 | 0.004419 | -27.636157 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 4 | 4 | 0.063574 | 0.007682 | [[4, 4, 4, 4, 5, 6, 7, 8, 9, 10, 12, 14, 14, 1... | 0.061765 | 0.007219 | -38.782883 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 5 | 5 | 0.049859 | 0.008158 | [[3, 3, 3, 3, 3, 3, 4, 4, 3, 3, 4, 6, 6, 6, 6,... | 0.048996 | 0.007822 | -38.706199 | [[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,... |
| 6 | 6 | 0.099783 | 0.008644 | [[4, 4, 4, 6, 7, 7, 8, 8, 8, 10, 13, 15, 15, 1... | 0.097077 | 0.008029 | -18.788862 | [[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,... |
| 7 | 7 | 0.080933 | 0.006014 | [[3, 3, 3, 3, 4, 4, 5, 6, 7, 8, 10, 11, 10, 10... | 0.079083 | 0.005661 | -18.898113 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,... |
| 8 | 8 | 0.086307 | 0.004760 | [[4, 3, 3, 3, 4, 4, 5, 6, 7, 8, 9, 9, 10, 10, ... | 0.084231 | 0.004460 | -13.418348 | [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
| 9 | 9 | 0.080437 | 0.009030 | [[3, 3, 3, 3, 4, 4, 5, 6, 7, 8, 9, 9, 10, 10, ... | 0.078341 | 0.008463 | -13.544913 | [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,... |
Quick demonstration of first-order error accumulation¶
As noted before, for $\gamma \Delta t$ large, first-order integration like we are doing here can accumulate substantial error. Calculating exactly how error will accumulate in an integrator for a dynamic process like this is beyond the scope here, and probably depends on a lot of details. E.g., the ordering of steps - in a given step, does the transmission update from S->I state occur before or after the infection update that sends agents from I->S? Do we use midpoint methods, timer countdowns, or take advantage of the unique memorylessness of the exponential distribution to simply remove a random fraction each time? All that is beyond scope here, but just want to demonstrate that the error in the equilibrium value $I(t \rightarrow \infty)$ becomes large when the mean infectious period $\frac{1}{\gamma}$ is on the same order as $\Delta t$
In [11]:
Copied!
# %%capture
gammas = [1 / infmean for infmean in [1, 1.5, 2, 2.5, 3, 6, 12, 18, 30, 45, 60, 90, 120, 180, 240, 300]]
betas = [3 * gamma for gamma in gammas]
NTESTS = len(gammas)
nticks = 3000
seeds = list(range(NTESTS))
pop = 1e5
final_expected = np.array([])
final_observed = np.array([])
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 20
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
for i, (seed, beta, gamma) in enumerate(zip(seeds, betas, gammas)):
parameters = PropertySet({"prng_seed": seed, "nticks": nticks, "verbose": True, "beta": beta, "inf_mean": 1 / gamma})
model = Model(scenario, parameters)
infdurdist = dists.exponential(scale=parameters.inf_mean)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run(label=f"SIS {i+1:2} of {NTESTS}, seed={seed:2}, beta={beta:.3f}, gamma={gamma:.5f}")
final_observed = np.append(final_observed, model.nodes.I[-1,0])
final_expected = np.append(final_expected, pop * (1 - gamma / beta))
# %%capture
gammas = [1 / infmean for infmean in [1, 1.5, 2, 2.5, 3, 6, 12, 18, 30, 45, 60, 90, 120, 180, 240, 300]]
betas = [3 * gamma for gamma in gammas]
NTESTS = len(gammas)
nticks = 3000
seeds = list(range(NTESTS))
pop = 1e5
final_expected = np.array([])
final_observed = np.array([])
scenario = grid(M=1, N=1, population_fn=lambda x,y: pop, origin_x=0, origin_y=0)
initial_infected = 20
scenario["S"] = scenario.population - initial_infected
scenario["I"] = initial_infected
for i, (seed, beta, gamma) in enumerate(zip(seeds, betas, gammas)):
parameters = PropertySet({"prng_seed": seed, "nticks": nticks, "verbose": True, "beta": beta, "inf_mean": 1 / gamma})
model = Model(scenario, parameters)
infdurdist = dists.exponential(scale=parameters.inf_mean)
model.components = [SIS.Susceptible(model), SIS.Infectious(model, infdurdist), SIS.Transmission(model, infdurdist)]
model.run(label=f"SIS {i+1:2} of {NTESTS}, seed={seed:2}, beta={beta:.3f}, gamma={gamma:.5f}")
final_observed = np.append(final_observed, model.nodes.I[-1,0])
final_expected = np.append(final_expected, pop * (1 - gamma / beta))
SIS 1 of 16, seed= 0, beta=3.000, gamma=1.00000: 100%|██████████| 3000/3000 [00:01<00:00, 1816.37it/s] SIS 2 of 16, seed= 1, beta=2.000, gamma=0.66667: 100%|██████████| 3000/3000 [00:01<00:00, 1822.22it/s] SIS 3 of 16, seed= 2, beta=1.500, gamma=0.50000: 100%|██████████| 3000/3000 [00:01<00:00, 1791.96it/s] SIS 4 of 16, seed= 3, beta=1.200, gamma=0.40000: 100%|██████████| 3000/3000 [00:01<00:00, 1783.50it/s] SIS 5 of 16, seed= 4, beta=1.000, gamma=0.33333: 100%|██████████| 3000/3000 [00:01<00:00, 1914.60it/s] SIS 6 of 16, seed= 5, beta=0.500, gamma=0.16667: 100%|██████████| 3000/3000 [00:01<00:00, 2145.71it/s] SIS 7 of 16, seed= 6, beta=0.250, gamma=0.08333: 100%|██████████| 3000/3000 [00:01<00:00, 2372.33it/s] SIS 8 of 16, seed= 7, beta=0.167, gamma=0.05556: 100%|██████████| 3000/3000 [00:01<00:00, 2414.61it/s] SIS 9 of 16, seed= 8, beta=0.100, gamma=0.03333: 100%|██████████| 3000/3000 [00:01<00:00, 2282.53it/s] SIS 10 of 16, seed= 9, beta=0.067, gamma=0.02222: 100%|██████████| 3000/3000 [00:01<00:00, 2495.26it/s] SIS 11 of 16, seed=10, beta=0.050, gamma=0.01667: 100%|██████████| 3000/3000 [00:01<00:00, 2527.17it/s] SIS 12 of 16, seed=11, beta=0.033, gamma=0.01111: 100%|██████████| 3000/3000 [00:01<00:00, 2563.65it/s] SIS 13 of 16, seed=12, beta=0.025, gamma=0.00833: 100%|██████████| 3000/3000 [00:01<00:00, 2569.93it/s] SIS 14 of 16, seed=13, beta=0.017, gamma=0.00556: 100%|██████████| 3000/3000 [00:01<00:00, 2420.63it/s] SIS 15 of 16, seed=14, beta=0.013, gamma=0.00417: 100%|██████████| 3000/3000 [00:01<00:00, 2557.22it/s] SIS 16 of 16, seed=15, beta=0.010, gamma=0.00333: 100%|██████████| 3000/3000 [00:01<00:00, 2612.99it/s]
In [ ]:
Copied!
plt.plot(gammas, np.abs(1 - final_observed / final_expected), "o")
plt.xlabel(r"$\gamma$")
plt.ylabel("$| 1 - \\frac{I(\\infty)_{obs}}{I(\\infty)_{exp}} |$")
plt.title(r"Error in equilibrium infected fraction increases roughly linearly in $\gamma \Delta t$")
plt.show()
plt.plot(gammas, np.abs(1 - final_observed / final_expected), "o")
plt.xlabel(r"$\gamma$")
plt.ylabel("$| 1 - \\frac{I(\\infty)_{obs}}{I(\\infty)_{exp}} |$")
plt.title(r"Error in equilibrium infected fraction increases roughly linearly in $\gamma \Delta t$")
plt.show()
Out[ ]:
Text(0.5, 1.0, 'Error in equilibrium infected fraction increases roughly linearly in $\\gamma \\Delta t$')
In [13]:
Copied!
# Validate for each row that mean(recovered / cases) is close to gamma
for idx, row in output.iterrows():
cases_arr = np.array(row["cases"][0])
recovered_arr = np.array(row["recovered"][0])
# Avoid division by zero
valid = cases_arr > 0
avg_ratio = np.mean(recovered_arr[valid] / cases_arr[valid])
print(f"Row {idx}: gamma={row['gamma']:.6f}, avg(recovered/cases)={avg_ratio:.6f}, diff={abs(avg_ratio - row['gamma']):.6f}")
# Validate for each row that mean(recovered / cases) is close to gamma
for idx, row in output.iterrows():
cases_arr = np.array(row["cases"][0])
recovered_arr = np.array(row["recovered"][0])
# Avoid division by zero
valid = cases_arr > 0
avg_ratio = np.mean(recovered_arr[valid] / cases_arr[valid])
print(f"Row {idx}: gamma={row['gamma']:.6f}, avg(recovered/cases)={avg_ratio:.6f}, diff={abs(avg_ratio - row['gamma']):.6f}")
Row 0: gamma=0.004783, avg(recovered/cases)=0.004250, diff=0.000533 Row 1: gamma=0.008704, avg(recovered/cases)=0.008376, diff=0.000329 Row 2: gamma=0.004881, avg(recovered/cases)=0.004815, diff=0.000065 Row 3: gamma=0.004713, avg(recovered/cases)=0.004469, diff=0.000245 Row 4: gamma=0.007682, avg(recovered/cases)=0.007357, diff=0.000325 Row 5: gamma=0.008158, avg(recovered/cases)=0.008077, diff=0.000081 Row 6: gamma=0.008644, avg(recovered/cases)=0.008481, diff=0.000163 Row 7: gamma=0.006014, avg(recovered/cases)=0.005878, diff=0.000136 Row 8: gamma=0.004760, avg(recovered/cases)=0.004952, diff=0.000192 Row 9: gamma=0.009030, avg(recovered/cases)=0.008840, diff=0.000190