<jupyter_start><jupyter_text>General<jupyter_code>import numpy as np
import random
import math
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from tqdm import tqdm


import os

try:
  os.mkdir("figures")
except OSError:
  pass

#plt.rcParams.update({'figure.dpi': 200})
#plt.rcParams.update({'font.size': 22})<jupyter_output><empty_output><jupyter_text>Deterministic solution<jupyter_code>def SIR_derivative(u, t, rate_infection, rate_recovery):
    S, I, R = u
    N = S + I + R
    dSdt = -rate_infection * S * I / N
    dIdt =  rate_infection * S * I / N - rate_recovery * I
    dRdt =  rate_recovery * I
    return dSdt, dIdt, dRdt

def SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times):
  u0 = S0, I0, R0
  
  # Integrate SIR equations over the time grid, t.
  sol = odeint(SIR_derivative, u0, times, args=(rate_infection, rate_recovery))
  S, I, R = sol.T
  return np.array(S), np.array(I), np.array(R)

def exponential_growth(growth_rate, I0, times):
  # Integrate exponential growth
  def derivative(I,t,growth_rate):
    return growth_rate*I
  
  sol = odeint(derivative, I0, times, args=(growth_rate,))
  I, = sol.T
  return I<jupyter_output><empty_output><jupyter_text>Kinetic Monte Carlo Simulation der SIR Dynamik<jupyter_code>def simulate(rate_infection, rate_recovery, S0, I0, R0, times):
  S = S0; I = I0; R = R0; 
  time = 0;

  # measurement variables
  trajectory_S = np.zeros(len(times))
  trajectory_I = np.zeros(len(times))
  trajectory_R = np.zeros(len(times))
  i_meas = 0
  if times[i_meas] < 0:
    print("error")
    return -1
  if times[i_meas]==0:
    trajectory_S[0] = S0
    trajectory_I[0] = I0
    trajectory_R[0] = R0
    i_meas = 1

  # kinetic Monte Carlo simulation loop
  while I > 0 and time < times[-1]:
    dt, event = new_event(rate_infection, rate_recovery, S, I, R)
    time += dt

    # measure before update of system, because time will always surpass the 
    # measurement time
    while time > times[i_meas]:
      trajectory_S[i_meas] = S
      trajectory_I[i_meas] = I
      trajectory_R[i_meas] = R
      i_meas += 1
      if i_meas > len(times)-1:
        break
    
    # update system
    if event == "recovery":
      I -= 1
      R += 1
    if event == "infection":
      S -= 1
      I += 1

  return trajectory_S, trajectory_I, trajectory_R
    
def new_event(rate_infection, rate_recovery, S, I, R):
  N = S + I + R

  # calculate rates
  total_rate_infection = rate_infection * I * S / N
  total_rate_recovery  = rate_recovery  * I
  total_rate = total_rate_infection + total_rate_recovery

  # draw random exponentially distributed inter event time
  dt = np.random.exponential(1/total_rate)

  # select event with weight determinde by rate 
  event = ""
  if np.random.rand() < total_rate_infection/total_rate:
    event = "infection"
  else:
    event = "recovery"

  return dt, event<jupyter_output><empty_output><jupyter_text>Example<jupyter_code>rate_recovery  = 1/7
rate_infection = 2/7
I0 = 10
N = int(1e5)
S0 = N-I0
R0 = 0
time_total = 21

np.random.seed(1000)

num_trajectories = 10000
times = np.arange(0, time_total, 1)
list_Is_10 = [np.zeros(num_trajectories) for t in range(len(times))]
for i in tqdm(range(num_trajectories)):
  S, I, R = simulate(rate_infection, rate_recovery, S0, I0, R0, times)
  for t in range(len(times)):
    list_Is_10[t][i] = I[t]

  if i < 100:
    plt.plot(times, I, color='gray', linewidth=0.5)

# deterministic solution
S_det, I_det, R_det = SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times)
plt.plot(times, I_det, color='black', label='deterministic', linewidth=2)

plt.xlabel('time')
plt.ylabel('number active cases')
plt.legend()
plt.savefig('figures/kineticMC_example_I0_10.png', dpi=200)
fig, axs = plt.subplots(len(times), 1, 
                        sharex=True, 
                        figsize=(6,6))

axs[-1].set_xlabel('number active cases')

for (i,t) in enumerate(range(len(times))):
  axs[i].hist(list_Is_10[t], bins=np.arange(0,200,1));
  axs[i].set_ylim(0,200)
  axs[i].spines["top"].set_visible(False)
  axs[i].spines["right"].set_visible(False)

