A Concrete Introduction to Probability (using Python)¶

This notebook covers the basics of probability theory, with Python 3 implementations. (You should have some background inprobability andPython.)

In 1814, Pierre-Simon Laplacewrote:

Probability ... is thus simply a fraction whose numerator is the number of favorable cases and whose denominator is the number of all the cases possible ... when nothing leads us to expect that any one of these cases should occur more than any other.

Laplace really nailed it, way back then! If you want to untangle a probability problem, all you have to do is be methodical about defining exactly what the cases are, and then careful in counting the number of favorable and total cases. We'll start being methodical by defining some vocabulary:

• Experiment:An occurrence with an uncertain outcome that we can observe.
  For example, rolling a die.
• Outcome:The result of an experiment; one particular state of the world. What Laplace calls a "case."
  For example:4.
• Sample Space:The set of all possible outcomes for the experiment.
  For example,{1, 2, 3, 4, 5, 6}.
• Event:A subset of possible outcomes that together have some property we are interested in.
  For example, the event "even die roll" is the set of outcomes{2, 4, 6}.
• Probability:As Laplace said, the probability of an event with respect to a sample space is the number of favorable cases (outcomes from the sample space that are in the event) divided by the total number of cases in the sample space. (This assumes that all outcomes in the sample space are equally likely.) Since it is a ratio, probability will always be a number between 0 (representing an impossible event) and 1 (representing a certain event).
  For example, the probability of an even die roll is 3/6 = 1/2.

This notebook will develop all these concepts; I also have asecond part that covers paradoxes in Probability Theory.

Code forP¶

P is the traditional name for the Probability function:

Read this as implementing Laplace's quote directly:"Probability is thus simply a fraction whose numerator is the number of favorable cases and whose denominator is the number of all the cases possible."

Warm-up Problem: Die Roll¶

What's the probability of rolling an even number with a single six-sided fair die?

We can define the sample spaceD and the eventeven, and compute the probability:

It is good to confirm what we already knew.

You may ask: Why does the definition ofP uselen(event & space) rather thanlen(event)? Because I don't want to count outcomes that were specified inevent but aren't actually in the sample space. Consider:

Here,len(event) andlen(space) are both 6, so if just divided, thenP would be 1, which is not right.The favorable cases are theintersection of the event and the space, which in Python is(event & space).Also note that I useFraction rather than regular division because I want exact answers like 1/3, not 0.3333333333333333.

Urn Problems¶

