1
00:00:01,050 --> 00:00:03,000
In this video, I wanna talk about

2
00:00:03,000 --> 00:00:04,440
computational simulation,

3
00:00:04,440 --> 00:00:06,690
so a code that would do a random

4
00:00:06,690 --> 00:00:09,570
or what we call stochastic simulation

5
00:00:09,570 --> 00:00:11,550
of these discrete dynamical systems.

6
00:00:11,550 --> 00:00:14,220
So usually with our
mathematical description,

7
00:00:14,220 --> 00:00:16,500
we're gonna be following the average state

8
00:00:16,500 --> 00:00:18,840
of possible future for that system.

9
00:00:18,840 --> 00:00:22,950
But we might sometimes be
interested in following the,

10
00:00:22,950 --> 00:00:25,320
so random evolution for some systems

11
00:00:25,320 --> 00:00:28,380
or the full distribution
of possible states.

12
00:00:28,380 --> 00:00:30,720
There are mathematical tools to do that.

13
00:00:30,720 --> 00:00:31,620
So if you're interested,

14
00:00:31,620 --> 00:00:34,050
you might wanna look into master equations

15
00:00:34,050 --> 00:00:36,630
which follow the entire
distribution of possible futures

16
00:00:36,630 --> 00:00:37,473
for a system.

17
00:00:38,520 --> 00:00:40,350
In this course, we're gonna focus on

18
00:00:40,350 --> 00:00:43,380
mean-field approximations
for the mathematical version,

19
00:00:43,380 --> 00:00:46,110
but we're gonna be using
computational simulations,

20
00:00:46,110 --> 00:00:50,043
you know, using code to
describe these systems.

21
00:00:51,510 --> 00:00:53,040
So to make it a little simpler,

22
00:00:53,040 --> 00:00:55,890
because I'm gonna have to
write pseudocode in this video,

23
00:00:57,976 --> 00:01:00,150
I wanna go back to a simpler example.

24
00:01:00,150 --> 00:01:04,450
So we'll be going back to the
example from pharmacokinetic

25
00:01:05,826 --> 00:01:07,440
of drugs going into an organ.

26
00:01:07,440 --> 00:01:12,300
And I remember last week I
phrased it as a social system,

27
00:01:12,300 --> 00:01:14,730
so people entering and leaving a building.

28
00:01:14,730 --> 00:01:16,440
The reason why I do this is just that

29
00:01:16,440 --> 00:01:20,180
I tend to anthropomorphize all models.

30
00:01:22,590 --> 00:01:24,900
So I often, like, put myself
in the shoes of the parts

31
00:01:24,900 --> 00:01:27,060
when I'm trying to think
through the mechanism

32
00:01:27,060 --> 00:01:29,960
and then it's easier for me
to use the language of people.

33
00:01:30,930 --> 00:01:33,483
So, if you remember,

34
00:01:35,610 --> 00:01:38,610
in the simple model we
have one compartment

35
00:01:38,610 --> 00:01:42,010
which represents the population N of t

36
00:01:42,930 --> 00:01:44,883
of people currently in a building.

37
00:01:46,560 --> 00:01:50,670
We have one arrow for people
entering that building.

38
00:01:50,670 --> 00:01:54,630
So let's say that they arrive at rate a,

39
00:01:54,630 --> 00:01:56,730
and to use the notation
from the previous video,

40
00:01:56,730 --> 00:02:00,710
that means that we go from
nothing and create one person.

41
00:02:03,120 --> 00:02:05,793
So we add one to N of t at this rate a.

42
00:02:08,820 --> 00:02:11,700
And then we might have some
other rate at which they leave,

43
00:02:11,700 --> 00:02:16,023
where we go from a person to nothing.

44
00:02:17,610 --> 00:02:20,130
So that's a very simple model, right?

45
00:02:20,130 --> 00:02:25,130
And in a discrete system,

46
00:02:25,140 --> 00:02:28,500
we might still have continuous
value for a and l, right?

47
00:02:28,500 --> 00:02:32,070
So maybe if a is the
average number of person

48
00:02:32,070 --> 00:02:34,323
that enter that building per week.

49
00:02:36,870 --> 00:02:39,450
Well, it doesn't have to be discrete,

50
00:02:39,450 --> 00:02:41,490
even in a discrete system,