plt.savefig('figures/kineticMC_example_I0_10_distributions.png', dpi=200)
rate_recovery  = 1/7
rate_infection = 2/7
I0 = 10
N = int(1e5)
S0 = N-I0
R0 = 0
time_total = 120

np.random.seed(1000)

num_trajectories = 1000
times = np.arange(0, time_total, 1)
mean_trajectory = np.zeros(time_total)
for i in tqdm(range(num_trajectories)):
  S, I, R = simulate(rate_infection, rate_recovery, S0, I0, R0, times)
  mean_trajectory += I

  if i%10==0:
    plt.plot(times, I, color='gray', linewidth=0.5)

mean_trajectory /= num_trajectories

# deterministic solution
S_det, I_det, R_det = SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times)
plt.plot(times, I_det, color='black', label='deterministic', linewidth=2)

# mean
plt.plot(times, mean_trajectory, color='blue', label='mean', linewidth=2)

plt.xlabel('time')
plt.ylabel('number active cases')
plt.legend()
plt.savefig('figures/kineticMC_example_I0_10_long.png', dpi=200)<jupyter_output>100%|██████████| 1000/1000 [21:12<00:00,  1.27s/it]<jupyter_text>Extreme example<jupyter_code>rate_recovery  = 1/7.
rate_infection = 1.1/7.
I0 = 1
N = int(1e5)
S0 = N-I0
R0 = 0
time_total = 21

np.random.seed(1000)

num_trajectories = 10000
times = np.arange(0, time_total, 1)
list_Is_1 = [np.zeros(num_trajectories) for t in range(len(times))]
for i in tqdm(range(num_trajectories)):
  S, I, R = simulate(rate_infection, rate_recovery, S0, I0, R0, times)
  for t in range(len(times)):
    list_Is_1[t][i] = I[t]

  if i < 100:
    plt.plot(times, I, color='gray', linewidth=0.5)

# deterministic solution
S_det, I_det, R_det = SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times)
plt.plot(times, I_det, color='black', label='deterministic', linewidth=2)

plt.xlabel('time')
plt.ylabel('number active cases')
plt.legend()
plt.savefig('figures/kineticMC_example_I0_1.png', dpi=200)
fig, axs = plt.subplots(len(times), 1, 
                        sharex=True, 
                        figsize=(6,6))

axs[-1].set_xlabel('number active cases')

for (i,t) in enumerate(range(len(times))):
  axs[i].hist(list_Is_1[t], bins=np.arange(0,15,1));
  axs[i].set_ylim(0,10000)
  axs[i].spines["top"].set_visible(False)
  axs[i].spines["right"].set_visible(False)

plt.savefig('figures/kineticMC_example_I0_1_distributions.png', dpi=200)<jupyter_output><empty_output><jupyter_text>Bayesian Inference SIRMetropolis parameter sampling of SIR model:P(parameter|daten) = P(daten|parameter) * P(parameter)Choices:* P(parameter):       rate_infection in [0,1]   rate_recovery in [0,1]      I0 in [1,100]* P(daten|parameter) = (I(data)-I(model)))^2 [least squares]<jupyter_code># data source
np.random.seed(2000)
times = np.arange(1,120,1)
S, data_I_long, R = simulate(2/7, 1/7, S0, 10, R0, times)

plt.plot(times, data_I_long, 'o', label='stochastic sample')

plt.ylim(0,17000)
plt.xlabel('time')
plt.ylabel('number active cases')
plt.legend()
#plt.show()
plt.savefig('figures/inference_example_SIR_data.png', dpi=200)


# likelihood is product of Gaussians with variance sqrt(data)
def log_likelihood(data_I, model_I):
  return -np.sum(np.power(model_I-data_I,2)/2/data_I)

def parameter_proposal(p, dp):
  return p + (-1 + 2 * np.random.rand()) *dp

# known parameters
S0 = int(1e5)
R0 = 0
times = np.arange(1,120,1)

# initial parameter
I0 = 1
rate_infection = 0.1
rate_recovery  = 0.1
S,model_I,R = SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times)
logL = log_likelihood(data_I_long, model_I)

SIR_samples_I0 = [I0]
SIR_samples_rate_infection = [rate_infection]
SIR_samples_rate_recovery = [rate_recovery]
SIR_samples_logL = [logL]

SIR_best_fit = (logL, I0, rate_infection, rate_recovery)