Around 1700, Jacob Bernoulli wrote about removing colored balls from an urn in his landmark treatiseArs Conjectandi, and ever since then, explanations of probability have relied onurn problems. (You'd think the urns would be empty by now.)

For example, here is a three-part problemadapted from mathforum.org:

An urn contains 23 balls: 8 white, 6 blue, and 9 red. We select six balls at random (each possible selection is equally likely). What is the probability of each of these possible outcomes:

1. all balls are red

2. 3 are blue, 2 are white, and 1 is red
3. exactly 4 balls are white

So, an outcome is a set of 6 balls, and the sample space is the set of all possible 6 ball combinations. We'll solve each of the 3 parts using ourP function, and also using basic arithmetic; that is,counting. Counting is a bit tricky because:

• We have multiple balls of the same color.
• An outcome is aset of balls, where order doesn't matter, not asequence, where order matters.

To account for the first issue, I'll have 8 different white balls labelled'W1' through'W8', rather than having eight balls all labelled'W'. That makes it clear that selecting'W1' is different from selecting'W2'.

The second issue is handled automatically by theP function, but if I want to do calculations by hand, I will sometimes first count the number ofpermutations of balls, then get the number ofcombinations by dividing the number of permutations byc!, wherec is the number of balls in a combination. For example, if I want to choose 2 white balls from the 8 available, there are 8 ways to choose a first white ball and 7 ways to choose a second, and therefore 8 × 7 = 56 permutations of two white balls. But there are only 56 / 2 = 28 combinations, because(W1, W2) is the same combination as(W2, W1).

We'll start by defining the contents of the urn:

Now we can define the sample space,U6, as the set of all 6-ball combinations. We useitertools.combinations to generate the combinations, and then join each combination into a string:

I don't want to print all 100,947 members of the sample space; let's just peek at a random sample of them:

Is 100,947 really the right number of ways of choosing 6 out of 23 items, or "23 choose 6", as mathematicianscall it? Well, we can choose any of 23 for the first item, any of 22 for the second, and so on down to 18 for the sixth. But we don't care about the ordering of the six items, so we divide the product by 6! (the number of permutations of 6 things) giving us:

$$23 ~\mbox{choose}~ 6 = \frac{23 \cdot 22 \cdot 21 \cdot 20 \cdot 19 \cdot 18}{6!} = 100947$$

Note that $23 \cdot 22 \cdot 21 \cdot 20 \cdot 19 \cdot 18 = 23! \;/\; 17!$, so, generalizing, we can write:

$$n ~\mbox{choose}~ c = \frac{n!}{(n - c)! \cdot c!}$$

And we can translate that to code and verify that 23 choose 6 is 100,947:

Now we're ready to answer the 4 problems:

Urn Problem 1: what's the probability of selecting 6 red balls?¶

Let's investigate a bit more. How many ways of getting 6 red balls are there?

Why are there 84 ways? Because there are 9 red balls in the urn, and we are asking how many ways we can choose 6 of them:

So the probabilty of 6 red balls is then just 9 choose 6 divided by the size of the sample space:

Urn Problem 2: what is the probability of 3 blue, 2 white, and 1 red?¶

We can get the same answer by counting how many ways we can choose 3 out of 6 blues, 2 out of 8 whites, and 1 out of 9 reds, and dividing by the number of possible selections:

Here we don't need to divide by any factorials, becausechoose has already accounted for that.

We can get the same answer by figuring: "there are 6 ways to pick the first blue, 5 ways to pick the second blue, and 4 ways to pick the third; then 8 ways to pick the first white and 7 to pick the second; then 9 ways to pick a red. But the order'B1, B2, B3' should count as the same as'B2, B3, B1' and all the other orderings; so divide by 3! to account for the permutations of blues, by 2! to account for the permutations of whites, and by 100947 to get a probability:

Urn Problem 3: What is the probability of exactly 4 white balls?¶

We can interpret this as choosing 4 out of the 8 white balls, and 2 out of the 15 non-white balls. Then we can solve it the same three ways:

Revised Version ofP, with more general events¶

To calculate the probability of an even die roll, I originally said

even = {2, 4, 6}

But that's inelegant—I had to explicitly enumerate all the even numbers from one to six. If I ever wanted to deal with a twelve or twenty-sided die, I would have to go back and changeeven. I would prefer to defineeven once and for all like this:

Now in order to makeP(even, D) work, I'll have to modifyP to accept an event as eitheraset of outcomes (as before), or apredicate over outcomes—a function that returns true for an outcome that is in the event:

Here we see howsuch_that, the neweven predicate, and the newP work:

Note:such_that is just like the built-in functionfilter, exceptsuch_that returns a set.

We can now define more interesting events using predicates; for example we can determine the probability that the sum of a three-dice roll is prime (using a definition ofis_prime that is efficient enough for smalln):

Card Problems¶

Consider dealing a hand of five playing cards. We can definedeck as a set of 52 cards, andHands as the sample space of all combinations of 5 cards:

Now we can answer questions like the probability of being dealt a flush (5 cards of the same suit):

Or the probability of four of a kind:

Fermat and Pascal: Gambling, Triangles, and the Birth of Probability¶

Pierre de Fermat 1654

Blaise Pascal] 1654

Consider a gambling game consisting of tossing a coin. Player H wins the game if 10 heads come up, and T wins if 10 tails come up. If the game is interrupted when H has 8 heads and T has 7 tails, how should the pot of money (which happens to be 100 Francs) be split?In 1654, Blaise Pascal and Pierre de Fermat corresponded on this problem, with Fermatwriting:

Dearest Blaise,