51
00:02:41,490 --> 00:02:44,587
'cause if we measure it
over time and we say,

52
00:02:44,587 --> 00:02:47,040
"Well, this is not a very popular business

53
00:02:47,040 --> 00:02:48,250
and there's only 1.5

54
00:02:51,810 --> 00:02:54,100
people who enter

55
00:02:58,830 --> 00:03:02,610
per week," right?

56
00:03:02,610 --> 00:03:04,560
So that's my arrival rate, a.

57
00:03:04,560 --> 00:03:05,430
Well, that's possible.

58
00:03:05,430 --> 00:03:08,670
It's discrete data but divided
by some amount of time,

59
00:03:08,670 --> 00:03:11,193
so I can get fractional people.

60
00:03:13,320 --> 00:03:15,750
And then I might wanna
think of it in terms of day,

61
00:03:15,750 --> 00:03:18,720
maybe my model, I only
wanna think about, you know,

62
00:03:18,720 --> 00:03:20,910
a morning or an afternoon
in that building.

63
00:03:20,910 --> 00:03:24,270
So I could re-normalize
this to any amount of time.

64
00:03:24,270 --> 00:03:26,020
So really I could say that this is

65
00:03:34,350 --> 00:03:36,600
1.5 people who enter per week.

66
00:03:36,600 --> 00:03:40,893
So it's 1.5 over seven,
let's say, people per day.

67
00:03:46,200 --> 00:03:50,163
So that's equal to 3/2,

68
00:03:51,210 --> 00:03:55,863
3/14 people per day.

69
00:03:57,240 --> 00:03:58,710
Trying to save space here.

70
00:03:58,710 --> 00:04:00,140
And that's roughly...

71
00:04:01,410 --> 00:04:02,970
Let me just...

72
00:04:02,970 --> 00:04:07,133
That's roughly 20% or 21%.

73
00:04:08,400 --> 00:04:12,390
It's probably like 21.5, 21.4.

74
00:04:12,390 --> 00:04:14,583
Let's do 20% approximate.

75
00:04:15,540 --> 00:04:16,950
Okay, so what that means

76
00:04:16,950 --> 00:04:20,430
is that if my increments of
times are days, let's say,

77
00:04:20,430 --> 00:04:22,530
in my discrete time model,

78
00:04:22,530 --> 00:04:26,070
on any given day, I have 21% chance

79
00:04:26,070 --> 00:04:28,023
of seeing someone enter, right?

80
00:04:29,250 --> 00:04:32,250
In these models, you can always
re-normalize time like that.

81
00:04:33,240 --> 00:04:35,580
And the reason I'm pointing it out here

82
00:04:35,580 --> 00:04:38,520
is that 1.5 people per week,

83
00:04:38,520 --> 00:04:39,990
well if our increments are week,

84
00:04:39,990 --> 00:04:41,850
well that's sort of annoying

85
00:04:41,850 --> 00:04:43,650
because what we wanna be able to say is,

86
00:04:43,650 --> 00:04:46,233
is there someone who enters, yes or no?

87
00:04:47,490 --> 00:04:50,070
And then, often, re-normalizing
it to smaller rate

88
00:04:50,070 --> 00:04:53,070
just makes it a little easier
to think computationally.

89
00:04:53,070 --> 00:04:56,490
So here, we can see that
if my increments of time

90
00:04:56,490 --> 00:04:58,440
my plus one in time is a day,

91
00:04:58,440 --> 00:05:01,713
well during that day, there is
a 21% chance that they enter.

92
00:05:02,550 --> 00:05:06,090
And the intuition I wanna give you here

93
00:05:06,090 --> 00:05:11,090
is how do you introduce
randomness like that

94
00:05:11,640 --> 00:05:14,400
yes or no, does someone enter on that day

95
00:05:14,400 --> 00:05:15,933
in a computer code?

96
00:05:18,480 --> 00:05:19,970
So on a given day...

97
00:05:27,180 --> 00:05:28,770
Well, throughout this course,

98
00:05:28,770 --> 00:05:31,950
the only thing I'm gonna
assume we know how to do

99
00:05:31,950 --> 00:05:35,350
is how to get a random
number between zero and one

100
00:05:37,380 --> 00:05:38,313
in a computer.