np.random.seed(1000)
num_updates = int(1e5)
dI0 = 0.1
drate = 0.01
for i in tqdm(range(num_updates)):
  #propose new parameters
  new_I0 = I0
  new_rate_infection = rate_infection
  new_rate_recovery = rate_recovery
  p_select = np.random.rand()
  if (p_select < 1/3.):
    new_I0 = parameter_proposal(I0, dI0)
  elif (p_select < 2/3.):
    new_rate_infection = parameter_proposal(rate_infection, drate)
  else:
    new_rate_recovery  = parameter_proposal(rate_recovery,  drate)

  # accept or reject new parameters
  if (1 < new_I0 < 100) \
   and (0 < new_rate_recovery < 1) \
   and (0 < new_rate_infection < 1):
    S,model_I,R = SIR_deterministic(new_rate_infection, new_rate_recovery, \
                                    S0, new_I0, R0, times)
    new_logL = log_likelihood(data_I_long, model_I)
    #accept with weight P(new)/P(old)
    if np.random.rand() < np.exp(new_logL - logL): 
      I0 = new_I0
      rate_infection = new_rate_infection
      rate_recovery = new_rate_recovery
      logL = new_logL

  # store best fit
  if logL > SIR_best_fit[0]:
    SIR_best_fit = (logL, I0, rate_infection, rate_recovery)

  SIR_samples_I0.append(I0)
  SIR_samples_rate_infection.append(rate_infection)
  SIR_samples_rate_recovery.append(rate_recovery)
  SIR_samples_logL.append(logL)
fig, axs = plt.subplots(4, 2, sharey=False, tight_layout=True, figsize=(6,6))

num_equilibration = int(5e4)

axs[0,0].plot(SIR_samples_rate_recovery)
axs[1,0].plot(SIR_samples_rate_infection)
axs[2,0].plot(SIR_samples_I0)
axs[3,0].plot(SIR_samples_logL)
axs[0,1].hist(SIR_samples_rate_recovery[num_equilibration:], 100, orientation="horizontal");
axs[1,1].hist(SIR_samples_rate_infection[num_equilibration:], 100, orientation="horizontal");
axs[2,1].hist(SIR_samples_I0[num_equilibration:], 100, orientation="horizontal");
axs[3,1].hist(SIR_samples_logL[num_equilibration:], 100, orientation="horizontal");

axs[3,0].set_xlabel("Monte Carlo time")
axs[3,1].set_xlabel("distributions")

axs[0,0].set_ylabel("rate recovery")
axs[1,0].set_ylabel("rate infection")
axs[2,0].set_ylabel("I0")
axs[3,0].set_ylabel("log(likelihood)")

plt.savefig('figures/inference_example_SIR_MC.png', dpi=200)
S0 = int(1e5)
R0 = 0
I0 = 10
times = np.arange(1,120,1)
S, model_truth, R = SIR_deterministic(2/7, 1/7, S0, I0, R0, times)
print("truth:",log_likelihood(data_I_long, model_truth), 10, 2/7, 1/7)

print("best fit:",SIR_best_fit)
logL, I0, rate_infection, rate_recovery = SIR_best_fit
S, model_I, R = SIR_deterministic(rate_infection, rate_recovery, S0, I0, R0, times)

samples_rate_infection = np.array(SIR_samples_rate_infection)
samples_rate_recovery = np.array(SIR_samples_rate_recovery)
samples_logL = np.array(SIR_samples_logL)

print("expectation:", 
      np.sum(samples_rate_infection[num_equilibration:]*np.exp(samples_logL[num_equilibration:]))/np.sum(np.exp(samples_logL[num_equilibration:])),
      np.sum(samples_rate_recovery[num_equilibration:]*np.exp(samples_logL[num_equilibration:])/np.sum(np.exp(samples_logL[num_equilibration:]))) )


plt.plot(times, data_I_long, 'o')
plt.plot(times, model_truth, label='truth', linewidth=4)
plt.plot(times, model_I, label = 'best fit', linewidth=4)

plt.ylim(0,17000)
plt.legend()
plt.xlabel('time')
plt.ylabel('number active cases')
plt.savefig('figures//inference_example_SIR_bestfit.png', dpi=200)

#credibility interval
indices_sorted = samples_logL.argsort()<jupyter_output>truth: -42931.56105855961 10 0.2857142857142857 0.14285714285714285
best fit: (-74.36513594743387, 4.801935021066452, 0.2836168932292803, 0.14236636995404373)
expectation: 0.2834310748078799 0.14228767142942236<jupyter_text>Bayesian Inference exponential growthMetropolis parameter sampling of exponential growth:P(parameter|daten) = P(daten|parameter) * P(parameter)Choices:* P(parameter):       growth_rate uniform in [-3,3]      I0 in [1,100]* P(daten|parameter) = (I(data)/I(model) - 1))^2 [least squares]<jupyter_code># random sample (I0=10, rate_infection=3/7, rate_recovery=1/7)
times = np.arange(0,21,1)
S_sample = np.array([99990,99985,99978,99974,99968,99960,99955,99946,99927,99911,99895,99876,99856,99824,99804,99764,99722,99682,99637,99589,99527])
I_sample = np.array([10,13,19,18,19,25,27,35,50,61,67,76,81,97,106,135,162,176,192,216,249])
S_det, I_det, R_det = SIR_deterministic(2/7, 1/7, int(1e5)-10, 10, 0, times)
I_noise = I_det + np.random.normal(0,np.sqrt(I_det),len(times))