As to the problem of how to divide the 100 Francs, I think I have found a solution that you will find to be fair. Seeing as I needed only two points to win the game, and you needed 3, I think we can establish that after four more tosses of the coin, the game would have been over. For, in those four tosses, if you did not get the necessary 3 points for your victory, this would imply that I had in fact gained the necessary 2 points for my victory. In a similar manner, if I had not achieved the necessary 2 points for my victory, this would imply that you had in fact achieved at least 3 points and had therefore won the game. Thus, I believe the following list of possible endings to the game is exhaustive. I have denoted 'heads' by an 'h', and tails by a 't.' I have starred the outcomes that indicate a win for myself.

h h h h *       h h h t *       h h t h *       h h t t *h t h h *       h t h t *       h t t h *       h t t tt h h h *       t h h t *       t h t h *       t h t tt t h h *       t t h t         t t t h         t t t t

I think you will agree that all of these outcomes are equally likely. Thus I believe that we should divide the stakes by the ration 11:5 in my favor, that is, I should receive (11/16)*100 = 68.75 Francs, while you should receive 31.25 Francs.

I hope all is well in Paris,

Your friend and colleague,

Pierre

Pascal agreed with this solution, andreplied with a generalization that made use of his previous invention, Pascal's Triangle. There's evena book about it.

We can solve the problem with the tools we have:

Our answer agrees with Pascal and Fermat; we're in good company!

Non-Equiprobable Outcomes: Probability Distributions¶

So far, we have made the assumption that every outcome in a sample space is equally likely. In real life, we often get outcomes that are not equiprobable. For example, the probability of a child being a girl is not exactly 1/2, and the probability is slightly different for a second child. Anarticle gives the following counts for two-child families in Denmark, whereGB means a family where the first child is a girl and the second a boy:

GG: 121801    GB: 126840BG: 127123    BB: 135138

We will introduce three more definitions:

• Frequency: a number describing how often an outcome occurs. Can be a count like 121801, or a ratio like 0.515.

• Distribution: A mapping from outcome to frequency for each outcome in a sample space.

• Probability Distribution: A distribution that has beennormalized so that the sum of the frequencies is 1.

We defineProbDist to take the same kinds of arguments thatdict does: either a mapping or an iterable of(key, val) pairs, and/or optional keyword arguments.

We also need to modify the functionsP andsuch_that to accept either a sample space or a probability distribution as the second argument.

Here is the probability distribution for Danish two-child families:

And here are some predicates that will allow us to answer some questions:

The above says that the probability of a girl is somewhere between 48% and 49%, but that it is slightly different between the first or second child.

The above says that the sex of the second child is more likely to be the same as the first child, by about 1/2 a percentage point.

More Urn Problems: M&Ms and Bayes¶

Here's another urn problem (or "bag" problem)from prolific Python/Probability authorAllen Downey:

The blue M&M was introduced in 1995. Before then, the color mix in a bag of plain M&Ms was (30% Brown, 20% Yellow, 20% Red, 10% Green, 10% Orange, 10% Tan). Afterward it was (24% Blue , 20% Green, 16% Orange, 14% Yellow, 13% Red, 13% Brown).A friend of mine has two bags of M&Ms, and he tells me that one is from 1994 and one from 1996. He won't tell me which is which, but he gives me one M&M from each bag. One is yellow and one is green. What is the probability that the yellow M&M came from the 1994 bag?

To solve this problem, we'll first represent probability distributions for each bag:bag94 andbag96:

Next, defineMM as the joint distribution—the sample space for picking one M&M from each bag. The outcome'yellow green' means that a yellow M&M was selected from the 1994 bag and a green one from the 1996 bag.

First we'll look at the "One is yellow and one is green" part:

