Newton-Pepys Problem

This article has an accompanying YouTube video.

In this article, I would like to demonstrate an approach for solving relatively simple programming problems, of the sort you may encounter in a 100-level programming course in college. This technique doesn’t scale well to large applications, but it does seem a good approach to teach to new programming students.

As a basis for this example, I am going to work through a simple solution to an exercise in the programming textbook that I currently use in my introductory programming course. The text is Introduction to Programming in Python by Sedgewick, Wayne, and Dondero. The exercise is 1.3.40, Pepys’s problem. It is as follows,

In 1693 Samuel Pepys asked Isaac Newton which is more likely: getting 1 at least once when rolling a fair die six times or getting 1 at least twice when rolling it twelve times. Compose a program that could have provided Newton with a quick answer. ¹

When looking at a problem, even one as “simple” as this, it is easy to get overwhelmed by details. For this approach, we’re going to take things slowly, one step at a time, in an effort to avoid becoming bogged down in the details.

The basic idea is we are going to immediately start writing code, with minimal planning. We will write many small programs, progressively improving them, until we reach a solution to the problem.

Rolling a Die

From reading the problem, it seems that we will be working with dice. We want to start with as simple a problem as possible, so as a first step let’s figure out how to write a program that “rolls a die” a single time, and displays the result.

For somebody with a small amount of programming experience (or who has actually read the textbook–the example a die-roll program is incredibly common), the solution here will be readily apparent. However, for many people the solution won’t be so obvious.

The core idea here is that we need to isolate what a die “does”, and figure out an analogous action that we can take in Python. Along the way, we will throw out many unnecessary details about dice–a process known as abstraction. I’ve found that this process of mapping a “real-world” action into a corresponding Python action is often the most difficult step for the beginning programmer.

A good way to approach this problem is to start listing out all the properties you can think of of the object that you want to model. So, for a die, we might have,

• has six numbered sides (1-6)
• is a cube
• numbers are represented by dots (usually)
• can be “rolled”, resulting in one side facing up
• each side is equally likely to face up after a roll (because the die is fair)
• can have rounded off corners, or sharp corners
• dots are often inset, and of a different color from most of the cube
• often made of plastic
• etc., the list goes on

Once we have this list, the next step is to try to select from it only those properties that are necessary to answer the question. Revisiting the problem statement, you can see that this question is concerned with the likelihood of rolling certain numbers. So, realistically, all that we care about from that list are,

• has six numbered sides (1-6)
• can be “rolled”, resulting in one side facing up
• each side is equally likely to face up after a roll (because the die is fair)

The other properties listed are properties of the die that aren’t relevant to the problem.

Restating these properties in different terms, we need to write some Python code that will randomly give us a number between 1 and 6 (inclusive), with each number being equally likely. There are a few different ways to do this, but we’ll use random.randint for the task.

% python
>>> import random
>>> help(random.randint)
Help on method randint in module random:

randint(a, b) method of random.Random instance
    Return random integer in range [a, b], including both end points.

So, our basic die-rolling program will simply call this function, specifying a range of $[1,6]$, and write the result to stdout.

# die.py
# simulates a single die roll, 
# writing the result to stdout
import random

roll = random.randint(1, 6)
print(roll)

As an example of using this program,

% ls
die.py
% python die.py
6
% python die.py
3
% python die.py
1

Multiple Die Rolls

Now that we have the actual die roll implemented, the next phase is to write a program that rolls the die multiple times. We need to roll it six times, and then twelve times. Let’s just focus on getting the first case, rolling it six times, working. We can look at the twelve roll case later.

As a first pass, let’s do it the brute force way. If die.py rolls a die a single time, then we can simply call this program six times in a row to get our six rolls.

% ls
die.py
% python die.py
6
% python die.py
1
% python die.py
3
% python die.py
4
% python die.py
3
% python die.py
2

If this process seems overly manual, that is because it is. One of the first rules of programming is not to make the user do work that the program could do for them. So, rather than requiring six calls to our program, we can write a new program that behaves in the same way, by simply copying and pasting the code from die.py six times, like so,

# roll6.py
# simulates six die rolls, 
# writing the results to stdout
import random

roll = random.randint(1, 6)
print(roll)

roll = random.randint(1, 6)
print(roll)

roll = random.randint(1, 6)
print(roll)

roll = random.randint(1, 6)
print(roll)

roll = random.randint(1, 6)
print(roll)

roll = random.randint(1, 6)
print(roll)

Writing this program greatly reduces the workload on the user, because they can now get all six rolls by calling just one program,

% ls
roll6.py
% python roll6.py
4
3
3
5
1
6

Of course, the program itself is still poorly written. Another general rule of programming is that you should avoid repeating yourself wherever possible. In this case, we have the same two lines of code duplicated six times. Whenever you see code duplicated sequentially like in roll6.py, that is a good indication that a loop might be in order. In this case, we will rewrite roll6.py to use a simple for-loop, which will vastly reduce its length,

