CSE 541 Numerical Methods

1
CSE 541Numerical Methods

• Monte Carlo Methods

2
Monte-Carlo Methods

• A Monte-Carlo method refers to any method that
  makes use of random numbers
• Numerical Analysis
• Simulation of natural phenomena
• Simulation of experimental apparatus

3
Monte-Carlo Techniques

• Sample Applications
• Integration
• System simulation
• Computer graphics modeling and rendering
• Physical phenomena - radiation transport
• Communications - bit error rates
• VLSI designs - tolerance analysis

4
Modelling Fractal Mountains

• Fractals represent controlled randomness to
  produce natural looking behavior

5
A Monte-Carlo Method to Approximate p

• Take a unit circle (centered at the origin) and a
  bounding box (square with side length 2)
• Acircle p
• Asquare 4
• Take a point (x, y) inside the square
• The probability ((x, y) -1 x 1 and -1 y
  1 is inside the circle) chance the point
  will be inside the circle

outside
inside
For example, what is the chance of a heads if you
flip a coin?
6
A Monte-Carlo Method to Approximate p

• OK, lets flip a coin
• Consider the following procedure
• Initialize Ncircle 0 and Nsquare 0
• Choose a random point (x, y)
• If (x, y) is inside the circle then increment
  Ncircle
• Increment Nsquare
• Repeat

Run this N times
7
A Monte-Carlo Method to Approximate p

• The probability a random point is inside the unit
  circle
• Choose N random points and let M Ncircle
• Approximate p

As N Þ a better approximation
8
Estimate p Using Monte-Carlo

• Results
• N 10,000 p 3.1388
• N 100,000 p 3.1452
• N 1,000,000 p 3.14164
• N 10,000,000 p 3.1422784

9
Choosing Random Numbers

• Does it matter how the N random numbers are
  chosen?
• Choose random numbers evenly across the space

p 4
p 4
p 0
10
Estimate p Using Monte-Carlo