plt.plot(times, I_sample, 'o', label='stochastic sample')
plt.plot(times, I_noise,  'o', label='deterministic plus Gaussian noise')

plt.xlabel('time')
plt.ylabel('number active cases')
plt.legend()
plt.savefig('figures/inference_example_exponential_data.png', dpi=200)
def log_likelihood(data_I, model_I):
  return -np.sum(np.power(model_I-data_I,2)/2/data_I)

def parameter_proposal(p, dp):
  return p + (-1 + 2 * np.random.rand())*dp

data_I = I_sample

I0 = 1
growth_rate = 0
model_I = exponential_growth(growth_rate, I0, times)
logL = log_likelihood(data_I, model_I)

samples_I0 = [I0]
samples_growth_rate = [growth_rate]
samples_logL = [logL]

best_fit = (logL, I0, growth_rate)

np.random.seed(1000)
num_updates = int(1e5)
dI0 = 0.1
drate = 0.01
for i in tqdm(range(num_updates)):
  #propose new parameters
  p_select = np.random.rand()
  new_I0 = I0
  new_growth_rate = growth_rate
  if p_select < 0.5:
    new_I0 = parameter_proposal(I0, dI0)
  else:
    new_growth_rate = parameter_proposal(growth_rate, drate)

  if (1 <new_I0 < 100) and (-3 < new_growth_rate < 3):
    model_I = exponential_growth(new_growth_rate, new_I0, times)
    new_logL = log_likelihood(data_I, model_I)
    if np.random.rand() < np.exp(new_logL-logL): # accept
      I0 = new_I0
      growth_rate = new_growth_rate
      logL = new_logL

  if logL > best_fit[0]:
    best_fit = (logL, I0, growth_rate)

  samples_I0.append(I0)
  samples_growth_rate.append(growth_rate)
  samples_logL.append(logL)
fig, axs = plt.subplots(3, 2, sharey=False, tight_layout=True, figsize=(6,6))

num_equilibration = int(10000)
axs[0,0].plot(samples_growth_rate)
axs[1,0].plot(samples_I0)
axs[2,0].plot(samples_logL)
axs[0,1].hist(samples_growth_rate[num_equilibration:], 100, orientation="horizontal");
axs[1,1].hist(samples_I0[num_equilibration:], 100, orientation="horizontal");
axs[2,1].hist(samples_logL[num_equilibration:], 100, orientation="horizontal");

axs[2,0].set_xlabel("Monte Carlo time")
axs[2,1].set_xlabel("distributions")

axs[0,0].set_ylabel("rate growth")
axs[1,0].set_ylabel("I0")
axs[2,0].set_ylabel("log(likelihood)")

plt.savefig('figures/inference_example_exponential_MC.png', dpi=200)
times=np.arange(0,21,1)
model_truth = exponential_growth(1/7, 10, times)
print("truth:",log_likelihood(data_I, model_truth), 10, 1/7)

print("best fit:",best_fit)
logL, I0, growth_rate = best_fit
model_I = exponential_growth(growth_rate, I0, times)

samples_growth_rate = np.array(samples_growth_rate)
samples_logL = np.array(samples_logL)

print("expectation:", 
      np.sum(samples_growth_rate[num_equilibration:]*np.exp(samples_logL[num_equilibration:]))/np.sum(np.exp(samples_logL[num_equilibration:]))
      )

plt.plot(times, data_I, 'o')
plt.plot(times, model_truth, label='truth', linewidth=4)
plt.plot(times, model_I, label = 'best fit', linewidth=4)

plt.xlabel('time')
plt.ylabel('number active infected')
plt.legend()
plt.savefig('figures/inference_example_exponential_bestfit.png', dpi=200)<jupyter_output><empty_output><jupyter_text>Export<jupyter_code>! zip figures.zip figures/*<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1412126365360717, 'fraction_numerical': 0.058211786251553405, 'mean_word_length': 5.302555647155812, 'pattern_counts': {'":': 2, '<': 52, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '197981', 'license': 'MIT', 'n_tokens_mistral': 6839, 'n_tokens_neox': 6169, 'n_words': 1125, 'path': 'lecture_MonteCarlo_SIR.ipynb', 'repo_name': 'zierenberg/lecture_MC_disease', 'size': '15289'}