In this video, I wanna talk about
computational simulation,
so a code that would do a random
or what we call stochastic simulation
of these discrete dynamical systems.
So usually with our mathematical description,
we're gonna be following the average state
of possible future for that system.
But we might sometimes be interested in following the,
so random evolution for some systems
or the full distribution of possible states.
There are mathematical tools to do that.
So if you're interested,
you might wanna look into master equations
which follow the entire distribution of possible futures
for a system.
In this course, we're gonna focus on
mean-field approximations for the mathematical version,
but we're gonna be using computational simulations,
you know, using code to describe these systems.
So to make it a little simpler,
because I'm gonna have to write pseudocode in this video,
I wanna go back to a simpler example.
So we'll be going back to the example from pharmacokinetic
of drugs going into an organ.
And I remember last week I phrased it as a social system,
so people entering and leaving a building.
The reason why I do this is just that
I tend to anthropomorphize all models.
So I often, like, put myself in the shoes of the parts
when I'm trying to think through the mechanism
and then it's easier for me to use the language of people.
So, if you remember,
in the simple model we have one compartment
which represents the population N of t
of people currently in a building.
We have one arrow for people entering that building.
So let's say that they arrive at rate a,
and to use the notation from the previous video,
that means that we go from nothing and create one person.
So we add one to N of t at this rate a.
And then we might have some other rate at which they leave,
where we go from a person to nothing.
So that's a very simple model, right?
And in a discrete system,
we might still have continuous value for a and l, right?
So maybe if a is the average number of person
that enter that building per week.
Well, it doesn't have to be discrete,
even in a discrete system,
'cause if we measure it over time and we say,
"Well, this is not a very popular business
and there's only 1.5
people who enter
per week," right?
So that's my arrival rate, a.
Well, that's possible.
It's discrete data but divided by some amount of time,
so I can get fractional people.
And then I might wanna think of it in terms of day,
maybe my model, I only wanna think about, you know,
a morning or an afternoon in that building.
So I could re-normalize this to any amount of time.
So really I could say that this is
1.5 people who enter per week.
So it's 1.5 over seven, let's say, people per day.
So that's equal to 3/2,
3/14 people per day.
Trying to save space here.
And that's roughly...
Let me just...
That's roughly 20% or 21%.
It's probably like 21.5, 21.4.
Let's do 20% approximate.
Okay, so what that means
is that if my increments of times are days, let's say,
in my discrete time model,
on any given day, I have 21% chance
of seeing someone enter, right?
In these models, you can always re-normalize time like that.
And the reason I'm pointing it out here
is that 1.5 people per week,
well if our increments are week,
well that's sort of annoying
because what we wanna be able to say is,
is there someone who enters, yes or no?
And then, often, re-normalizing it to smaller rate
just makes it a little easier to think computationally.
So here, we can see that if my increments of time
my plus one in time is a day,
well during that day, there is a 21% chance that they enter.
And the intuition I wanna give you here
is how do you introduce randomness like that
yes or no, does someone enter on that day
in a computer code?
So on a given day...
Well, throughout this course,
the only thing I'm gonna assume we know how to do
is how to get a random number between zero and one
in a computer.
And you know, there's a lot of literature behind this
because, you know, computers are never generating
truly random number,
so we call them pseudo-random number generators.
We're just gonna take that for granted.
If you wanna have a discussion about that,
we can meet on Wednesday morning during student hours,
but right now we're just gonna take that for granted.
So the one thing we know how to do
is generate number on the interval between zero and one.
So if I call a function, say, RAND1,
which I think is probably a thing in MATLAB
to get one number between zero and one,
I call it, I might get, you know, 17%.
I call it again, I might get 90%.
And I can call it again and again
and sort of fill this line with a uniform distribution.
So I'm uniformly sampling numbers between zero and one.
Okay, so on a given day,
I may able to do this, generate one number.
Well, if what you wanna do
is let people enter with probability 21%,
you can simply say,
"If on that day the random number that I generate
is smaller than 21%,"
which is only gonna happen, right,
so 21% is some threshold value here,
well that's only gonna happen 21% of the time.
21% of the time I get a number that's smaller,
so to the left of my red line,
and 79% of the time I get a number that's higher, right?
So that's 79%, 21%, if I did my math right.
So if random, this random number, is below 21%,
then I can say N of t goes to N of t plus one, right?
And then if my rate of departure,
the rate at which they leave, l,
is, you know, on an average day,
every person that's in the building...
I don't know which kind of building this is,
maybe some very non-popular hotel.
On a given day, everyone that is in the building
leaves with some probability,
let's call it l equals to 10%.
Then I would just have a separate loop
that goes over every individual.
And now with probability 10%,
all of these individuals would leave
and I would have N minus one to my population, right?
So that's sort of the intuition.
Once we have rates, it's really easy
to go over the parts in the system
and decide which transition actually do happen.
So it sounds really easy,
but there's a few problem that can happen.
So let's say, let's make it a little more interesting.
So we have our building,
and then we say that people enter
or arrive at some number a, some rate a,
and they leave at some rate l,
but now let's say that this building
that's not super popular is actually a clone factory.
And I'm saying that just because I wanna introduce
a new mechanism.
So while in this building, people can clone themselves,
and I'm gonna write it like that.
So here, we went from, you know, nothing to something,
which maybe I should write p for person
since my variable here is N of t,
but hopefully that's not too confusing.
Here, we go from a contribution to N to nothing.
Well, here we go from one person to two people
at a rate c, right?
So people, while they're in this clone factory,
they get information and they get a marketing pitch
about why you should clone yourself,
and with some rate c,
they're gonna choose to clone themselves.
Okay, so let's say I wanna model, you know,
a week in this weird clone factory.
So I'm gonna say for every day...
I'm gonna stick to MATLAB code for now.
So t is gonna be the label of that day,
which is gonna go from day one in increments of delta t,
just to be general;
delta t here is one day;
up to maximum, capital T.
Capital T would be seven days for my week, right?
So for every day, I'm gonna ask, draw a random number.
If that random number is smaller than a delta t,
so the rate per day at which people arrive,
then my N of t,
which I would write like that...
I guess here I would have a little product.
Trying to write MATLAB,
but let's not focus too much on that.
And then here I would go from N of t plus one,
so someone just arrived.
And I think you have to end your loops in MATLAB.
Cool. So that's my arrivals.
Okay. What other mechanism do I have?
That's the first arrow.
The other arrow is the departure.
So here, the departures would be,
for everyone, so let's call it small n,
so we've labeled people with a number from one to N_t.
And all of these nice peoples,
well they get their marketing pitch
but with some probability,
and I get that by drawing a random number,
some probability l delta t, they're gonna choose to leave.
And if they leave,
I'll have to decrease this total number of people.
You see, like this is a pseudocode.
Here, I would have, you know, n of next day or something
to avoid having any like interference between my update here
and my end value here,
but that doesn't matter too much right now.
Let's just call it pseudocode.
And these are departures.
All right, how do I model
or how do I code this new mechanism now?
I end my departures and I got a third arrow that I added
for people cloning themselves.
Okay, well, actually I don't need to close
my departure right now
because everyone in my loop here, n equal one to N_t,
everyone here, this is everyone,
can choose to clone themselves, right?
And what we say is that a similar process
and with probability c delta t now,
they're gonna choose to clone themselves.
And cloning, you know, simply adds someone.
End. End.
Okay, so now we have some sort of a linear mechanism
for how people appear.
And I could like close this here
and I would have this random week
in the life of this cloning factory.
So this was cloning.
Okay.
Well, there's a few problems you can see occurring,
other than the possible interference
because of the name of the variables.
Oops, sorry about that.
Well, one of the problem is, you know,
what if one of these rate is such that, say, the arrival...
Right, let's say the company gets real popular,
what happens if a delta t is greater than one?
If, on average, in any one day I expect 10 people,
here I'm still incrementing them by one every day
because for sure my random number
is gonna be greater than one here.
But how do I actually get 10 people?
Well, that's why I had this general delta t here,
because we might wanna re-normalize our rates
so that this doesn't happen, right?
That's one thing we have to think about.
The other thing here
is in my departure and cloning mechanism,
is that order matters.
It matters because I'm testing departure first.
And anyone who left, who chooses to leave
shouldn't be allowed to clone themselves, right?
So really what I should have here is else.
Oops.
If you don't leave, you can choose to clone yourself.
But even then,
the fact that you're testing for departure first
means that you're less likely than you would want
to clone yourself.
Even worse, in the previous case,
you could have like people cloning themselves
and departing at the same time.
So you can have these, like, simultaneous events, right?
So I think I'm gonna pause this video now.
I'm gonna end it on these problems,
really to see how to make sure that, like,
the order in which we do our different mechanisms
doesn't matter
and to make sure that we don't have any problem
dealing with high rates,
so probabilities greater than one.
We really need to better understand the statistics
behind these process.
So I'm just gonna drop a word here.
Everything that we use here,
all of these mechanisms are what we call Poisson process.
And these Poisson processes
simply mean that something happens at a given rate
regardless of what happened beforehand.
And we really need to better understand
these Poisson processes if we wanna make sure
that we don't make any mistake with our code.
So I'm gonna end this video here.
In the next one, we're gonna look at the mathematics
of Poisson processes and then discuss what that means
for our computational implementation, all right?
I'll see you on the next one.