MA3264 Mathematical Modelling Lecture 6

1
MA3264 Mathematical ModellingLecture 6

• Monte Carlo Simulation

2
Motivation
Elevator Service Should we
assign elevators to even / odd floors?
use express elevators?
Traffic Control Where should we locate
one way streets?
traffic lights?
Experimental Method May be disruptive since
upset customers may take stairs!
frustrated drivers may take sidewalks!
Therefore decision makers require simulation.
3
Motivation
Simulation can be
Analogue wind tunnels
Digital computes using mathematical models
Mathematical Models
Deterministic formuli such as
Stochastic Markov Chains, etc.
Computations
Analytic directly compute variables or
probabilities
Problem very often intractible
Monte Carlo uses random number generator
4
Random Number Generation
gtgt m4 n 2 gtgt rand(m,n) ans 0.9501
0.8913 0.2311 0.7621 0.6068 0.4565
0.4860 0.0185

• gives m x n matrix whose entries are
• independent random numbers uniformly
• distributed over the interval 0,1

gtgt x rand(10,1) gtgt x rand(10,1) gtgt y
zeros(10,1) gtgt plot(x,y,'r') gtgt hold on gtgt x
rand(10,1) gtgt plot(x,y,'bo')
5
Random Number Generation
gtgt x rand(10,1) gtgt y rand(10,1) gtgt
plot(x,y,'r') gtgt hold on gtgt x rand(10,1) gtgt y
rand(10,1) gtgt plot(x,y,'bo')
gtgt x rand(100,1) gtgt y rand(100,1) gtgt
plot(x,y,'r')
6
Random Number Generation
7
Random Number Generation
8
Area Under a Curve
gtgt s 0.0011 gtgt fs s.2 gtgt
plot(s,fs,'ko') gtgt title('f(s) s2') gtgt count
0 gtgt for j 110000 if y(j) lt x(j)2 count
count 1 end end gtgt area_estimate
count/10000 area_estimate 0.3387
Question What is the area under the curve f(s)
s2 ?
9
Area Under a Curve
gtgt for m 1100 z rand(1000000,2) count
sum(z(,2) lt z(,1).2) area_estimate(m)
count/1000000 end
gtgt z rand(1000000,2) gtgt count sum(z(,2) lt
z(,1).2) gtgt area_estimate count/1000000 area_
estimate 0.3329 gtgt z rand(1000000,2) gtgt
count sum(z(,2) lt z(,1).2) gtgt area_estimate
count/1000000 area_estimate 0.3330 gtgt z
rand(1000000,2) gtgt count sum(z(,2) lt
z(,1).2) gtgt area_estimate count/1000000 area_
estimate 0.3328
Question How does accuracy depend on points
used ?
10
Volume Inside of a Solid Object
Problem Compute the volume of a solid that is
formed by intersecting a ball of radius 1 and a
cube of diameter 3/2 that has the same center as
the ball?
Solution Clearly the volume of the cube 27/8
3.750 Therefore is suffices to compute the
fraction of random points chosen inside the cube
that are inside the sphere. The following MATLAB
commands do the computation
gtgt N 10000000 gtgt count 0 gtgt xyz
1.5rand(N,3)-0.75 gtgt radii sqrt(xyz(,1).2
xyz(,2).2 xyz(,3).2) gtgt count sum(radii
lt 1) gtgt fraction count/10000000 fraction
0.9212 gtgt volume_estimate (27/8)fraction volume
_estimate 3.1092
11
Generating Random Numbers
Middle Square Method, invented by Ulam,
Metropolis
http//en.wikipedia.org/wiki/John_von_Neumann
1. Start with a four digit number , called
the seed.
2. Square it to obtain an eight digit number (add
leading zeros if necessary).
3. Take the middle four digits as the next random
number.
3084 5110 1121 2566 5843 1406 9768 4138 1230
5129 3066 4003 240 576 3317 24 5 0
gtgt x(1) 3084 gtgt for k 117 s x(k)2 s6
mod(s,1000000) r6 s6 - mod(s6,100) x(k1)
r6/100 end
Question Whats wrong?
Question Where to go?
12
Generating Random Numbers
gtgt x(1) 962947 gtgt for k 1100000 s
x(k)2 s9 mod(s,1000000000) r9 s9 -
mod(s9,1000) x(k1) r9/1000 end hist(x/1000000
)
Question What did we do?
gtgt x rand(100000,1) gtgt hist(x)
Question What method is better?
13
Generating Random Numbers
Linear Congruence Method, invented by
http//en.wikipedia.org/wiki/Derrick_Henry_Lehmer
1. Start with a four digit number , called
the seed.
2. Apply the transformation
Here a, b, c are three required positive
integers. The random integers are between 0 and
c-1.
gtgt a 15625 gtgt b 22221 gtgt c 231 gtgt x(1)
389274 gtgt for k 1100000 x(k1)mod(ax(k)b,
c) end x x/(c-1) hist(x)
14
Generating Random Numbers
It is often required to compute random numbers
that are not uniformly distributed over an
interval.
Example When a dice is rolled it produces a
random number in the set 1,2,3,4,5,6, it is
uniformly distributed over that (discrete) set of
points. The following MATLAB commands simulates
100 throws of a dice and computes how many times
each value occured.
gtgt for k 1100 gtgt y rand gtgt x(k)
round(6y1/2) gtgt end gtgt for val 16 gtgt
numbers(val) sum(x val) gtgt end gtgt
numbers numbers 17 18 19 18 17 11
15
Generating Random Numbers
Queueing theory models elevator (lift) and bus
service.
http//en.wikipedia.org/wiki/Queueing_theory
The time between people consecutively arriving
at an elevator or a bus stop is a random
variable that is exponentially distributed over
the interval
This means that
where
is the expected value of
Question How can we generate this random variable
using (only) a uniform random number generator?
16
Generating Random Numbers
Theorem If c gt 0 and y is uniformly distributed
over 0,1 then x - c ln(y) is exponentially
distributed over
Proof
The MATLAB code generates 100000 random samples
gtgt c 1 gtgt for k 1100000 y(k) rand x(k)
-clog(y(k)) end hist(x,100)
histogram of
histogram of
17
Generating Random Numbers
Theorem If is uniformly distributed over 0,1
and is independent of and is
exponentially distributed with expected value 2
then and and independent
and normally distributed with mean 0 and variance
1
Proof The conditions on x and y are satisfied iff
18
Generating Random Numbers
The MATLAB commands generate 100000 samples
gtgt for k 150000 u(k) rand r(k)
-2log(rand) x(2k-1) sqrt(r(k))cos(2piu(k))
x(2k) sqrt(r(k))sin(2piu(k)) end hist(x,1
00)
19
Simulating Choosing m from n People
The MATLAB commands below do this
people are called 1, 2, ,n c vector
containing chosen people r vector of people
not yet chosen r 1n for k 1m compute
random integer between 1 and n-k1 j
round((n-k1)rand0.5) c (k) r(j) r(j)
r(n-k1) r r(1n-k) end
20
Suggested ReadingProblems in Textbook
5.1 Area Under a Curve, p 177-181, Prob 5.1
5.2 Generating Random Numbers, p 177-181, Prob 5.2
5.3 Simulating Prob. Behavior, p 186-190, Prob 5.3
Recommended Websites
http//en.wikipedia.org/wiki/Monte_Carlo_method
http//en.wikipedia.org/wiki/Random_number_generat
or
http//en.wikipedia.org/wiki/Queueing_theory
21
Tutorial 6 Due Week 6-10 October
Page 181. Problem 1. Imagine that every week you
purchase lottery tickets until you win a prize.
Design and run a simulation program to compute
the average number of lottery tickets that you
purchase per week.
Page 181. Problem 3. Choose a number of random
points so that your estimated value of pi is
accurate to about 4 significant digits
Page 190. Problem 1 (not project 1)
22
Homework 2 Due Friday 10 October
1. Write a computer program to simulate waiting
times at a bus station during rush hours using
the assumptions and steps below. Run the program
200 times to generate a histogram of the waiting
times that people experience and a histogram of
the total number of bus arrivals. Explain in
detail your logic and algorithms.
a. People start to arrive at the bus station at
500 at the average rate of 8 per minute. They
stop arriving at 700. The times between
consecutive arrivals is exponentially
distributed. For each simulation, first compute
an array containing random samples of the times
(in order) that people arrive between 5-7PM.
b. Buses arrive promptly every 10 minutes
starting at 510 and continue until the last
passenger is picked up. Each bus arrives empty
and picks up exactly 60 people.
c. Each time a bus arrives the people waiting
scramble to board the bus. Simulate this by
choosing 60 people randomly, from among those
waiting, during each bus arrival.