9.8 Monte Carlo Simulation
This section under major construction.
In 1953 Enrico Fermi, John Pasta, and Stanslaw Ulam created the first "computer experiment" to study a vibrarting atomic lattice. Nonlinear system couldn't be analyzed by classical mathematics.
Simulation = analytic method that imitates a physical system. Monte Carlo simulation = use randomly generated values for uncertain variables. Named after famous casino in Monaco. At essentially each step in the evolution of the calculation, Repeat several times to generate range of possible scenarios, and average results. Widely applicable brute force solution. Computationally intensive, so use when other techniques fail. Typically, accuracy is proportional to square root of number of repetitions. Such techniques are widely applied in various domains including: designing nuclear reactors, predicting the evolution of stars, forecasting the stock market, etc.
Generating random numbers. The math library function Math.random generate a pseudo-random number greater than or equal to 0.0 and less than 1.0. If you want to generate random integers or booleans, the best way is to use the library Random. Program RandomDemo.java illustrates how to use it.
Note that you should only create a new Random object once per program. You will not get more "random" results by creating more than one. For debugging, you may wish to produce the same sequence of pseudo-random number each time your program executes. To do this, invoke the constructor with a long argument.
The pseudo-random number generator will use 1234567 as the seed. Use SecureRandom for cryptographically secure pseudo-random numbers, e.g., for cryptography or slot machines.
Linear congruential random number generator. With integer types we must be cognizant of overflow. Consider a * b (mod m) as an example (either in context of a^b (mod m) or linear congruential random number generator: Given constants a, c, m, and a seed x[0], iterate: x = (a * x + c) mod m. Park and Miller suggest a = 16807, m = 2147483647, c = 0 for 32-bit signed integers. To avoid overflow, use Schrage's method.
Exercise: compute cycle length.
Library of probability functions. OR-Objects contains many classic probability distributions and random number generators, including Normal, F, Chi Square, Gamma, Binomial, Poisson. You can download the jar file here. Program ProbDemo.java illustrates how to use it. It generate one random value from the gamma distribution and 5 from the binomial distribution. Note that the method is called getRandomScaler and not getRandomScalar.
Queuing models. M/M/1, etc. A manufacturing facility has M identical machines. Each machine fails after a time that is exponentially distributed with mean 1 / μ. A single repair person is responsible for maintaining all the machines, and the time to fix a machine is exponentially distributed with mean 1 / λ. Simulate the fraction of time in which no machines are operational.
Diffusion-limited aggregation.
Diffuse = undergo random walk. The physical process diffusion-limited aggregation (DLA) models the formation of an aggregate on a surface, including lichen growth, the generation of polymers out of solutions, carbon deposits on the walls of a cylinder of a Diesel engine, path of electric discharge, and urban settlement.The modeled aggregate forms when particles are released one at a time into a volume of space and, influenced by random thermal motion, they diffuse throughout the volume. There is a finite probability that the short-range attraction between particles will influence the motion. Two particles which come into contact with each other will stick together and form a larger unit. The probability of sticking increases as clusters of occupied sites form in the aggregate, stimulating further growth. Simulate this process in 2D using Monte Carlo methods: Create a 2D grid and introduce particles to the lattice through a launching zone one at a time. After a particle is launched, it wanders throughout with a random walk until it either sticks to the aggregate or wanders off the lattice into the kill zone. If a wandering particle enters an empty site next to an occupied site, then the particle's current location automatically becomes part of the aggregate. Otherwise, the random walk continues. Repeat this process until the aggregate contains some pre-determined number of particles. Reference: Wong, Samuel, Computational Methods in Physics and Engineering, 1992.
Program DLA.java simulates the growth of a DLA with the following properties. It uses the helper data type Picture.java. Set the initial aggregate to be the bottom row of the N-by-N lattice. Launch the particles from a random cell in top row. Assume that the particle goes up with probability 0.15, down with probability 0.35, and left or right with probability 1/4 each. Continue until the particles stick to a neighboring cell (above, below, left, right, or one of the four diagonals) or leaves the N-by-N lattice. The preferred downward direction is analogous to the effect of a temperature gradient on Brownian motion, or like how when a crystal is formed, the bottom of the aggregate is cooled more than the top; or like the influence of a gravitational force. For effect, we color the particles in the order they are released according to the rainbow from red to violet. Below are three simulations with N = 176; here is an image with N = 600.
| 
| 
|
Brownian motion. Brownian motion is a random process used to model a wide variety of physical phenomenon including the dispersion of ink flowing in water, and the behavior of atomic particles predicted by quantum physics. (more applications). Fundamental random process in the universe. It is the limit of a discrete random walk and the stochastic analog of the Gaussian distribution. It is now widely used in computational finance, economics, queuing theory, engineering, robotics, medical imaging, biology, and flexible manufacturing systems. First studied by a Scottish botanist Robert Brown in 1828 and analyzed mathematically by Albert Einstein in 1905. Jean-Baptiste Perrin performed experiments to confirm Einstein's predictions and won a Nobel Prize for his work. An applet to illustrate physical process that may govern cause of Brownian motion.
Simulating a Brownian motion. Since Brownian motion is a continuous and stochastic process, we can only hope to plot one path on a finite interval, sampled at a finite number of points. We can interpolate linearly between these points (i.e., connect the dots). For simplicitly, we'll assume the interval is from 0 to 1 and the sample points t0, t1, ..., tN are equally spaced in this interval. To simulate a standard Brownian motion, repeatedly generate independent Gaussian random variables with mean 0 and standard deviation sqrt(1/N). The value of the Brownian motion at time i is the sum of the first i increments.
| 
| 
|
Geometric Brownian motion. A variant of Brownian motion is widely used to model stock prices, and the Nobel-prize winning Black-Scholes model is centered on this stochastic process. A geometric Brownian motion with drift μ and volatility σ is a stochastic process that can model the price of a stock. The parameter μ models the percentage drift. If μ = 0.10, then we expect the stock to increase by 10% each year. The parameter σ models the percentage volatility. If σ = 0.20, then the standard deviation of the stock price over one year is roughly 20% of the current stock price. To simulate a geometric Brownian motion from time t = 0 to t = T, we follow the same procedure for standard Brownian motion, but multiply the increments, instead of adding them, and incorporate the drift and volatility parameters. Specifically, we multiply the current price by by (1 + μΔt + σsqrt(Δt)Z), where Z is a standard Gaussian and Δt = T/N Start with X(0) = 100, σ = 0.04.
Black-Scholes formula. Move to here?
Ising model. The motions of electrons around a nucleus produce a magnetic field associated with the atom. These atomic magnets act much like conventional magnets. Typically, the magnets point in random directions, and all of the forces cancel out leaving no overall magnetic field in a macroscopic clump of matter. However, in some materials (e.g., iron), the magnets can line up producing a measurable magnetic field. A major achievement of 19th century physics was to describe and understand the equations governing atomic magnets. The probability that state S occurs is given by the Boltzmann probability density function P(S) = e-E(S)/kT / Z, where Z is the normalizing constant (partition function) sum e-E(A)/kT over all states A, k is Boltzmann's constant, T is the absolute temperature (in degrees Kelvin), and E(S) is the energy of the system in state S.
Ising model proposed to describe magnetism in crystalline materials. Also models other naturally occurring phenomena including: freezing and evaporation of liquids, protein folding, and behavior of glassy substances.
Ising model. The Boltzmann probability function is an elegant model of magnetism. However, it is not practical to apply it for calculating the magnetic properties of a real iron magnet because any macroscopic chunk of iron contains an enormous number atoms and they interact in complicated ways. The Ising model is a simplified model for magnets that captures many of their important properties, including phase transitions at a critical temperature. (Above this temperature, no macroscopic magnetism, below it, systems exhibits magnetism. For example, iron loses its magnetization around 770 degrees Celsius. Remarkable thing is that transition is sudden.) reference
First introduced by Lenz and Ising in the 1920s. In the Ising model, the iron magnet is divided into an N-by-N grid of cells. (Vertex = atom in crystal, edge = bond between adjacent atoms.) Each cell contains an abstract entity known as spin. The spin si of cell i is in one of two states: pointing up (+1) or pointing down (-1). The interactions between cells is limited to nearest neighbors. The total magnetism of the system M = sum of si. The total energy of the system E = sum of - J si sj, where the sum is taken over all nearest neighbors i and j. The constant J measures the strength of the spin-spin interactions (in units of energy, say ergs). [The model can be extended to allow interaction with an external magnetic field, in which case we add the term -B sum of sk over all sites k.] If J > 0, the energy is minimized when the spins are aligned (both +1 or both -1) - this models ferromagnetism. if J < 0, the energy is minimized when the spins are oppositely aligned - this models antiferromagnetism.
Given this model, a classic problem in statistical mechanics is to compute the expected magenetism. A state is the specification of the spin for each of the N^2 lattice cells. The expected magnetism of the system E[M] = sum of M(S) P(S) over all states S, where M(S) is the magnetism of state S, and P(S) is the probability of state S occurring according to the Boltzmann probability function. Unfortunately, this equation is not amenable to a direct computational because the number of states S is 2N*N for an N-by-N lattice. Straightforward Monte Carlo integration won't work because random points will not contribute much to sum. Need selective sampling, ideally sample points proportional to e-E/kT. (In 1925, Ising solved the problem in one dimension - no phase transition. In a 1944 tour de force, Onsager solved the 2D Ising problem exactly. His solution showed that it has a phase transition. Not likely to be solved in 3D - see intractability section.)
Metropolis algorithm. Widespread usage of Monte Carlo methods began with Metropolis algorithm for calculation of rigid-sphere system. Published in 1953 after dinner conversation between Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller. Widely used to study equilibrium properties of a system of atoms. Sample using Markov chain using Metropolis' rule: transition from A to B with probability 1 if Δ E <= 0, and with probability e-ΔE/kT if Δ E > 0. When applied to the Ising model, this Markov chain is ergodic (similar to Google PageRank requirement) so the theory underlying the Metropolis algorithm applies. Converges to stationary distribution.
Program Cell.java, State.java, and Metropolis.java implements the Metropolis algorithm for a 2D lattice. Ising.java is a procedural programming version. "Doing physics by tossing dice." Simulate complicated physical system by a sequence of simple random steps.
Measuring physical quantities.
Measure magnetism, energy, specific heat when system has thermalized
(the system has reached a thermal equilibrium with its surrounding
environment at a common temperature T).
Compute the average energy Phase transition.
Phase transition occurs when temperature Tc is
2 / ln(1 + sqrt(2)) = 2.26918).
Tc is known as the Curie temperature.
Plot magnetization M (average of all spins) vs. temperature (kT = 1 to 4).
Discontinuity of slope is signature of second order phase transition.
Slope approaches infinity.
Plot energy (average of all spin-spin interactions)
vs. temperature (kT = 1 to 4). Smooth curve through
phase transition.
Compare against
exact solution.
Critical temperature for which algorithm dramatically slows down.
Below are the 5000th sample trajectory
for J/kT = 0.4 (hot / disorder) and 0.47 (cold / order).
The system becomes magnetic as temperature decreases; moreover,
as temperature decreases the probability that neighboring sites
have the same spin increasing (more clumping).
Experiments.
Exact solution for Ising model known for 1D and 2D; NP-hard for 3d
and nonplanar graphs.
Models phase changes in binary alloys and spin glasses.
Also models neural networks, flocking birds, and beating heart cells.
Over 10,000+ papers published using the Ising model.
That is, can you find a real number x and
an integer N for which r equals N?
Solution: No, it can't happen in IEEE floating point
arithmetic. The roundoff error will not cause the result
to be N, even if Math.random returns 0.9999999999....
However, this method does not produce integers
uniformly at random because floating point numbers are
not evenly distributed. Also, it involves casting and multiplying,
which are excessive.
Remark: although the boundary value problem above can be solved analytically,
numerical simulations like the one above are useful when the region has a more
complicated shape or needs to be repeated for different boundary conditions.
simulated annealing
Otherwise repeat.
has the desired distribution. Use Exercise xyz from Section 3 to compute
standard normal deviates.
The color intensities Ip and Is are used to determine
the saturation in the HSB color format.
Program Rainbow.java simulates this process.
Q + A
Exercises
Creative Exercises
double x = Math.random();
int r = (int) (x * N);
N = ceil(ln U / ln (1 - p))
where U is a variable with the uniform distribution. Use the
Math library methods Math.ceil, Math.log, and
Math.random.
public static int poisson(double c) {
double t = 0.0;
for (int x = 0; true; x++) {
t = t - Math.log(Math.random()) / c; // sum exponential deviates
if (t > 1.0) return x;
}
}
(x1)2 + ... + (xN)2 ≤ 1
( x1/r, x2/r, ..., xN/r ), where r = sqrt((x1)2 + ... + (xN)2)
Ip = 1/2 (s(1-s)2 + p(1-p)2)
Is = 1/2 (s2(1-s)2 + p2(1-p)2)
p = (sin(θi-θr)/sin(θi+θr))2
r = (tan(θi-θr)/tan(θi+θr))2