# roll6v2.py
# simulates six die rolls, 
# writing the results to stdout
import random

for _ in range(6):
    roll = random.randint(1, 6)
    print(roll)

Note that here the loop variable is called _, in keeping with the convention of using _ as the name of a variable which will not be used.

There is another modification that we can make to improve this program. Rather than writing $6$ directly in the call to range in our for loop, let’s define a variable instead. It’s generally a good idea to avoid so called “magic numbers” in programs. Magic numbers are numerical constants that appear inside of expressions, such as the $6$ here.

# roll6v3.py
# simulates six die rolls, 
# writing the results to stdout
import random

roll_count = 6
for _ in range(roll_count):
    roll = random.randint(1, 6)
    print(roll)

While we are removing magic numbers, we have an opportunity to generalize our program a bit. There is nothing special about the number 6 here, other than it being the specific value associated with the problem we are trying to solve. You may want to roll a die any number of times. So, why not allow the user to specify this number? We can use a command-line argument to do just this. The result is,

# roll.py
# simulates die rolls, writing the results
# to stdout. The number of times to roll is
# accepted as a command-line argument.
import random
import sys

roll_count = int(sys.argv[1])
for _ in range(roll_count):
    roll = random.randint(1, 6)
    print(roll)

This program will let us roll the die as many times as we like. For example,

% ls
roll.py
% python roll.py 3
5
3
1
% python roll.py 6
5
5
4
2
3
5

Counting the Ones

Now that we have a program that will roll a die a specified number of times, we need to count up the number of times that a one appears in the output. The problem wants us to determine whether it is more likely to get a single one in six rolls, or a pair of ones in twelve rolls, so we need to have a way of counting the ones.

To accomplish this task, we can define a counter variable above the loop, and increment it each time a one appears. This will require an if-statement. This will result in count_rolls.py,

# count_rolls.py
# simulates die rolls, counting the number of times
# that a 1 is rolled. The number of ones is written
# to stdout.
import random
import sys

roll_count = int(sys.argv[1])
one_count = 0
for _ in range(roll_count):
    roll = random.randint(1, 6)
    print(roll)
    if roll == 1: # a 1 was rolled
        one_count += 1

print(one_count)

Of course, you may notice that we’ve let another magic number creep into our program. So let’s rectify that by making the number we’re looking for (one in this specific problem) a command-line argument.

# count_rollsv2.py
# simulates die rolls, counting the number of times
# that a 1 is rolled. The number of ones is written
# to stdout.
import random
import sys

roll_count = int(sys.argv[1])
target = int(sys.argv[2])
target_count = 0 # we won't always be looking for 1, so let's name this
                 # more generically too
for _ in range(roll_count):
    roll = random.randint(1, 6)
    print(roll)
    if roll == target:
        target_count += 1

print(target_count)

We now have a program that will report the number of appearances of a specified target in the requested number of rolls. For example,

% ls
count_rollsv2.py
% python count_rollsv2.py 6 1
6
6
3
1
3
6
1

In this run, we are counting the number of ones that appear in six rolls. The last number in the output corresponds to this count, and the others are the rolls themselves. It is helpful at this stage to keep the program writing the individual rolls to stdout, as this will allow us to test it. We do need to make sure that the count displayed at the end is accurate, after all!

Calculating the Probability

As our next step, we need to actually calculate probabilities. Again, we are still mostly focused on the six roll case (we’ll handle twelve later), so we need to determine the probability of our target number appearing at least once in our output.

First, we need to actually define what we mean by probability. For our purposes, the following definition will work quite nicely, $$ \text{probability} = \frac{\text{number of times “target” appears at least once}} {\text{number of tests we run}} $$

The only problem is that, for this to work, we need to run our program many times.

Having Python execute our count_rollsv2.py program repeatedly is pretty easy: we can create a new program called pepys.py which contains the code from count_rollsv2.py wrapped in another loop, like so:

# pepys.py
# Repeatedly roll a die a specified number of times and count
# how frequently a target value appears.

import random
import sys

test_count = 100000 # magic numbers bad!
for _ in range(test_count):
    roll_count = int(sys.argv[1])
    target = int(sys.argv[2])
    target_count = 0
    for _ in range(roll_count):
        roll = random.randint(1, 6)
        if roll == target:
            target_count += 1

    print(target_count)

This code will execute count_rollsv2.py 100000 times in a row. Note that the print(roll) statement has been removed at this point, to keep the amount of output down. The selection of 100000 for the number of tests was arbitrary, any sufficiently large value will do.

We can now start modifying this, in order to get our desired output together. From the definition of probability stated above, we need the number of tests that we run and the number of times target_count is at least 1. The first of these is the number of times that the outer loop runs (test_count), and the second is going to require us to add some code.

