Transcript SimSET: a Simulation System for Emission Tomography

Introduction to Monte Carlo
Simulation
Robert L. Harrison
University of Washington Medical Center
Seattle, Washington, USA
Supported in part by PHS grants CA42593 and CA126593
What is Monte Carlo?
• Monte Carlo
methods are a class
of algorithms that
rely on repeated
random sampling to
compute their
results. (Wikipedia)
An example: estimating an
integral
1
• Randomly generate
points in square.
100
400
1000points.
points.
•• Count
the points
• 278
69
under
curve.
curve.
717under
under
the curve.
Estimate
0.69.
0.695.
0.717.
•• Divide
the=number
under curve by total
number of points.
0
0
0
0
1
1
A little history…
• Buffon’s needle.
• Student / William Gosset.
• Manhattan project.
Buffon’s needle
• Georges-Louis Leclerc,
Comte de Buffon (1707 –
1788).
• A party game using a
lined piece of paper and
a pin.
• Estimates the value of 
 = 2*drops/hits.
http://www.shodor.org/interactivate/activities/buffon
William Gosset (1876-1937)
• Supervised the experimental brewery at
Guiness.
• Developed methods to draw
conclusions from small sample size
experiments.
• Published under the penname A.
Student.
• Student t-test.
Manhattan project
• Development of first nuclear weapons.
• Needed to extrapolate knowledge of a
single neutron’s behavior to the
behavior of a system of neutrons.
• Experimentation was to dangerous,
time-consuming, costly.
Manhattan project…
• John von Neumann (1903-1957) and
Stanislaus Ulam (1909-1984).
• Developed particle tracking Monte Carlo
simulation.
Why do people simulate?
Why do people simulate?
• Low cost way to experiment.
• Easy to repeat/modify.
• Study systems to complex to solve
analytically.
• Examine hidden quantities.
Disadvantages
•
•
•
•
Time consuming.
Doesn’t give exact solution.
Is the model correct?
Bugs, bugs, bugs.
Alternatives
• Analytic solution.
• Experiment.
Mechanics of simulation
•
•
•
•
General method.
Random number generator.
Sampling from a distribution.
Accelerating a simulation.
General Monte Carlo
Method
• Model a physical system as a (series of)
probability density function(s).
• Repeatedly sample from the pdf(s).
• Tally the statistics of interest.
• (There are many Monte Carlo
simulations that don’t fit the above
exactly.)
Random numbers
• Generally we sample from a uniform
distribution of integers over [0,n].
• Convert to a uniform distribution on the
interval (0,0), [0,1) or [0,1].
– (e.g., divide integer value by n or n+1.)
• Lists of truly random numbers.
– E.g., time between radioactive decays.
• Pseudo-random number generators.
– Calculate a randoms number!?
Pseudo-random number
generators
• Usually called random number
generators (RNG).
• (Early - bad) Middle squares.
• (Still around - bad) Linear congruential
generator.
• (Good) Mersenne twister.
Middle square
• Von Neumann, 1946.
To generate a sequence of 10 digit integers,
start with one, and square it and then take the
middle 10 digits from the result as the next
number in the sequence
2
• e.g. 5772156649 = 33317792380594909291
the next number is
• Not a recommended method!
Linear congruential
generator
• Lehmer, 1948.
In+1 = (aIn + c) mod m
• I0 = Starting value (seed)
• a, c, and m are specially chosen
constants.
Linear congruential
generator…
• Values used by
some built-in
random functions
are very poor.
• Good values for a,
c, and m on-line.
• Even with good
values, there can be
problems.
Three dimensional sampling from
a linear congruential generator
with poorly chosen values,
RANDU (Wikipedia).
Mersenne twister
• Used in MatLab, Maple.
• Does very well on most tests of
randomness.
• Ridiculously long period.
• Reasonably fast.
Mersenne twister
• Used in MatLab, Maple.
• Does very well on most tests of
randomness.
• Ridiculously long period.
• Reasonably fast.
• Quite complicated to
implement/describe.
Sampling from a distribution
• We have a uniform distribution, how do
we sample from
– An exponential distribution?
– Other continuous distributions?
– A Poisson distribution?
Exponential distribution
• Time between
radioactive decays.
• Distance a photon
travels before an
interaction.
• Mean 1/, variance
1/
Probability density function (pdf)
(Wikipedia)
Cumulative distribution function (cdf)
Sampling from the
exponential distribution
• Sample u from the
uniform distribution
on [0,1).
Sampling from the
exponential distribution
• Sample u from [0,1).
• Locate u on the yaxis of the
u
exponential
distribution cdf.
Sampling from the
exponential distribution
• Sample u from [0,1).
• Locate u on y-axis.
• Find the x value that
u
corresponds to u.
x
Sampling from the
exponential distribution
• Sample u from [0,1).
• Locate u on y-axis.
• Find x = F-1(u).
u
x
Sampling from an invertible
continuous distribution
• Sample u from the unit interval.
• Find x = F-1(u).
Other sampling methods for
continuous distributions
• Approximate F-1(u), e.g., as a piecewise
linear function.
– Table lookups are fast!
– Approximate => bias.
• Acceptance-rejection method.
– AKA rejection method, accept-reject
method, rejection sampling.
Acceptance-rejection
•
To sample from a
pdf f(x):
1. Choose pdf g(x)
and constant c such
that
c*g(x) ≥ f(x) for all x
www.mathworks.com
Acceptance-rejection
•
To sample from f(x):
1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate a random
number v from g(x).
3. Generate a random
number u from the
uniform distribution
on (0,1).
v
Acceptance-rejection
•
To sample from f(x):
1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate v from g(x).
3. Generate u from
uniform distribution.
4. If c*u ≤ f(v)/g(v)
accept v
5. else reject v and go
to 2.
Acceptance-rejection
•
To sample from f(x):
1. Choose g(x) and c:
c*g(x) ≥ f(x) for all x.
2. Generate v from g(x).
3. Generate uniform u.
4. If c*u ≤ f(v)/g(v)
accept v
5. else reject, go to 2.
•
The average number
of iterations is c.
Acceptance-rejection
• Sampling from a normal distribution.
– While the cdf is invertible, the resulting formula is
computationally inefficient.
– Generally sampling is done using the Box-Muller
or Ziggurat algorithms.
– Use a library function!
• Acceptance-rejection also be used to remove
the bias from an approximation/table lookup
sampling method.
Sampling from a normal
distribution
• Gaussian, bell curve.
• Everything is normal.
– Central limit theorem:
the sum of a large
number of independent
random variables is
approximately normally
distributed.
pdf
• Generally sampled
using the acceptancerejection method.
• Use a library function!
cdf
(Wikipedia)
Poisson distribution
• Discrete probability
distribution.
• {0,1,2…} (nonnegative integer)
result.
Probability mass function
(Wikipedia)
• Mean = variance = 
Cumulative density function
Poisson distribution
• Distribution for: how
many light bulbs will
burn out in a month?
– Large building.
– Maintenance
department immediately
replaces burned out
bulb.
– Choose  = the average
number that burn out in
a month
How many
mathematicians
does it take to
change a light bulb?
(www.clipartof.com)
In earlier work,
Haynor[1] has
shown that one
mathematician can
change a light bulb.
Poisson distribution
• Other examples:
– Number of radioactive
decays.
– Number of bugs
discovered in a
program.
– Number of calls to a
call center.
– Number of scintillation
photons.
How many
mathematicians
does it take to
change a light bulb?
None. It's left to the
reader as an
(www.clipartof.com)
exercise.
Sampling from Poisson
distributions
•
•
Use a library function!
Use the exponential
distribution.
1. Let T be the time period,
set t = 0, k = 0.
2. Sample tk from the
exponential distribution.
3. t = t + tk, k = k + 1;
4. If t < T, go to 2.
5. Return (k - 1).
•
For  large approximate
using a normal
distribution.
How many quantum
physicists does it
take to change a
light bulb?
(www.clipartof.com)
They can't. If they
know where the
socket is, they
cannot locate the
new bulb.
Keeping track of results
• Tally (score, bin, histogram) the
quantities we are interested in, e.g.:
– Number of light bulbs replaced by month;
– Time to failure for a product.
Building a simulation
• Often a single ‘event’ in a simulation will
require sampling from many different
pdfs, some many times, and be tallied in
many ways.
Building a simulation
• For example, a nuclear physics simulation
might model
–
–
–
–
the time between decays
the decay process
the paths of photons and particles
the energy deposited during tracking and any
secondary particles created thereby
– the detection electronics…
• and tally
– energy, position, and type of detected particles.
Take away
• Simulation a key tool for analytically
intractable problems.
• Currently the Mersenne Twister is a good
choice of random number generator.
• Random sampling for functions:
–
–
–
–
cdfs are analytically invertible;
use library functions if available;
acceptance-rejection method;
piecewise approximation.
• Four key distributions: uniform, exponential,
normal, Poisson.
References
J.T. Bushberg, The essential physics of medical imaging, Lippincott Williams & Wilkins, 2002.
K.P. George et al, Brain Imaging in Neurocommunicative Disorders, in Medical speech-language pathology: a
practitioner's guide, ed. A.F. Johnson, Thieme, 1998.
D.E. Heron et al, FDG-PET and PET/CT in Radiation Therapy Simulation and Management of Patients Who Have
Primary and Recurrent Breast Cancer, PET Clin, 1:39–49, 2006.
E.G.A. Aird and J. Conway, CT simulation for radiotherapy treatment planning, British J Radiology, 75:937-949, 2002.
R. McGarry and A.T. Turrisi, Lung Cancer, in Handbook of Radiation Oncology: Basic Principles and Clinical
Protocols, ed. B.G. Haffty and L.D. Wilson, Jones & Bartlett Publishers, 2008.
R. Schmitz et al, The Physics of PET/CT Scanners, in PET and PET/CT: a clinical guide, ed. E. Lin and A. Alavi,
Thieme, 2005.
W.P. Segars and B.M.W. Tsui, Study of the efficacy of respiratory gating in myocardial SPECT using the new 4-D
NCAT phantom, IEEE Transactions on Nuclear Science, 49(3):675-679, 2002.
V.W. Stieber et al, Central Nervous System Tumors, in Technical Basis of Radiation Therapy: Practical Clinical
Applications, ed. S.H. Levitt et al, Springer, 2008.
P. Suetens, Fundamentals of medical imaging, Cambridge University Press, 2002.
www.wofford.org/ecs/ScientificProgramming/MonteCarlo/index.htm
www.shodor.org/interactivate/activities/buffon
www.impactscan.org/slides/impactcourse/introduction_to_ct_in_radiotherapy