6.3 Monte Carlo Simulations
Key terms: stochastic process, Monte Carlo simulation, random walk, expected value
6.3.1 Basic Concepts
A stochastic process is one that is governed at least partially by chance. Familiar examples include traffic patterns, weather systems, stock market behavior, and the spread of infectious diseases. Scientists in many disciplines use computer simulations to study processes, events, and systems that are too complicated to analyze mathematically. A Monte Carlo simulation is a program that estimates numerical properties of a stochastic process by sampling random inputs and averaging the results over many trials.
To illustrate the idea with a concrete example, consider a simple game of chance: you roll three dice and win if one of the rolled numbers equals the sum of the other two. What is the probability of winning? If you played the game many times, you might start to get a rough sense of the probability, but for an accurate approximation you would need to play millions of times. The percentage of wins over millions of trials would closely approximate the actual probability of winning. This is exactly how a Monte Carlo simulation for this problem would work.
Listing 6.3.1 is a Monte Carlo simulation that estimates the probability of winning the game just described. If the game were simulated a small number of times, say 100 or 1000, the output would not be consistent across different executions of the program. But with 100 million trials, the output is consistently between 20.83% and 20.84%, suggesting that the true probability is within this narrow range.
Listing 6.3.1 - SumGame.java
package chap06.sect3;
import java.util.concurrent.ThreadLocalRandom;
/**
* Estimates the probability of winning the Sum Game by Monte Carlo simulation. In this game, the
* player rolls three dice and wins if one of the rolled numbers equals the sum of the other two.
*
* @author Drue Coles
*/
public class SumGame {
public static void main(String[] args) {
final int numGames = 100_000_000;
int wins = 0;
System.out.printf("Simulating %,d trials of the Sum Game... %n", numGames);
for (int i = 0; i < numGames; i++) {
if (playerWins()) {
wins++;
}
}
double probability = (double) wins / numGames * 100;
System.out.printf("Estimated probability of winning: %.3f%% %n", probability);
}
/**
* Simulates the Sum Game.
*
* @return true if the player wins and false otherwise
*/
private static boolean playerWins() {
ThreadLocalRandom rand = ThreadLocalRandom.current();
int a = rand.nextInt(1, 7);
int b = rand.nextInt(1, 7);
int c = rand.nextInt(1, 7);
return a == b + c || b == a + c || c == a + b;
}
}
Output 6.3.1
Simulating 100,000,000 trials of the Sum Game...
Estimated probability of winning: 20.827%
6.3.2 Case Study: Random Walks
A random walk (in two dimensions) is the path traced by a point in the plane repeatedly moving one unit in a random direction. Listing 6.3.2 is a Monte Carlo simulation that estimates the expected length of a random walk from the center of a circle to its boundary. The expected length is an example of an expected value in probability theory — intuitively, the average outcome of repeated trials of a probabilistic experiment.
As a simple illustration, suppose you flip a coin and win a dollar (heads) or nothing (tails). Your expected profit is 50 cents because if you played the game many times, your average profit would approach 50 cents. After a million plays, the average might be 49.8357 cents; after a billion plays, it would likely be even closer. Similarly, the expected length of a random walk is the number of steps that would be taken on average over many independent random walks.
Listing 6.3.2 - RandomWalk.java
package chap06.sect3;
import java.util.Scanner;
import java.util.concurrent.ThreadLocalRandom;
/**
* Estimates by Monte Carlo simulation the expected number of steps for a random walk that starts
* at the center of a circle of given radius and continues until it exits the circle.
*
* @author Drue Coles
*/
public class RandomWalk {
public static void main(String[] args) {
Scanner in = new Scanner(System.in);
System.out.print("Enter radius of circle: ");
int radius = in.nextInt();
final int numWalks = 100_000;
int numSteps = 0; // total over all walks
// Monte Carlo simulation
System.out.printf("Simulating %,d random walks... %n", numWalks);
for (int i = 0; i < numWalks; i++) {
numSteps += randomWalk(radius);
}
double expectedSteps = (double) numSteps / numWalks;
System.out.printf("Expected length of random walk: %,d steps. %n", (int) expectedSteps);
}
/**
* Simulates a random walk that starts at the center of a circle of given radius and continues
* until it exits the circle.
*
* @return the number of steps
*/
private static int randomWalk(int radius) {
int steps = 0;
int x = 0; // x-coordinate of current position
int y = 0; // y-coordinate of current position
ThreadLocalRandom rand = ThreadLocalRandom.current();
// take random steps while the point remains in the circle
while (x * x + y * y <= radius * radius) {
switch (rand.nextInt(4)) { // move one unit in a random cardinal direction
case 0 -> x++;
case 1 -> x--;
case 2 -> y++;
case 3 -> y--;
}
steps++;
}
return steps;
}
}
Output 6.3.2
Enter radius of circle: 100
Simulating 100,000 random walks...
Expected length of random walk: 10,127 steps.
6.3.3 Case Study: Approximating π
It is perhaps surprising that Monte Carlo simulations can be used to approximate a purely geometric quantity such as π. This is shown in Listing 6.3.3; see the class documentation for the underlying idea.
Listing 6.3.3 - PiApproximator.java
package chap06.sect3;
import java.util.concurrent.ThreadLocalRandom;
/**
* Approximates the value of π using a Monte Carlo simulation. A random point (x, y) is chosen from
* the unit square. The probability that the point lies inside the quarter circle of radius 1 is
* π/4. Therefore, π ≈ 4 × (fraction of points inside the quarter circle).
*
* @author Drue Coles
*/
public class PiApproximator {
public static void main(String[] args) {
final int numPoints = 100_000_000;
int pointsInCircle = 0;
ThreadLocalRandom rand = ThreadLocalRandom.current();
System.out.printf("Generating %,d random points from the unit square... %n", numPoints);
for (int i = 0; i < numPoints; i++) {
// random point (x, y) in unit square
double x = rand.nextDouble();
double y = rand.nextDouble();
if (x * x + y * y < 1) { // inside quarter circle
pointsInCircle++;
}
}
double piApprox = 4.0 * pointsInCircle / numPoints;
// A double-precision floating-point number has about 16 decimal digits of precision, but
// the last digit may be unreliable due to binary rounding, so 15 is specified.
String label1 = "The value of π estimated by Monte Carlo simulation";
String label2 = "Double-precision floating-point value nearest to π";
System.out.printf("%s: %.15f %n", label1, piApprox);
System.out.printf("%s: %.15f %n", label2, Math.PI);
}
}
Output 6.3.3
Generating 100,000,000 random points from the unit square...
The value of π estimated by Monte Carlo simulation: 3.141841600000000
Double-precision floating-point value nearest to π: 3.141592653589793
This technique converges more slowly than the Leibniz series approximation (Listing 6.2.1).
6.3.4 Case Study: Craps Probability
Listing 6.3.4a is a Monte Carlo simulation for approximating the probability of winning at Craps. The game logic is implemented as a helper method that returns a Boolean value indicating whether the player wins or loses.
Listing 6.3.4a - CrapsProbabilityCalculator.java
package chap06.sect3;
import java.util.concurrent.ThreadLocalRandom;
/**
* Estimates by Monte Carlo simulation the probability of winning the game of Craps.
*
* @author Drue Coles
*/
public class CrapsProbabilityCalculator {
public static void main(String[] args) {
final int numGames = 100_000_000;
int wins = 0;
System.out.printf("Simulating %,d trials of the game of Craps... %n", numGames);
for (int i = 0; i < numGames; i++) {
if (playerWins()) {
wins++;
}
}
double probability = (double) wins / numGames * 100;
System.out.printf("Estimated probability of winning: %.3f %n", probability);
}
/**
* Simulates the game of Craps.
*
* @return true if player wins, false otherwise
*/
public static boolean playerWins() {
final int comeOutRoll = rollDice();
if (comeOutRoll == 7 || comeOutRoll == 11) {
return true;
}
if (comeOutRoll == 2 || comeOutRoll == 3 || comeOutRoll == 12) {
return false;
}
int roll = rollDice();
while (roll != comeOutRoll && roll != 7) {
roll = rollDice();
}
return roll == comeOutRoll;
}
/**
* Simulates rolling a pair of dice.
*
* @return the sum of numbers rolled
*/
public static int rollDice() {
ThreadLocalRandom rand = ThreadLocalRandom.current();
int die1 = rand.nextInt(1, 7);
int die2 = rand.nextInt(1, 7);
return die1 + die2;
}
}
Output 6.3.4a
Simulating 100,000,000 trials of the game of Craps...
Estimated probability of winning: 49.298
Now consider extending the program to also calculate the expected number of rolls in a game. This
introduces a puzzle: it may seem that the playCraps method will now need to return two values
(the outcome of the game and the number of rolls), but in Java a method can only return a single
value. The program could instead run two separate Monte Carlo simulations, one for the probability
of winning and the other for the expected length of a game, but ideally each iteration of the game
would provide a data point for both calculations.
One way to solve this puzzle would be to return some kind of object that encapsulates the two
values. For example, the Boolean-valued outcome and the number of rolls could be combined as a
string; the caller (main) could extract the two pieces of information as substrings. Listing
6.3.4b works along these lines but uses a single number to encode the two values: it returns an
int whose absolute value is the number of rolls and whose sign (positive or negative) indicates
the outcome (win or lose).
Listing 6.3.4b - CrapsProbabilityCalculator2.java
package chap06.sect3;
import java.util.concurrent.ThreadLocalRandom;
/**
* Estimates by Monte Carlo simulation the probability of winning the game of Craps. The program
* also outputs the expected number of rolls per game and greatest number of rolls observed.
*
* @author Drue Coles
*/
public class CrapsProbabilityCalculator2 {
public static void main(String[] args) {
final int numGames = 100_000_000;
int wins = 0;
int rolls = 0; // total across all games
int maxNumRolls = 0; // maximum in a single game
System.out.printf("Simulating %,d trials of the game of Craps... %n", numGames);
for (int i = 0; i < numGames; i++) {
int result = playCraps();
if (result > 0) {
wins++;
}
int positiveResult = Math.abs(result);
rolls += positiveResult;
maxNumRolls = Math.max(maxNumRolls, positiveResult);
}
double probability = (double) wins / numGames * 100;
double expectedLength = (double) rolls / numGames;
System.out.printf("Estimated probability of winning: %.3f%% %n", probability);
System.out.printf("Expected number of rolls per game: %.3f %n", expectedLength);
System.out.printf("Greatest number of rolls: %d %n", maxNumRolls);
}
/**
* Plays the game of Craps.
*
* @return number of rolls (positive = player wins, negative = player loses)
*/
private static int playCraps() {
final int comeOutRoll = rollDice();
if (comeOutRoll == 7 || comeOutRoll == 11) {
return 1; // positive for winning
}
if (comeOutRoll == 2 || comeOutRoll == 3 || comeOutRoll == 12) {
return -1; // negative for losing
}
int numRolls = 1;
do {
int roll = rollDice();
numRolls++;
if (roll == comeOutRoll) { // player wins
return numRolls;
}
if (roll == 7) { // player loses
return -numRolls;
}
} while (true);
}
/**
* Rolls a pair of dice.
*
* @return the sum of the numbers rolled
*/
private static int rollDice() {
ThreadLocalRandom rand = ThreadLocalRandom.current();
int die1 = rand.nextInt(1, 7);
int die2 = rand.nextInt(1, 7);
return die1 + die2;
}
}
Output 6.3.4b
Simulating 100,000,000 trials of the game of Craps...
Estimated probability of winning: 49.293%
Expected number of rolls per game: 3.375
Greatest number of rolls: 63
6.3.5 Case Study: Business Decisions
Monte Carlo simulations have important applications in biology, business, engineering, physics, and other areas. In business, they can be used to reason about the time customers may spend waiting for a service under various assumptions about service time, arrival rate, and other factors that vary unpredictably. Listing 6.3.5 illustrates the idea with a simple queueing model. The class documentation describes the scenario and the expected value being estimated. The code is self-explanatory except for a minor technical point discussed below.
Listing 6.3.5 - WaitingForMassage.java
package chap06.sect3;
import java.util.Scanner;
import java.util.concurrent.ThreadLocalRandom;
/**
* Estimates by Monte Carlo simulation the expected maximum number of people waiting for a massage
* while a single therapist provides massages of fixed duration on a first-come first-served
* basis to people arriving at random intervals.
*
* @author Drue Coles
*/
public class WaitingForMassage {
public static void main(String[] args) {
Scanner in = new Scanner(System.in);
System.out.print("Time for each massage (in minutes): ");
final int duration = in.nextInt();
System.out.print("Average time between arrivals (in minutes): ");
final int arrivalInterval = in.nextInt();
System.out.print("Period of operation (in minutes): ");
final int period = in.nextInt();
System.out.println();
final int trials = 1_000_000;
int totalMaxWaiting = 0;
System.out.printf("Simulating %,d operational periods... %n", trials);
for (int i = 0; i < trials; i++) {
totalMaxWaiting += simulateMassageQueue(duration, arrivalInterval, period);
}
double average = (double) totalMaxWaiting / trials;
System.out.printf("Expected maximum number of people waiting: %.2f %n", average);
}
/**
* Simulates one operational period.
*
* @param duration time in minutes for each massage
* @param arrivalInterval average time in minutes between arrivals
* @param period total operation time in minutes
* @return maximum number of people waiting for a massage at one time
*/
private static int simulateMassageQueue(int duration, int arrivalInterval, int period) {
int peopleWaiting = 0;
int maxPeopleWaiting = 0;
int timeRemaining = 0; // time remaining for current massage
ThreadLocalRandom rand = ThreadLocalRandom.current();
// If there are n minutes on average between arrivals, then there is a 1/n chance of an
// arrival during any given minute.
final double arrivalProbability = 1.0 / arrivalInterval;
// simulate minute-by-minute operation
for (int t = 0; t < period; t++) {
// deduct one minute from time remaining for current massage
if (timeRemaining > 0) {
timeRemaining--;
}
// start new massage if none is in progress and somebody is waiting
if (timeRemaining == 0 && peopleWaiting > 0) {
peopleWaiting--;
timeRemaining = duration;
}
// check for new arrival
if (rand.nextDouble() < arrivalProbability) {
peopleWaiting++;
maxPeopleWaiting = Math.max(peopleWaiting, maxPeopleWaiting);
}
}
return maxPeopleWaiting;
}
}
Output 6.3.5
Time for each massage (in minutes): 10
Average time between arrivals (in minutes): 12
Period of operation (in minutes): 180
Simulating 1,000,000 operational periods...
Expected maximum number of people waiting: 3.20
In the simulateMassageQueue method, the for loop repeats once for each minute of service. In the
body of the loop, it must be determined whether a new customer is arriving for service at the
current time. This is a random event, but the average time between arrivals is known. A k-minute
average is modeled by assuming a probability of 1/k that a new arrival occurs in any given minute.
Specifically, in the if statement checking for a new arrival, nextDouble returns a random
double between 0 and 1, which is compared with the assumed arrival rate of
1.0 / arrivalInterval.