101
00:05:39,420 --> 00:05:42,960
And you know, there's a lot
of literature behind this

102
00:05:42,960 --> 00:05:45,360
because, you know, computers
are never generating

103
00:05:45,360 --> 00:05:46,650
truly random number,

104
00:05:46,650 --> 00:05:50,130
so we call them pseudo-random
number generators.

105
00:05:50,130 --> 00:05:51,720
We're just gonna take that for granted.

106
00:05:51,720 --> 00:05:53,820
If you wanna have a discussion about that,

107
00:05:53,820 --> 00:05:57,600
we can meet on Wednesday
morning during student hours,

108
00:05:57,600 --> 00:05:59,670
but right now we're just
gonna take that for granted.

109
00:05:59,670 --> 00:06:01,370
So the one thing we know how to do

110
00:06:03,900 --> 00:06:07,590
is generate number on the
interval between zero and one.

111
00:06:07,590 --> 00:06:10,683
So if I call a function, say, RAND1,

112
00:06:11,730 --> 00:06:14,430
which I think is probably
a thing in MATLAB

113
00:06:14,430 --> 00:06:16,380
to get one number between zero and one,

114
00:06:17,400 --> 00:06:21,960
I call it, I might get, you know, 17%.

115
00:06:21,960 --> 00:06:25,193
I call it again, I might get 90%.

116
00:06:26,460 --> 00:06:30,040
And I can call it again and again

117
00:06:30,900 --> 00:06:33,840
and sort of fill this line
with a uniform distribution.

118
00:06:33,840 --> 00:06:38,343
So I'm uniformly sampling
numbers between zero and one.

119
00:06:39,720 --> 00:06:41,580
Okay, so on a given day,

120
00:06:41,580 --> 00:06:44,193
I may able to do this,
generate one number.

121
00:06:45,120 --> 00:06:46,980
Well, if what you wanna do

122
00:06:46,980 --> 00:06:51,090
is let people enter with probability 21%,

123
00:06:51,090 --> 00:06:52,537
you can simply say,

124
00:06:52,537 --> 00:06:57,537
"If on that day the random
number that I generate

125
00:06:58,050 --> 00:07:01,140
is smaller than 21%,"

126
00:07:01,140 --> 00:07:03,600
which is only gonna happen, right,

127
00:07:03,600 --> 00:07:06,513
so 21% is some threshold value here,

128
00:07:09,240 --> 00:07:11,700
well that's only gonna
happen 21% of the time.

129
00:07:11,700 --> 00:07:14,880
21% of the time I get a
number that's smaller,

130
00:07:14,880 --> 00:07:17,070
so to the left of my red line,

131
00:07:17,070 --> 00:07:21,390
and 79% of the time I get a
number that's higher, right?

132
00:07:21,390 --> 00:07:25,923
So that's 79%, 21%, if
I did my math right.

133
00:07:27,390 --> 00:07:32,390
So if random, this random
number, is below 21%,

134
00:07:32,910 --> 00:07:37,910
then I can say N of t goes
to N of t plus one, right?

135
00:07:41,490 --> 00:07:45,900
And then if my rate of departure,

136
00:07:45,900 --> 00:07:48,120
the rate at which they leave, l,

137
00:07:48,120 --> 00:07:50,580
is, you know, on an average day,

138
00:07:50,580 --> 00:07:52,170
every person that's in the building...

139
00:07:52,170 --> 00:07:53,850
I don't know which kind
of building this is,

140
00:07:53,850 --> 00:07:56,790
maybe some very non-popular hotel.

141
00:07:56,790 --> 00:07:59,670
On a given day, everyone
that is in the building

142
00:07:59,670 --> 00:08:01,533
leaves with some probability,

143
00:08:02,970 --> 00:08:06,053
let's call it l equals to 10%.

144
00:08:11,220 --> 00:08:13,470
Then I would just have a separate loop

145
00:08:13,470 --> 00:08:15,480
that goes over every individual.

146
00:08:15,480 --> 00:08:20,400
And now with probability 10%,

147
00:08:20,400 --> 00:08:23,280
all of these individuals would leave

148
00:08:23,280 --> 00:08:28,280
and I would have N minus
one to my population, right?

