Computational Methods in Physics PHYS 3437

1
Computational Methods in Physics PHYS 3437

• Dr Rob Thacker
• Dept of Astronomy Physics (MM-301C)
• thacker_at_ap.smu.ca

2
Todays Lecture

• Introduction to Monte Carlo methods
• Background
• Integration techniques

3
Introduction

• Monte Carlo refers to the use of random numbers
  to model random events that may model a
  mathematical of physical problem
• Typically, MC methods require many millions of
  random numbers
• Of course, computers cannot actually generate
  truly random numbers
• However, we can make the period of repetition
  absolutely enormous
• Such pseudo-random number generators are based on
  truncation of numbers of their significant digits
• See Numerical Recipes, p 266-280 (2nd edition
  FORTRAN)

4

• Anyone who considers arithmetical methods of
  producing random digits is, of course, in a state
  of sin.

John von Neumann
5
History of numerical Monte Carlo methods

• Another contribution to numerical methods related
  to research at Los Alamos
• Late 1940s scientists want to follow paths of
  neutrons following various sub-atomic collision
  events
• Ulam von Neumann suggest using random sampling
  to estimate this process
• 100 events can be calculated in 5 hours on ENIAC
• The method is given the name Monte Carlo by
  Nicholas Metropolis
• Explosion of inappropriate use in the 1950s gave
  the technique a bad name
• Subsequent research illuminated when the method
  was appropriate

6
Terminology

• Random deviate a distribution of numbers
  choosen uniformly between 0,1
• Normal deviate numbers chosen randomly between
  (-8,8) weighted by a Gaussian

7
Background to MC integration

• Suppose we have a definite integral
• Given a good set of N sample points xi we can
  estimate the integral as

Sample points e.g. x3 x9
Each sample point yields an element of the
integral of width (b-a)/N and height f(xi)
f(x)
a
b
8
What MC integration really does

• While the previous explanation is a reasonable
  interpretation of the way MC integration works,
  the most popular explanation
  is below

Height given by random sample of f(x)
Average
a
b
9
Mathematical Applications

• Lets formalize this just a little bit
• Since by the mean value theorem
• We can approximate the integral by calculating
  (b-a)ltfgt, and we can calculate ltfgt by averaging
  many values of f(x)
• Where xi?a,b and the values are chosen randomly

10
Example

• Consider evaluating
• Lets take N1000, then evaluate f(x)ex with
  x?0,1 at 1000 random points
• For this set of points define
• I1(b-a)ltfgtN,1ltfgtN,1 since b-a1
• Next choose 1000 different x?0,1 and create a
  new estimate I2ltfgtN,2
• Next choose another 1000 different x?0,1 and
  create a new estimate I3ltfgtN,3

11
Distribution of the estimates

• We can carry on doing this, say 10,000 times at
  which point well have 10,000 values estimating
  the integral, and the distribution of these
  values will be a normal distribution
• The distribution of the all of the IN integrals
  constrains the errors we would expect on a single
  IN estimate
• This is the Central Limit Theorem, for any given
  IN estimate the sum of the random variables
  within it will converge toward a normal
  distribution
• Specifically, the standard deviation will be the
  estimate of the error in a single IN estimate
• The mean, x0, will approach e-1

1
1/e
x0-sN
x0sN
x0
12
Calculating sN

• The formula for the standard deviation of N
  samples is
• If there is no deviation in the data then RHS is
  zero
• Given some deviation as N?8, the RHS will settle
  to some constant value gt 0 (in this case
  0.2420359)
• Thus we can write

A rough measure of how good a random number
generator is how well does a histogram of the
10,000 estimates fit to a Gaussian.
13
Add mc.ps plot
Increasing the number of integral estimates makes
the distribution closer and closer to the
infinite limit.
1000 samples per I integration Standard
deviation is 0.491/v1000
14
Resulting statistics

• For data that fits a Gaussian, the theory of
  probability distribution functions asserts that
• 68.3 of the data (ltfgtN) will fall within sN of
  the mean
• 95.4 of the data (19/20) will fall within 2sN
  of the mean
• 99.7 of the data will fall within 3sN etc
• Interpretation of poll data
• These results will be accurate to 4, (19 times
  out of 20)
• The 4 corresponds to 2s ? s0.02
• Since s1/sqrt(N) ? N2500
• This highlights one of the difficulties with
  random sampling, to improve the result by a
  factor of 2 we must increase N by a factor of 4!

15
Why would we use this method to evaluate
integrals?

• For 1D it doesnt make a lot of sense
• Taking h1/N then composite trapezoid rule
  errorh21/N2N-2
• Double N, get result 4 times better
• In 2D, we can use an extension of the trapezoid
  rule to use squares
• Taking h1/N1/2 then error?h2?N-1
• In 3D we get h1/N1/3 then error?h2?N-2/3
• In 4D we get h1/N1/4 then error?h2?N-1/2

16
MC beneficial for Ngt4

• Monte Carlo methods always have sNN-1/2
  regardless of the dimension
