0

I'm trying to make an SIR model for disease in python.

import random
import numpy as np
import math
import matplotlib.pyplot as plt
seed=random.seed()
nums=1000
fract=0.01
recovr=0.2
num_infect=1
numd=40
reps=40
data=np.full(numd,1)
sdata=np.full(numd,1)
remov=0
remove=0
for k in range(reps):
    remove=0                    #remove per day
    P=np.full(nums,1)
    n=0
    while n<math.trunc(fract*nums):
        if P[random.randint(0,nums-1)]!=0:
            n+=1
    n=0
    while n<1:
        if P[random.randint(0,nums-1)]==1:
            P[random.randint(0,nums-1)]=2
            n=1
    for j in range(numd):
        for i in range(nums):
            if P[i]==2:
                for ii in range(num_infect):
                    if P[random.randint(0,nums-1)]==1:
                        P[random.randint(0,nums-1)]=2
                        remov+=1
            if random.random()<recovr:
                if P[random.randint(0,nums-1)]==2:
                    P[random.randint(0,nums-1)]=0
        c,s=0,0
        for i in range(nums):
            if P[i]==2:
                c+=1
            if P[i]==1:
                s+=1
        data[j]=(data[j]+c)
        sdata[j]=(sdata[j]+s)

for i in range(numd):              #average after repeating reps times
    data[i]=data[i]//numd
    sdata[i]=sdata[i]//numd
num_infected=remov/reps
print(num_infected)

So the maximum number of people that could be infected should be 1000 but for some reason num_infected comes out as number around 1700. Like it's very confusing because it should add 1 onto the counter (remov) for every person that gets the disease but for some reason it goes past 1000 even though I set the maximum number of people to 1000 and no way to be reinfected after becoming immune. It does remov/reps because the code repeats 40 times and that finds the average.

Jet
  • 1
  • 1
  • 1
    It is a little hard to follow what you are doing and this demonstrates the pitfalls of using meaningless variable names. I would look at all the places where you do something like `if P[random.randint(0,nums-1)]==2: P[random.randint(0,nums-1)]=0` Note that these people are not the same people. so a sick person is *not necessarily* recovering here.. – JonSG Mar 09 '22 at 16:27
  • Just a note expanding on the above comment; code is meant to be readable and understandable by humans. Your code has no helpful comments or variable names. It should be obvious, given comments and variable names, what any given block of code does, and why. It's going to make it much harder for you or others to debug your code when you basically have to run it to figure out what it does or is supposed to do. – Random Davis Mar 09 '22 at 16:33
  • It's common for mathematicians and physicists to use short variable names. If for example the Greek letter ρ is used in a formula, then the variable name `rho` is often used. A programmer prefers a name like `density` (if that's what ρ denotes). – md2perpe Mar 09 '22 at 16:39
  • Run through your code with a [debugger](https://stackoverflow.com/questions/25385173/what-is-a-debugger-and-how-can-it-help-me-diagnose-problems). – Alias Cartellano Mar 09 '22 at 17:53

0 Answers0