149
00:08:28,680 --> 00:08:30,210
So that's sort of the intuition.

150
00:08:30,210 --> 00:08:32,430
Once we have rates, it's really easy

151
00:08:32,430 --> 00:08:34,260
to go over the parts in the system

152
00:08:34,260 --> 00:08:37,413
and decide which transition
actually do happen.

153
00:08:41,880 --> 00:08:42,990
So it sounds really easy,

154
00:08:42,990 --> 00:08:45,780
but there's a few problem that can happen.

155
00:08:45,780 --> 00:08:48,480
So let's say, let's make it
a little more interesting.

156
00:08:51,960 --> 00:08:53,283
So we have our building,

157
00:08:55,890 --> 00:08:59,610
and then we say that people enter

158
00:08:59,610 --> 00:09:02,670
or arrive at some number a, some rate a,

159
00:09:02,670 --> 00:09:04,293
and they leave at some rate l,

160
00:09:05,160 --> 00:09:06,690
but now let's say that this building

161
00:09:06,690 --> 00:09:10,563
that's not super popular is
actually a clone factory.

162
00:09:14,310 --> 00:09:16,140
And I'm saying that just
because I wanna introduce

163
00:09:16,140 --> 00:09:17,280
a new mechanism.

164
00:09:17,280 --> 00:09:20,310
So while in this building,
people can clone themselves,

165
00:09:20,310 --> 00:09:22,110
and I'm gonna write it like that.

166
00:09:22,110 --> 00:09:27,000
So here, we went from, you
know, nothing to something,

167
00:09:27,000 --> 00:09:29,700
which maybe I should write p for person

168
00:09:29,700 --> 00:09:32,130
since my variable here is N of t,

169
00:09:32,130 --> 00:09:34,080
but hopefully that's not too confusing.

170
00:09:34,080 --> 00:09:38,400
Here, we go from a
contribution to N to nothing.

171
00:09:38,400 --> 00:09:43,400
Well, here we go from
one person to two people

172
00:09:43,680 --> 00:09:45,390
at a rate c, right?

173
00:09:45,390 --> 00:09:48,240
So people, while they're
in this clone factory,

174
00:09:48,240 --> 00:09:51,330
they get information and
they get a marketing pitch

175
00:09:51,330 --> 00:09:52,830
about why you should clone yourself,

176
00:09:52,830 --> 00:09:54,180
and with some rate c,

177
00:09:54,180 --> 00:09:56,730
they're gonna choose to clone themselves.

178
00:09:56,730 --> 00:09:58,920
Okay, so let's say I
wanna model, you know,

179
00:09:58,920 --> 00:10:01,323
a week in this weird clone factory.

180
00:10:03,480 --> 00:10:07,680
So I'm gonna say for every day...

181
00:10:07,680 --> 00:10:10,440
I'm gonna stick to MATLAB code for now.

182
00:10:10,440 --> 00:10:14,832
So t is gonna be the label of that day,

183
00:10:14,832 --> 00:10:18,210
which is gonna go from day
one in increments of delta t,

184
00:10:18,210 --> 00:10:19,140
just to be general;

185
00:10:19,140 --> 00:10:21,240
delta t here is one day;

186
00:10:21,240 --> 00:10:23,820
up to maximum, capital T.

187
00:10:23,820 --> 00:10:26,703
Capital T would be seven
days for my week, right?

188
00:10:27,690 --> 00:10:32,690
So for every day, I'm gonna
ask, draw a random number.

189
00:10:35,070 --> 00:10:38,973
If that random number is
smaller than a delta t,

190
00:10:40,710 --> 00:10:45,710
so the rate per day at
which people arrive,

191
00:10:45,960 --> 00:10:48,150
then my N of t,

192
00:10:48,150 --> 00:10:49,910
which I would write like that...

193
00:10:50,942 --> 00:10:51,930
I guess here I would
have a little product.

194
00:10:51,930 --> 00:10:52,980
Trying to write MATLAB,

195
00:10:52,980 --> 00:10:55,290
but let's not focus too much on that.

196
00:10:55,290 --> 00:10:58,260
And then here I would
go from N of t plus one,

197
00:10:58,260 --> 00:11:00,240
so someone just arrived.

198
00:11:00,240 --> 00:11:02,690
And I think you have to
end your loops in MATLAB.

