Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
Week 4: aim to cover
Monte Carlo integration, floating point numbers, roundo↵ error (Lecture 7)
Confidence intervals, MC integration, roundo↵ error (Lab 4) Error propagation (Lecture 8)
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
Monte Carlo integration
It is frequently necessary to compute definite integrals. The functions being integrated are often quite complicated and it may not be possible to find an indefinite integral in closed form. Since what is wanted is a numerical result, computers are used to find numerical approximations to the definite integral. In the simplest case, we are given a function f (x ) and two limits of integration a and b and the object is to approximate
Zb
f (x)dx.
a
There are methods that are completely deterministic — every time you run the program you will get exactly the same answer.
There is another class of methods that use random numbers to compute the value of definite integrals — they are called Monte Carlo methods, after the famous Monaco casino. They are uncompetitive for simple one-dimensional integrals, but prove to be superior for complicated integrals over many variables.
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
Hit-and-Miss method
One method — so-called hit-and-miss Monte Carlo — uses the relation of the integral to the area under the graph of f (x). Suppose we want the
value of Z b
f(x)dx
a
where f (x) 0 and we know the maximum value of f (x) over [a,b]: 0 f (x) M
Then the rectangle just including all the graph of f (x ) over [a,b] has area (b a)M. In other words, the fraction of the rectangle lying under the
c u r v e y = f ( x ) i s j u s t R ab f ( x ) d x
(b a)M
In the hit-and-miss method, we generate points uniformly in the rectangle (generate an x-coordinate from U(a,b), generate a y-coordinate from U(0,M), then form the point (x,y)) and count the fraction of them that lie under the curve y = f (x).
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
Bu↵on’s needle (1733)
If we randomly toss a needle of length l = 1 onto a floor with floorboards of width d = 2, what is the probability the needle crosses over a crack between floorboards?
The needle crosses a crack if y < l sin ✓ where y is the distance of the 2
centre of the needle to the nearest crack, and ✓ is the angle of the upper part of the needle relative to the cracks, measured from the positive direction.Since0
a
In the simple Monte Carlo method, we generate values of x from U(a,b) and calculate the sample mean of the set of values
f ̄ ⌘ 1 X f ( x i ) ni
Our estimate of tZhe integral is then just b a times the sample mean.
b
f(x)dx =(b a)
Both MC methods extend easily to higher-dimensional integrals.
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
A comparison
Here is a comparison for a simple 1-dimensional integral, so not really typical of the applications that MC integration is used for. We
approximate
Z1
exp( x) dx = 1 exp( 1) ⇡ 0.63212
I =
First using both MC methods and a method using quasirandom numbers.
0
10 -1
10 -2
10 -3
10 -4
10 -5
102 103
104 105
Errors for MC and QR integration
n-1/2
log(n)/n
HitMiss
MC
QR
n
error
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
A comparison
Then we compare simple MC integration with 2 classical methods — the trapezoid rule and Simpson’s rule.
100
10 -2
10 -4
10 -6
10 -8 10 -10
10 -12 10 -14 10 -16 10 -18
10 -20
102 103
104 105
MC and classical quadrature
MC
trap
simp
n-1/2
n-2
n-4
n
error
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
In higher dimensions
The various integration methods di↵er in how they extend to higher dimensions. The classical methods require many points in each dimension, which makes them very expensive as the dimension d grows. Stated another way, the error falls very slowly with the number of evaluation points once the dimension gets larger than, say, 10.
Method MC
QR trapezoid Simpson’s rule
Error
n 1/2 log(n)d n 1
n 2/d n 4/d
Table: Asymptotic error bound for large n, d
We see that MC methods have an error bound independent of the
dimension, and win out for d > 8 or so.
This is why for many complicated processes that are equivalent to high-dimensional integrals (like pricing exotic stock options), Monte Carlo methods are the method of choice.
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
Overview of MC integration
Advantages:
error decay rate independent of dimension can estimate the statistical error
easy to program and parallelize
Disadvantages:
error decays slowly for low dimensions
sometimes the variance is very large, so need special methods to reduce it
sometimes the estimates are correlated, so more care required to estimate error
Numerical Methods & Scientific Computing: lecture notes
Stochastic simulation
Statistical errors
More advanced topics: not in MAST30028
Variance reduction methods e.g. importance sampling Markov chain methods
Monte Carlo methods for optimization
Numerical Methods & Scientific Computing: lecture notes
Errors
How numbers are stored in a computer
To understand one ubiquitous source of error in scientific computing (roundo↵ error), have to understand how numbers are stored.
In scientific problems, numbers can vary greatly in magnitude =) we try to store numbers with a fixed relative error (precision), not absolute error. Floating Point numbers were invented to do this.
We describe the current standard: IEEE 754 used by most computers (perhaps not completely) for double precision numbers
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Floating point numbers
To represent numbers with greatly varying size, we use scientific notation
In the computer, we do the same but store both the exponent and fraction in base 2 (binary representation).
10112 representstheinteger23+21+20 =8+0+2+1=1110
=) we approximate real numbers by a nearby rational number with a finite binary expansion
Example
6.023 ⇥ 1023, 1.055 ⇥ 10 34
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Floating point numbers
x=±(1.f1f2···ft)⇥2e =±(1+f)2e
f isthefractionormantissa: 0f <1,eachbinarydigitfi =0or1
e is the exponent, a binary integer
Example
125.75=64+32+16+8+4+1+0.5+0.25=1111101.112 = 1.1111101112 ⇥ 26 = 1.1111101112 ⇥ 21102
so 1 + f = 1.1111101112, e = 1102
e determines the range of numbers that can be represented
f determines the number of significant digits carried ! the precision of numbers that can be represented
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Double precision numbers
In double precision, use 64 bits (8 bytes) to store each number 1 for sign ± 11
11 for exponent ! 2 = 2048 possible values for e
52 for fraction =) t = 52
=) to store 106 double precision numbers takes 8 Mbytes.
To get the exponent of the number, subtract 1023 from the 11-bit integer 2 [1, 2046]! e 2 [ 1022, 1023] (so we get equal numbers of ‘small’ and ‘large’ numbers).
The extreme values e = 1023, 1024 are used for special cases or error flags.
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Overflow
Using the largest exponent e = 1023 and the largest 1 + f = 1.11 · · · 12 we get the largest number representable ⇡ 21024 ⇡ 10308
called realmax in MATLAB
Any number larger than this causes Overflow
and is stored as infinity, with e = 1024, f = 0 called Inf in MATLAB
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Underflow
Using the smallest exponent e = 1022 and the smallest
1 + f = 1.00 · · · 02 we get the smallest (non-zero, normalized) number representable realmin= 2 1022 ⇡ 10 308
If e = 1023, use denormalized numbers (x = ±f 2 1022)
=) smallest non-zero number representable = 2 1074 ⇡ 5 ⇥ 10 324 Any number smaller than this gives Underflow and ! 0
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Exceptions
Finally
Results of illegal operations (divide by zero, etc. ) use e = 1024, f 6= 0
and are stored as NaN (Not a Number) in MATLAB.
Both NaN and Inf propagate in subsequent calculations, once they have been produced.
So you can diagnose the problem!
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Machine numbers
Only numbers 2 [realmin, realmax] with a binary fraction expansion that terminates in less than 53 digits + the denormalized numbers can be represented exactly as a Floating Point number ! the set of machine numbers , denoted by F
125.75 = 1.1111101112 ⇥ 21102 which needs only 9 binary digits in the fraction =) 125.75 is a machine number
all the rest will be approximated with an error ! Roundo↵ error
If rounding errors vanished, 95% of numerical analysis would
remain.
Trefethen’s Maxims
Example
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
The structure of the Floating Point number system
Since there are 252 numbers in [1,2), 252 numbers in [2,4), 252 numbers in [4,8) etc.
The machine numbers are nonuniformly distributed
The granularity of the machine numbers is specified by the gap between 1 and the next machine number 1 + 2 t , called machine epsilon "M Hence IEEE double precision numbers have machine epsilon
2 52 ⇡ 2 ⇥ 10 16, given by eps in MATLAB.
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Roundo↵ error
If you try to produce a non-machine number x, must store it as the nearest machine number produced by rounding — we denote this number fl(x).
The default rounding mode: round to nearest, then even
=) maximumerror="M/2⇥2e =u⇥2e where unit roundo↵ u = "M /2.
=) max. relative error
=u⇥2e/(1+f)⇥2e =u/(1+f)2(u/2,u]u
unit roundo↵ determines the precision with which floating point numbers are stored
So
u fl(x) x u x
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Floating point precision
Since t = 52, IEEE double precision numbers have precision 2 53 ⇡ 10 16
16 decimal digits of precision
Can also have IEEE single precision numbers with unit roundo↵ 2 24 ⇡ 10 7
7 decimal digits of precision e.g. ANSI C float, MATLAB single
We make this error just in storing the numbers in a computer.
Floating point precision sets a lower bound to the level of (relative) data error
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Floating point arithmetic
From the error on storing a number
u fl(x) x u
we get
x fl(x)=x(1+ ), | |u 8x2R
i.e. the nearest machine number fl(x) to x is at most a factor 1±u away. Since arithmetic operations on machine numbers produce results that usually are not machine numbers so must get rounded, we get a model for Floating point arithmetic (ignoring underflow/overflow)
fl(xopy)=(xopy)(1+ ), | |u8x,y2F,op=+, ,⇥,÷
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
Summary
1 Some numbers (machine numbers) can be represented exactly as a Floating Point number; most can’t
2 there’s a smallest and largest Floating Point number
3 Floating Point numbers have a precision u — the source of inevitable roundo↵ error
Numerical Methods & Scientific Computing: lecture notes
Errors
Floating point numbers
End of Lecture 7