• Comparing to the 4D convergence behaviour we see
  that MC integration becomes practical at this
  point
• It wouldnt make any sense for 3D though
• For anything higher than 4D (e.g. 6D,9D which are
  possible!) MC methods tend to be the only way of
  doing these calculations
• MC methods also have the useful property of being
  comparatively immune to singularities, provided
  that
• The random generator doesnt hit the singularity
• The integral does indeed exist!

17
Importance sampling

• In reality many integrals have functions that
  vary rapidly in one part of the number line and
  more slowly in others
• To capture this behaviour with MC methods
  requires that we introduce some way of putting
  our points where we need them the most
• We really want to introduce a new function into
  the problem, one that allows us to put the
  samples in the right places

18
General outline

• Suppose we have two similar functions g(x)
  f(x), and g(x) is easy to integrate, then

19
General outline cont

• The integral we have derived
  has some nice properties
• Because g(x)f(x) (i.e. g(x) is a reasonable
  approximation of f(x) that is easy to integrate)
  then the integrand should be approximately 1
• and the integrand shouldnt vary much!
• It should be possible to calculate a good
  approximation with a fairly small number of
  samples
• Thus by applying the change of variables and
  mapping our sample points we get a better answer
  with fewer samples

20
Example

• Lets look at integrating f(x)ex again on 0,1
• MC random samples are 0.23,0.69,0.51,0.93
• Our integral estimate is then

21
(No Transcript)
22
Apply importance sampling

• We first need to decide on our g(x) function, as
  a guess let us take g(x)1x
• Well it isnt really a guess we know this is
  the first two terms of the Taylor expansion of
  ex!
• y(x) is thus given by
• For end points we get y(0)0, y(1)3/2
• Rearrange y(x) to give x(y)

23
Set up integral evaluate samples

• The integral to evaluate is now
• We must now choose ys on the interval 0,3/2

y
0.345 1.038
1.035 1.211
0.765 1.135
1.395 1.324
Close to 1 because g(x)f(x)
24
Evaluate

• For the new integral we have
• Clearly this technique of ensuring the integrand
  doesnt vary too much is extremely powerful
• Importance sampling is particularly important in
  multidimensional integrals and can add 1 or 2
  significant figures of accuracy for a minimal
  amount of effort

25
Rejection technique

• Thus far weve look in detail at the effect of
  changing sample points on the overall estimate of
  the integral
• An alternative approach may be necessary when you
  cannot easily sample the desired region which
  well call W
• Particularly important in multi-dimensional
  integrals when you can calculate the integral for
  a simple boundary but not a complex one
• We define a larger region V that includes W
• Note you must also be able to calculate the size
  of V easily
• The sample function is then redefined to be zero
  outside the volume, but have its normal value
  within it

26
Rejection technique diagram
f(x)
V
W
Region we want to calculate
Area of Wintegral of region V multiplied
by fraction of points falling below f(x) within V
Algorithm just count the total number of points
calculated the number in W!
27
Better selection of points sub-random sequences

• Choosing N points using a uniform deviate
  produces an error that converges as N-0.5
• If we could choose points better we could make
  convergence faster
• For example, using a Cartesian grid of points
  leads to a method that converges as N-1
• Think of Cartesian points as avoiding one
  another and thus sampling a given region more
  completely
• However, we dont know a priori how fine the grid
  should be
• We want to avoid short range correlations
  particles shouldnt be too close to one another
• A better solution is to choose points that
  attempt to maximally avoid one another

28
A list of sub-random sequences

• Tore-SQRT sequences
• Van der Corput Halton sequences
• Faure sequence
• Generalized Faure sequence
• Nets (t,s)-sequences
• Sobol sequence
• Niederreiter sequence
• Well look very briefly at Halton Sobol
  sequences, both of which are covered in detail in
  Numerical Recipes

Many to choose from!
29
Haltons sequence

• Suppose in 1d we obtain the jth number in
  sequence, denoted Hj, via
• (1) write j as a number in base b, where b is
  prime
• e.g. 17 in base 3 is 122
• (2) Reverse the digits and place a radix point in
  front
• e.g. 0.221 base 3
• It should be clear why this works, adding an
  additional digit makes the mesh of numbers
  progressively finer
• For a sequence of points in n dimensions
  (xi1,,xin) we would typically use the first n
  primes to generate separate sequences for each of
  the xij components

30
2d Haltons sequence example
Pairs of points constructed from base 3 5
Halton sequences
31
Sobol (1967) sequence

• Useful method described in Numerical Recipes as
  providing close to N-1 convergence rate
• Algorithms are also available at www.netlib.org

From Num. Recipes
32
Summary

• MC methods are a useful way of numerically
  integrating systems that are not tractable by
  other methods
• The key part of MC methods is the N-0.5
  convergence rate
• Numerical integration techniques can be greatly
  improved using importance sampling
• If you cannot write down a function easily then
  the rejection technique can often be employed

33
Next Lecture

• More on MC methods simulating random walks