199
00:11:04,590 --> 00:11:06,213
Cool. So that's my arrivals.

200
00:11:12,120 --> 00:11:14,700
Okay. What other mechanism do I have?

201
00:11:14,700 --> 00:11:16,230
That's the first arrow.

202
00:11:16,230 --> 00:11:18,363
The other arrow is the departure.

203
00:11:20,310 --> 00:11:22,650
So here, the departures would be,

204
00:11:22,650 --> 00:11:25,503
for everyone, so let's call it small n,

205
00:11:26,670 --> 00:11:31,593
so we've labeled people with
a number from one to N_t.

206
00:11:32,760 --> 00:11:34,350
And all of these nice peoples,

207
00:11:34,350 --> 00:11:36,270
well they get their marketing pitch

208
00:11:36,270 --> 00:11:38,790
but with some probability,

209
00:11:38,790 --> 00:11:40,860
and I get that by drawing a random number,

210
00:11:40,860 --> 00:11:45,720
some probability l delta t,
they're gonna choose to leave.

211
00:11:45,720 --> 00:11:46,803
And if they leave,

212
00:11:51,810 --> 00:11:56,703
I'll have to decrease this
total number of people.

213
00:11:57,600 --> 00:11:59,322
You see, like this is a pseudocode.

214
00:11:59,322 --> 00:12:02,820
Here, I would have, you know,
n of next day or something

215
00:12:02,820 --> 00:12:07,650
to avoid having any like
interference between my update here

216
00:12:07,650 --> 00:12:09,870
and my end value here,

217
00:12:09,870 --> 00:12:12,020
but that doesn't matter
too much right now.

218
00:12:13,260 --> 00:12:15,142
Let's just call it pseudocode.

219
00:12:15,142 --> 00:12:16,392
And these are departures.

220
00:12:22,200 --> 00:12:23,870
All right, how do I model

221
00:12:23,870 --> 00:12:28,200
or how do I code this new mechanism now?

222
00:12:28,200 --> 00:12:31,860
I end my departures and I got
a third arrow that I added

223
00:12:31,860 --> 00:12:33,633
for people cloning themselves.

224
00:12:35,910 --> 00:12:37,920
Okay, well, actually I don't need to close

225
00:12:37,920 --> 00:12:39,060
my departure right now

226
00:12:39,060 --> 00:12:44,060
because everyone in my loop
here, n equal one to N_t,

227
00:12:44,430 --> 00:12:48,990
everyone here, this is everyone,

228
00:12:48,990 --> 00:12:51,390
can choose to clone themselves, right?

229
00:12:51,390 --> 00:12:56,390
And what we say is that a similar process

230
00:12:58,200 --> 00:13:01,053
and with probability c delta t now,

231
00:13:02,310 --> 00:13:04,200
they're gonna choose to clone themselves.

232
00:13:04,200 --> 00:13:08,793
And cloning, you know,
simply adds someone.

233
00:13:10,200 --> 00:13:11,160
End. End.

234
00:13:12,240 --> 00:13:16,050
Okay, so now we have some
sort of a linear mechanism

235
00:13:16,050 --> 00:13:18,120
for how people appear.

236
00:13:18,120 --> 00:13:19,680
And I could like close this here

237
00:13:19,680 --> 00:13:22,230
and I would have this random week

238
00:13:22,230 --> 00:13:24,123
in the life of this cloning factory.

239
00:13:25,140 --> 00:13:26,140
So this was cloning.

240
00:13:28,440 --> 00:13:29,273
Okay.

241
00:13:29,273 --> 00:13:32,430
Well, there's a few problems
you can see occurring,

242
00:13:32,430 --> 00:13:34,770
other than the possible interference

243
00:13:34,770 --> 00:13:37,320
because of the name of the variables.

244
00:13:37,320 --> 00:13:38,470
Oops, sorry about that.

245
00:13:40,350 --> 00:13:42,450
Well, one of the problem is, you know,

246
00:13:42,450 --> 00:13:47,450
what if one of these rate is
such that, say, the arrival...

247
00:13:48,930 --> 00:13:51,420
Right, let's say the
company gets real popular,

248
00:13:51,420 --> 00:13:56,420
what happens if a delta
t is greater than one?