Now we can answer the question: given that we got a yellow and a green (but don't know which comes from which bag), what is the probability that the yellow came from the 1994 bag?

So there is a 74% chance that the yellow comes from the 1994 bag.

Answering this question was straightforward: just like all the other probability problems, we simply create a sample space, and useP to pick out the probability of the event in question, given what we know about the outcome.But in a sense it is curious that we were able to solve this problem with the same methodology as the others: this problem comes from a section titledMy favorite Bayes's Theorem Problems, so one would expect that we'd need to invoke Bayes Theorem to solve it. The computation above shows that that is not necessary.

Of course, wecould solve it using Bayes Theorem. Why is Bayes Theorem recommended? Because we are asked about the probability of an event given the evidence, which is not immediately available; however the probability of the evidence given the event is.

Before we see the colors of the M&Ms, there are two hypotheses,A andB, both with equal probability:

A: first M&M from 94 bag, second from 96 bagB: first M&M from 96 bag, second from 94 bagP(A) = P(B) = 0.5

Then we get some evidence:

E: first M&M yellow, second green

We want to know the probability of hypothesisA, given the evidence:

P(A | E)

That's not easy to calculate (except by enumerating the sample space). But Bayes Theorem says:

P(A | E) = P(E | A) * P(A) / P(E)

The quantities on the right-hand-side are easier to calculate:

P(E | A) = 0.20 * 0.20 = 0.04P(E | B) = 0.10 * 0.14 = 0.014P(A)     = 0.5P(B)     = 0.5P(E)     = P(E | A) * P(A) + P(E | B) * P(B)          = 0.04     * 0.5  + 0.014    * 0.5   =   0.027

And we can get a final answer:

P(A | E) = P(E | A) * P(A) / P(E)          = 0.04     * 0.5  / 0.027          = 0.7407407407

You have a choice: Bayes Theorem allows you to do less calculation at the cost of more algebra; that is a great trade-off if you are working with pencil and paper. Enumerating the state space allows you to do less algebra at the cost of more calculation; often a good trade-off if you have a computer. But regardless of the approach you use, it is important to understand Bayes theorem and how it works.

There is one important question that Allen Downey does not address:would you eat twenty-year-old M&Ms?😨

Newton's Answer to a Problem by Pepys¶

Isaac Newton 1693

Samuel Pepys 1693

This paper explains how Samuel Pepys wrote to Isaac Newton in 1693 to pose the problem:

Which of the following three propositions has the greatest chance of success?

1. Six fair dice are tossed independently and at least one “6” appears.
2. Twelve fair dice are tossed independently and at least two “6”s appear.
3. Eighteen fair dice are tossed independently and at least three “6”s appear.

Newton was able to answer the question correctly (although his reasoning was not quite right); let's see how we can do. Since we're only interested in whether a die comes up as "6" or not, we can define a single die and the joint distribution overn dice as follows:

Now we are ready to determine which proposition is more likely to have the required number of sixes:

We reach the same conclusion Newton did, that the best chance is rolling six dice.

Simulation¶

Sometimes it is inconvenient to explicitly define a sample space. Perhaps the sample space is infinite, or perhaps it is just very large and complicated, and we feel more confident in writing a program tosimulate one pass through all the complications, rather than try toenumerate the complete sample space.Random sampling from the simulationcan give an accurate estimate of the probability.

Simulating Monopoly¶

Considerproblem 84 from the excellentProject Euler, which asks for the probability that a player in the game Monopoly ends a roll on each of the squares on the board. To answer this we need to take into account die rolls, chance and community chest cards, and going to jail (from the "go to jail" space, from a card, or from rolling doubles three times in a row). We do not need to take into account anything about buying or selling properties or exchanging money or winning or losing the game, because these don't change a player's location. We will assume that a player in jail will always pay to get out of jail immediately.

A game of Monopoly can go on forever, so the sample space is infinite. But even if we limit the sample space to say, 1000 rolls, there are $21^{1000}$ such sequences of rolls (and even more possibilities when we consider drawing cards). So it is infeasible to explicitly represent the sample space.

But it is fairly straightforward to implement a simulation and run it for, say, 400,000 rolls (so the average square will be landed on 10,000 times). Here is the code for a simulation:

And the results:

I'll show a histogram of the squares, with a dotted red line at the average:

Another way to see the results:

There is one square far above average:JAIL, at a little over 6%. There are four squares far below average: the three chance squares,CH1,CH2, andCH3, at around 1% (because 10 of the 16 chance cards send the player away from the square), and the "Go to Jail" square, square number 30 on the plot, which has a frequency of 0 because you can't end a turn there. The other squares are around 2% to 3% each, which you would expect, because 100% / 40 = 2.5%.

The Central Limit Theorem / Strength in Numbers Theorem¶

So far, we have talked of anoutcome as being a single state of the world. But it can be useful to break that state of the world down into components. We call these componentsrandom variables. For example, when we consider an experiment in which we roll two dice and observe their sum, we could model the situation with two random variables, one for each die. (Our representation of outcomes has been doing that implicitly all along, when we concatenate two parts of a string, but the concept of a random variable makes it official.)

TheCentral Limit Theorem states that if you have a collection of random variables and sum them up, then the larger the collection, the closer the sum will be to anormal distribution (also called aGaussian distribution or abell-shaped curve). The theorem applies in all but a few pathological cases.

As an example, let's take 5 random variables reprsenting the per-game scores of 5 basketball players, and then sum them together to form the team score. Each random variable/player is represented as a function; calling the function returns a single sample from the distribution:

And here is a function to sample a random variablek times, show a histogram of the results, and return the mean:

The two top-scoring players have scoring distributions that are slightly skewed from normal:

The next two players have bi-modal distributions; some games they score a lot, some games not:

The fifth "player" (actually the sum of all the other players on the team) looks like this:

Now we define the team score to be the sum of the five players, and look at the distribution:

Sure enough, this looks very much like a normal distribution. The Central Limit Theorem appears to hold in this case. But I have to say "Central Limit" is not a very evocative name, so I propose we re-name this as theStrength in Numbers Theorem, to indicate the fact that if you have a lot of numbers, you tend to get the expected result.

Conclusion¶

We've had an interesting tour and met some giants of the field: Laplace, Bernoulli, Fermat, Pascal, Bayes, Newton, ... even Mr. Monopoly and The Count.

The conclusion is: be explicit about what the problem says, and then methodical about defining the sample space, and finally be careful in counting the number of outcomes in the numerator and denominator. Easy as 1-2-3.

Appendix: Continuous Sample Spaces¶

Everything up to here has been about discrete, finite sample spaces, where we canenumerate all the possible outcomes.

But I was asked aboutcontinuous sample spaces, such as the space of real numbers. The principles are the same: probability is still the ratio of the favorable cases to all the cases, but now instead ofcounting cases, we have to (in general) compute integrals to compare the sizes of cases.Here we will cover a simple example, which we first solve approximately by simulation, and then exactly by calculation.

The Hot New Game Show Problem: Simulation¶

Oliver Roeder posedthis problem in the 538Riddler blog:

Two players go on a hot new game show calledHigher Number Wins. The two go into separate booths, and each presses a button, and a random number between zero and one appears on a screen. (At this point, neither knows the other’s number, but they do know the numbers are chosen from a standard uniform distribution.) They can choose to keep that first number, or to press the button again to discard the first number and get a second random number, which they must keep. Then, they come out of their booths and see the final number for each player on the wall. The lavish grand prize — a case full of gold bullion — is awarded to the player who kept the higher number. Which number is the optimal cutoff for players to discard their first number and choose another? Put another way, within which range should they choose to keep the first number, and within which range should they reject it and try their luck with a second number?

We'll use this notation:

• A,B: the two players.
• A,B: the cutoff values they choose: the lower bound of the range of first numbers they will accept.
• a,b: the actual random numbers that appear on the screen.

For example, if playerA chooses a cutoff ofA = 0.6, that means thatA would accept any first number greater than 0.6, and reject any number below that cutoff. The question is: What cutoff,A, should playerA choose to maximize the chance of winning, that is, maximize P(a >b)?

First, simulate the number that a player with a given cutoff gets (note thatrandom.random() returns a float sampled uniformly from the interval [0..1]):

Now compare the numbers returned with a cutoff ofA versus a cutoff ofB, and repeat for a large number of trials; this gives us an estimate of the probability that cutoffA is better than cutoffB:

Now define a function,top, that considers a collection of possible cutoffs, estimate the probability for each cutoff playing against each other cutoff, and returns a list with theN top cutoffs (the ones that defeated the most number of opponent cutoffs), and the number of opponents they defeat:

We get a good idea of the top cutoffs, but they are close to each other, so we can't quite be sure which is best, only that the best is somewhere around 0.60. We could get a better estimate by increasing the number of trials, but that would consume more time.

The Hot New Game Show Problem: Exact Calculation¶

More promising is the possibility of makingPwin(A, B) an exact calculation. But before we get toPwin(A, B), let's solve a simpler problem: assume that both playersA andB have chosen a cutoff, and have each received a number above the cutoff. What is the probability thatA gets the higher number? We'll call thisPhigher(A, B). We can think of this as a two-dimensional sample space of points in the (a,b) plane, wherea ranges from the cutoffA to 1 andb ranges from the cutoff B to 1. Here is a diagram of that two-dimensional sample space, with the cutoffsA=0.5 andB=0.6:

The total area of the sample space is 0.5 × 0.4 = 0.20, and in general it is (1 -A) · (1 -B). What about the favorable cases, whereA beatsB? That corresponds to the shaded triangle below:

The area of a triangle is 1/2 the base times the height, or in this case, 0.4² / 2 = 0.08, and in general, (1 -B)² / 2. So in general we have:

Phigher(A, B) = favorable / totalfavorable     = ((1 - B) ** 2) / 2 total         = (1 - A) * (1 - B)Phigher(A, B) = (((1 - B) ** 2) / 2) / ((1 - A) * (1 - B))Phigher(A, B) = (1 - B) / (2 * (1 - A))

And in this specific case we have:

 A = 0.5; B = 0.6 favorable         = 0.4 ** 2 / 2 = 0.08 total             = 0.5 * 0.4    = 0.20 Phigher(0.5, 0.6) = 0.08 / 0.20  = 0.4

But note that this only works when the cutoffA ≤B; whenA >B, we need to reverse things. That gives us the code:

We're now ready to tackle the full game. There are four cases to consider, depending on whetherA andB gets a first number that is above or below their cutoff choices:

For example, the first row of this table says that the event of both first numbers being above their respective cutoffs has probability (1 -A) · (1 -B), and if this does occur, then the probability ofA winning is Phigher(A,B).We're ready to replace the old simulation-basedPwin with a new calculation-based version:

That was a lot of algebra. Let's define a few tests to check for obvious errors:

Let's repeat the calculation with our new, exactPwin:

It is good to see that the simulation and the exact calculation are in rough agreement; that gives me more confidence in both of them. We see here that 0.62 defeats all the other cutoffs, and 0.61 defeats all cutoffs except 0.62. The great thing about the exact calculation code is that it runs fast, regardless of how much accuracy we want. We can zero in on the range around 0.6:

This says 0.618 is best, better than 0.620. We can get even more accuracy:

So 0.61803 is best. Does that numberlook familiar? Can you prove that it is what I think it is?

To understand the strategic possibilities, it is helpful to draw a 3D plot ofPwin(A, B) for values ofA andB between 0 and 1:

What does thisPringle of Probability show us? The highest win percentage forA, the peak of the surface, occurs whenA is around 0.5 andB is 0 or 1. We can confirm that, finding the maximumPwin(A, B) for many different cutoff values ofA andB:

SoA could win 62.5% of the time if onlyB would chose a cutoff of 0. But, unfortunately forA, a rational playerB is not going to do that. We can ask what happens if the game is changed so that playerA has to declare a cutoff first, and then playerB gets to respond with a cutoff, with full knowledge ofA's choice. In other words, what cutoff shouldA choose to maximizePwin(A, B), given thatB is going to take that knowledge and pick a cutoff that minimizesPwin(A, B)?

And what if we run it the other way around, whereB chooses a cutoff first, and thenA responds?

In both cases, the rational choice for both players in a cutoff of 0.61803, which corresponds to the "saddle point" in the middle of the plot. This is astable equilibrium; consider fixingB = 0.61803, and notice that ifA changes to any other value, we slip off the saddle to the right or left, resulting in a worse win probability forA. Similarly, if we fixA = 0.61803, then ifB changes to another value, we ride up the saddle to a higher win percentage forA, which is worse forB. So neither player will want to move from the saddle point.

The moral for continuous spaces is the same as for discrete spaces: be careful about defining your space; count/measure carefully, and let your code take care of the rest.