• double x, y, pi
• const int m_nMaxSamples 100000000
• int count 0
• for (int k 0 k
• x 2.0 drand48() 1.0 // Map to the range
  -1,1
• y 2.0 drand48() 1.0
• if (xx yy
• pi 4.0 (double) count / (double)m_nMaxSamples

11
Random Numbers

• Definition of random from Merriam-Webster
• Main Entry randomFunction adjectiveDate
  15651 a lacking a definite plan, purpose, or
  pattern b made, done, or chosen at random random passages from the book2 a relating to,
  having, or being elements or events with definite
  probability of occurrence b
  being or relating to a set or to an element of a
  set each of whose elements has equal probability
  of occurrence also
  characterized by procedures designed to obtain
  such sets or elements

12
Random Numbers

• Random numbers
• A set of numbers that have nothing to do with the
  other numbers in the sequence
• In a uniform distribution of random numbers in
  the range 0,1 , every number has the same
  chance of turning up
• 0.00001 is just as likely as 0.5000

13
Random Sampling

• Use some chaotic system (Balls in a barrel
  Lotto)
• Use a process that is inherently random
• Radioactive decay
• Thermal noise
• Cosmic ray arrival
• Tables of a few million random numbers
• Hooking up a random machine to a computer

14
Computer Generated Random Numbers

• Definition of algorithm from dictionary.com
• algorithm
• n a precise rule (or set of rules) specifying
  how to solve some problem.

15
Random vs. Pseudo-random

• Random numbers have no defined sequence or
  formulation. Thus, for any n random numbers, each
  appears with equal probability.
• Computer limitations
• Algorithms A precise rule (or set of rules)
  specifying how to solve some problem. Implies one
  value depends on the others.
• Machine precision If we restrict ourselves to
  the set of 32-bit integers, then our numbers will
  start to repeat after some very large n. The
  numbers thus clump within this range and around
  these integers.
• pseudo-random adj. Of, relating to, or being
  random numbers generated by a definite, nonrandom
  computational process.
• Computer algorithms are restricted to generating
  pseudo-random numbers

16
John Von Neumann (circa 1946)

• An algorithm to randomly generate 10 digit
  integers
• Start with a 10 digit integer
• Square it and take middle 10 digits from the
  answer
• Example 57721566492 33317792380594909291
• Need to map these to a particular range, e.g. 0,
  1
• Smaller numbers lump together and create short
  loops
• For example, 4 digit integers
• 61002 37210000
• 21002 04410000
• 41002 16810000
• 51002 65610000

17
Pseudo-random Number Generators

• Most pseudo-random number generators use the
  following two tools
• Modulo arithmetic
• Large prime numbers

18
Linear Congruential Methods

• Lehmer (circa 1948)
• where a, c ³ 0 , m I0,a,c
• Advantage Very fast
• Problems
• Poor choice of the constants can lead to very
  poor sequences
• The relationship will repeat with a period no
  greater than m

19
The RANDU Generator (IBM, circa 1960s)

• Algorithm
• This generator was later found to have a serious
  problem

2D Distribution
20
Prime Numbers

• The remainder when a number is divided by a prime
  is random
• The number 231 1 is a prime number, a 75, c
  0, m 231 1
• This algorithm calculates some number and then
  takes the lower-order bits (mod operator)
• The integers in the sequence are mapped to 0,
  1, the xis
• I0 is called the seed of the process

21
The Algorithm Drand48

• Let a 2736731631558, c 138, and m 248
• Uses higher-bit integers to generate 48-bit
  random numbers

22
Other Ideas

• Why use integers?
• Was more efficient on old computers
• Why use one seed?
• Has a greater period than m
• The RANMAR generator (CERN Library) uses
• 103 initial seeds and a period of 1043
• The ultimate random number generator?

23
Properties of Pseudo-random Numbers

• These algorithms generate periodic sequences
  (hence not random). A random number generated may
  match our initial seed.
• The restriction to quantized numbers (a
  finite-set), leads to problems in
  high-dimensional space
• The individual digits in the random number may
  not be independent. There may be a higher
  probability that a 3 will follow a 5.

24
Property 2 Higher Dimensional Space

• Quantized numbers lead to problems in
  high-dimensional space
• Many points end up to be co-planar
• For example, in two dimensions numbers may
  cluster to lines

25
Property 2 Higher Dimensional Space

• The Marsaglia Effect (1968)
• Random numbers fall mainly in planes (or
  hyperplanes)
• 3D distribution of RANDU

Problems seen when observed at the right angle
26
The C Standard Library Functions

• include
• Function rand (type man rand on the Unix
  command line)
• Only results in 16-bit integers
• Has a periodicity of 231 though
• Rather poor pseudo-random number generator
• Function random (type man random)
• Results in 32-bit integers
• Uses rand() to build an initial table
• Has a periodicity of around 269
• Slightly better pseudo-random number generator
• Function drand48 (type man drand48)
• Uses double precision
• Pretty good routine
• May not be as portable

27
Initializing With Seeds

• You can set the seed with srand, srandom or
  srand48
• A seed value will generate the same sequence of
  random numbers
• Deterministic program using the same seed will
  produce the same program output on each run
• Non-deterministic program using a different seed
  on each run will produce different program
  outputs
• For example, call the seed method with the
  lower-order bits of the system clock

28
Deterministic Program

• A procedurally generated mountain
• If we re-generate the fractal mountain every time
  the camera moves, we dont want the mountain to
  change

29
Mapping Random Numbers

• Most computer library support for random numbers
  only provides random numbers over a fixed range
• Random integers from zero to some maximum M
• Random floating-point or double-precision numbers
  mapped to the range 0 1
• You need to map this to your desired range
• Our approximation to p example requires numbers
  in the range -1 1
• Use translate and scale operations (remember
  Gaussian Quadrature)

30
Monte-Carlo Techniques

• Problem What is the probability that 10 dice
  throws add up exactly to 32?
• Exact Way. Calculate this exactly by counting all
  possible ways of making 32 from 10 dice
• We may not always be able to compute an exact
  way
• Approximate (Lazy) Way. Simulate throwing the
  dice (say 500 times), count the number of times
  the results add up to 32, and divide this by 500.
  
• Lazy Way can get quite close to the correct
  answer quite quickly

31
1D Function Integration

• Exact way Analytical Integration
• Approximate way Quadrature
• Approximate way Monte-Carlo
• Randomly sample the area enclosed by x Î a, b
  and y Î 0, max (p(x))

32
1D Function Integration

• Intuitively

33
Monte-Carlo And Integration

• A Monte-Carlo approach can be advantageous for
  integration if
• The function is defined over a complicated domain
• Integration in higher dimensions

34
Complicated Domain

• Define two domains
• D Complicated domain
• D Simple domain and a superset of D
• Pick random points over D
• Count
• N points over D
• N points over D

35
Complex and Higher Dimensions

• Higher dimensional shapes can be complex
• Need a very fine grid for standard quadrature to
  get samples in Region R
• Monte-Carlo methods can make use of all points in
  the simpler domain

36
Higher Dimensions

• Our Standard Quadrature approaches compute a
  numerical value of a definite integral by
• where points xi are uniformly spaced

37
Higher Dimensions

• Consider integral in d dimensions
• The error with N sampling points is

38
Monte-Carlo Error

• From probability theory, the Monte Carlo error
  decreases with sample size N as
• independent of dimension d