249
00:13:57,630 --> 00:14:00,270
If, on average, in any one
day I expect 10 people,

250
00:14:00,270 --> 00:14:03,210
here I'm still incrementing
them by one every day

251
00:14:03,210 --> 00:14:05,130
because for sure my random number

252
00:14:05,130 --> 00:14:07,503
is gonna be greater than one here.

253
00:14:08,790 --> 00:14:10,980
But how do I actually get 10 people?

254
00:14:10,980 --> 00:14:13,320
Well, that's why I had
this general delta t here,

255
00:14:13,320 --> 00:14:15,600
because we might wanna
re-normalize our rates

256
00:14:15,600 --> 00:14:17,550
so that this doesn't happen, right?

257
00:14:17,550 --> 00:14:19,550
That's one thing we have to think about.

258
00:14:20,610 --> 00:14:22,000
The other thing here

259
00:14:23,430 --> 00:14:26,640
is in my departure and cloning mechanism,

260
00:14:26,640 --> 00:14:28,113
is that order matters.

261
00:14:33,060 --> 00:14:36,270
It matters because I'm
testing departure first.

262
00:14:36,270 --> 00:14:40,170
And anyone who left, who chooses to leave

263
00:14:40,170 --> 00:14:43,620
shouldn't be allowed to
clone themselves, right?

264
00:14:43,620 --> 00:14:48,620
So really what I should have here is else.

265
00:14:49,200 --> 00:14:50,033
Oops.

266
00:14:55,080 --> 00:14:57,810
If you don't leave, you can
choose to clone yourself.

267
00:14:57,810 --> 00:14:59,730
But even then,

268
00:14:59,730 --> 00:15:02,520
the fact that you're
testing for departure first

269
00:15:02,520 --> 00:15:05,220
means that you're less
likely than you would want

270
00:15:05,220 --> 00:15:06,483
to clone yourself.

271
00:15:08,340 --> 00:15:09,720
Even worse, in the previous case,

272
00:15:09,720 --> 00:15:11,490
you could have like
people cloning themselves

273
00:15:11,490 --> 00:15:12,650
and departing at the same time.

274
00:15:12,650 --> 00:15:15,873
So you can have these, like,
simultaneous events, right?

275
00:15:17,130 --> 00:15:21,150
So I think I'm gonna pause this video now.

276
00:15:21,150 --> 00:15:23,370
I'm gonna end it on these problems,

277
00:15:23,370 --> 00:15:27,330
really to see how to make sure that, like,

278
00:15:27,330 --> 00:15:31,560
the order in which we do
our different mechanisms

279
00:15:31,560 --> 00:15:33,000
doesn't matter

280
00:15:33,000 --> 00:15:36,000
and to make sure that we
don't have any problem

281
00:15:36,000 --> 00:15:37,200
dealing with high rates,

282
00:15:37,200 --> 00:15:39,240
so probabilities greater than one.

283
00:15:39,240 --> 00:15:42,450
We really need to better
understand the statistics

284
00:15:42,450 --> 00:15:44,190
behind these process.

285
00:15:44,190 --> 00:15:46,200
So I'm just gonna drop a word here.

286
00:15:46,200 --> 00:15:48,030
Everything that we use here,

287
00:15:48,030 --> 00:15:51,723
all of these mechanisms are
what we call Poisson process.

288
00:15:52,710 --> 00:15:54,000
And these Poisson processes

289
00:15:54,000 --> 00:15:56,760
simply mean that something
happens at a given rate

290
00:15:56,760 --> 00:15:59,040
regardless of what happened beforehand.

291
00:15:59,040 --> 00:16:00,780
And we really need to better understand

292
00:16:00,780 --> 00:16:02,850
these Poisson processes
if we wanna make sure

293
00:16:02,850 --> 00:16:04,950
that we don't make any
mistake with our code.

294
00:16:04,950 --> 00:16:07,100
So I'm gonna end this video here.

295
00:16:07,100 --> 00:16:09,570
In the next one, we're gonna
look at the mathematics

296
00:16:09,570 --> 00:16:12,210
of Poisson processes and
then discuss what that means

297
00:16:12,210 --> 00:16:15,150
for our computational
implementation, all right?

298
00:16:15,150 --> 00:16:16,600
I'll see you on the next one.