Additionally, there is some code that we don’t need to have inside of the for loop. The values of target and roll_count do not change as the program runs, so there is no need to refresh their values at the start of every test. Let’s move those outside of the loop so as not to waste a bunch of effort.

# pepysv2.py
# Repeatedly roll a die a specified number of times and determine
# the probability of a target value appearing at least once in
# a given set of rolls.

import random
import sys

roll_count = int(sys.argv[1])
target = int(sys.argv[2])
test_count = 100000

# The number of tests that have met the threshold of target 
# appearing at least once
threshold_met = 0

for _ in range(test_count):
    target_count = 0
    for _ in range(roll_count):
        roll = random.randint(0, 6)
        if roll == target:
            target_count += 1

    if target_count >= 1: # if the threshold is met
        threshold_met += 1

# calculate and display the probability
probability = threshold_met / test_count
print(probability)

Running this program a few times for our desired six rolls per test with a target of one will give,

% ls
pepysv2.py
% python pepysv2.py 6 1
0.66586
% python pepysv2.py 6 1
0.66579
% python pepysv2.py 6 1
0.66685

You may, however, have noticed that a magic number slipped into our code: the threshold. The program has this if statement: if target_count >= 1:, with a hard-coded 1. Let’s fix this, making the threshold a command-line argument as well.

# pepysv3.py
# Repeatedly roll a die a specified number of times and determine
# the probability of a target value appearing at least once in
# a given set of rolls.

import random
import sys

roll_count = int(sys.argv[1])
target = int(sys.argv[2])
threshold = int(sys.argv[3])
test_count = 100000

# The number of tests that have met the threshold of target 
# appearing at least once
threshold_met = 0

for _ in range(test_count):
    target_count = 0
    for _ in range(roll_count):
        roll = random.randint(0, 6)
        if roll == target:
            target_count += 1

    if target_count >= threshold:
        threshold_met += 1

# calculate and display the probability
probability = threshold_met / test_count
print(probability)

We could also make test_count a command-line argument, but our selection of $100000$ works rather well, and we already have three arguments. Let’s leave this as a fixed number for now. I’ll address this number later, in a different article on Monte Carlo analysis.

Answering the Question

Okay, we’re now ready to answer the question. “But wait?” you may well ask, “what about the second case with the twelve rolls? We haven’t addressed that yet.” And that’s true, I have kept pushing that topic off. But, in the process of setting up our program with command-line arguments, we actually solved that problem already! We can address the question at hand, whether you’re more likely to get a single one in six rolls, or a pair in twelve, using two calls to the program,

% ls
pepysv3.py
% pepysv3.py 6 1 1
0.66688
% pepysv3.py 12 1 2
0.61773

So, as it turns out, the first case is slightly more likely than the second!

Conclusion

In this article, we walked through developing the solution to a simple programming problem using an iterative technique. We began by casting the problem in the simplest form possible (roll a die one time), and gradually added complexity until we arrived at a program that solved the actual problem.

At each step, we began by solving a specific question (roll a die once, roll a die six times, etc.) in the most concrete manner possible, even if that meant duplicating code. Once we had a working program, we began to modify it, gradually generalizing it using variables, command-line arguments, and loops. Once it was in a reasonable state, we added the next level of complexity to the problem.

Along the way, we also saw how focusing on a single specific problem, and generalizing it, can automatically yield the solution to other, related problems, without even needing to think about them in advance.

The problem that we were attempting to solve was,

In 1693 Samuel Pepys asked Isaac Newton which is more likely: getting 1 at least once when rolling a fair die six times or getting 1 at least twice when rolling it twelve times. Compose a program that could have provided Newton with a quick answer. ¹

and we addressed it by solving the following, smaller problems, in order.

1. Roll a die a single time and display the result.
2. Roll a die six times and display each result.
3. Roll a die an arbitrary number of times and display each result.
4. Roll a die an arbitrary number of times, and count the number of appearances of the number one.
5. Roll a die an arbitrary number of times, and count the number of appearances of an arbitrary target.
6. Repeat the solution from (5) 100000 times.
7. Repeat the solution from (5) 100000 times, and count the number of times that the target appears one time. Calculate the probability from this.
8. Repeat the solution from (5) 100000 times, and count the number of times that the target appears an arbitrary number of times. Calculate the probability from this.

The final solution to our problem was reached by calling the solution to (8) two times–once for six rolls, and once for twelve.

This approach is a bit laborious, and not well suited for large and complicated problems. However, it does address several “pain points” for beginning programmers. It encourages the examination of specific cases, which are usually easier to get a handle on than trying to generalize from the outset. It also encourages identifying and tackling the simplest possible variation of the problem, and then gradually adding complexity, which can help prevent you from becoming overwhelmed.