0:59
The particular kind of model we use will be SIR.
And it takes a little while to explain it and then we still we still also have to
talk little bit about how the Julia code works but then we can plug the results.
So SRI, the name of the model stands for
three kinds of people involved in the disease.
The first kind are the susceptibles.
These are people who might still get ill.
Initially there are S at node, that time equals zero,
but then S time goes on, this changes.
And so we write S(t) for
the number of people who are susceptible to the disease at time t.
At any given time as the disease progresses,
there maybe some infected people, which we symbolize with i.
And of course initially S(t) equal zero.
There are some people who are ill because otherwise the epidemic would not start.
And then finally people previously infected but
no longer infected called removed.
We assume that they remain non infectious from that time onward.
So this is not a disease like flu where people can get re-infected all the time.
Unfortunately, remove does not mean that they have recovered.
The illness isn't a joke.
So it's important to notice that S, I and R do not have to be whole numbers.
All we want is for S, I, and R to be close to what is observed.
And by close we just mean in some few percent.
If we can get a model that gets within a few percent of the values
that are observed we have done brilliantly well.
Epidemiology is not an exact science.
Getting the numbers within a few percent unfortunately is only possible with
hindsight usually.
We will do this with the SIR model for our able epidemic.
Okay so let's get the equations.
What we do is that we assume that time goes into discrete steps.
And the distance from one time to another,
the duration of time form one step to the next will always be dt.
That means that T sub I which is the time of the I steps is t sub zero.
Which is zero, plus I times dt, so
t sub I is just I times dt.
So code from t sub I to t sub I plus one,
we want to compute values for S of t I plus one.
I at D sub I plus one and RD sub I plus one.
To predict these we will only use the values at the previous time step so
that's s and t sub I and S t sub I and R t sub I plus one.
And they are what we use, so we need a formula to compute S and
t sub I plus one using these.
And similarly for I.
To compute I t sub I plus one using those three.
And so our model is a set of three equations, and they all have this form.
The new value is the old value, plus the gains, minus the losses.
So that's a very obvious form in which models are derived.
And all of the effort of modeling goes into working out how you might compute
at every time step gains and the losses.
And for susceptibles,
of course, the only thing that can happen is that a susceptible becomes infected.
The simplest SIR models use the law of mass action.
And what it says is that all susceptibles are exactly equivalent and
they all have an equal chance of meeting an infected.
Similarly all infecteds are exactly equal,
and they all have an exactly equal chance of meeting a susceptible.
And so the number of meetings between infected and
susceptibles is just proportional to the product S times I.
That's the goal of mass.
Now the number of infections that will result from meetings is not equal to SMI.
It's not equal equal to S times I.
It's just proportional to it.
We will say that this is a constant,
that the probability that it's acceptable is infected when
4:53
an infected person is exactly the same for all people at all times.
And so then, lambda SI is the rate at which susceptibles become infected.
This is a rate per day, because we are using days as our time variable.
And so the actual loss over a time step dt in terms of day.
So say dt is one tenth of a day.
Then actually only a tenth of lambda SI should get infected.
And so lambda SI times dt is the rate at which
the lambda SI times dt is the number of new infections that have generated.
During this time step dt.
The new value is the old value minus the losses.
Now we go on to infected people.
So infected people stop being infected by being removed or being recovered.
And again we say that all infected people have the same chance of recovery
at every day.
So some constant fraction is removed every day.
But this is the fraction gamma.
And so therefore gamma I is the number of infecteds that get removed per day.
And so gamma I dt, that's the loss of infecteds.
The kind of infecteds we already know, because every bit of
population that gets lost from susceptible gets added to infected.
And so here is the new infecteds, the previous level of infected
plus all the ones that are newly infected minus the ones that are now being removed.
Because it's an SIR model, once you're removed you stay there, so
all that happens is that there's a gain term.
And this gain term is of course the same as the last term from the infecteds, and
so that's our third equation.
And our SIR model is just three equations for S, I, and R.
6:43
So now we go on to writing the Julia code.
What I choose to do is to write a single function, which takes S, I and
R all at once, and gives back the new values of S, I and R all at once.
So these functions of course have to know the values of the parameters.
And the way we're going to pause these values to Julia is not to write them
as part of the argument of the function.
But just to declare them at the level where the function is created,
and then the function will see them on the inside of the function.
This is quite surprising.
Many experienced programmers would not expect to see this.
The way it works is the following.
So, if we assign to the variable b the value a, that is no, I'll be here.
Then the function f of x is created so when I run this, that executes.
Then this line executes and it creates a function x.
Now, this point the value of b.
Is actually known.
So when we run this, we get 56 as the outcome.
If we change that to six, of course that becomes 48.
But if we change that to four, then that becomes 24.
So makes our function simple to write.
We have a population vector with three elements.
One, two and three and they will be the susceptibles, the infected, and
the removeds.
And then they'll create the new s is the old susceptibles minus lambda
times the susceptibles times the infected times dt.
8:54
Now we compiles and we set our parameter dt, which is the time state length.
And lambda and
gamma which are the epidemic parameters which specify some sort of input argument.
And recall that this is a way of using one assignment to set three different values.
So we can set S, I, and R.
S becomes 1000, I ten, R becomes 20.
Notice that only 1000 is a float, the other two are integers.
In fact, what will happen is that Julia will convert ten and 20 to integers.
It's not the best form of programming, but Julia forgives us that, and
that's part of why it's easy and quick to write you the app code.
If we want to optimize the running of the Julia code then,
we might find in the end that, this is one of the reasons why our code is running
slower than we might wanted, right.
So, what does it do?
It gives us an output of 975,
this is susceptible as minor lambda times infected.
So, the number ten that become 10,000
then gam is one tenth that means 1000.
Okay so the function works and so we can go back to the top and
now we discuss how we look through all the times.
We want that function to run from start time to the end time and
that can't be infinite so we define an end time t final.
We will use 610 day to fit data from Wikipedia.
So what that means is that there are n steps.
Tf divided by dt, but now you've gotta watch out.
Computer arithmetic is not infinitely precise.
So if we compute n times dt, so we've defined n this way.
And now we calculate Tf as n times dt,
this should be exactly true mathematically, algebraically.
But, in a computer we cannot be sure that this is the case.
So if we wanted to, we tested this is an equality might not work exactly.
So what I would like to do instead is to use time as a variable,
and to keep explicit track of this number I.
Then I know exactly how many steps I want to take, at the end of it,
the time I reach might not be exact the Tf that I started with.
But it will be extremely close to it.
Actually, the take home message is that you should not rely on exact equality,
even if you can reasonably expect it.
Computer code that relies on exact equality is much more prone to error,
and it really requires professional programmers to use.
We know exactly how many values we are going to need,
then we can create an array to hold all of them.
If we take n steps, we need n plus one, that is of course we need a value for
T equals zero and then the values n steps on from there.
So total number of n plus one values.
So here's the concept code, now this code doesn't actually do anything.
The reason it doesn't do anything is that I only have comment lines inside the loop.
This is just to say, okay final 610 days, so
the n steps is an integer, that's why I call the round function.
And round with type integer means that it's going to round it to integer types.
So that gives me an integer number of steps.
And in the results I will hold in a variable, three columns n plus one rows.
And also a vector of times, which will just be n plus one float values.
We need to specify our initial values and our initial result values.
Right.
Now, we are ready to run the model.
So I just cut and paste this code here, initializing everything.
Set all the values we need.
We need to set lambda, we need to set gam, we need to set dt, we need to set tfinal,
we need to set S not, we need to set I not, we need to set R not.
None of these can be taken for granted.
And then we initialize.
We calculate the number of steps.
We initialize a vect.
We initialize on the right, or results and for the time vector,
we initialize this array.
So we've created the array here and then we initialize its first row.
We're not going to compute this row.
We're also not going to compute this time.
Having done that, we can now just take n steps,
now if we want to take a step so we have one.
The next value is going to go to two.
And then if the step is two the next value goes to three and
in the end we end up with the n plus one steps that we put in our array.