<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <div style="text-align
Step1:  Read this as implementing Laplace's quote directly
Step2:  It is good to confirm what we already knew.
You may ask
Step3:  Here, len(event) and len(space) are both 6, so if just divided, then P would be 1, which is not right.
The favorable cases are the intersection of the event and the space, which in Python is (event &amp; space).
Also note that I use Fraction 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 treatise Ars Conjectandi, and ever since then, explanations of probability have relied on urn problems. (You'd think the urns would be empty by now.) 
<center><a href="https
Step4:  Now we can define the sample space, U6, as the set of all 6-ball combinations.  We use itertools.combinations to generate the combinations, and then join each combination into a string
Step5:  I don't want to print all 100,947 members of the sample space; let's just peek at a random sample of them
Step6:  Is 100,947 really the right number of ways of choosing 6 out of 23 items, or  "23 choose 6", as  mathematicians call 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
Step7:  Now we're ready to answer the 4 problems
Step8:  Let's investigate a bit more. How many ways of getting 6 red balls are there?
Step9:  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
Step10:  So the probabilty of 6 red balls is then just 9 choose 6 divided by the size of the sample space
Step11:  Urn Problem 2
Step12:  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
Step13:  Here we don't need to divide by any factorials, because choose has already accounted for that. 
We can get the same answer by figuring
Step14:  Urn Problem 3
Step15:  Revised Version of P, with more general events
To calculate the probability of an even die roll, I originally said
even = {2, 4, 6}
But that's inelegant&mdash;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 change  even.  I would prefer to define even once and for all like this
Step17:  Now in order to make P(even, D) work, I'll have to modify P to accept an event as either
a set of outcomes (as before), or a predicate over outcomes&mdash;a function that returns true for an outcome that is in the event
Step18:  Here we see how such_that, the new even predicate, and the new P work
Step19:  Note
Step20:  Card Problems
Consider dealing a hand of five playing cards. We can define deck as a set of 52 cards, and Hands as the sample space of all combinations of 5 cards
Step21:  Now we can answer questions like the probability of being dealt a flush (5 cards of the same suit)
Step22:  Or the probability of four of a kind
Step23:  Non-Equiprobable Outcomes
Step26:  We also need to modify the functions P and such_that to accept either a sample space  or a probability distribution as the second argument.
Step27:  Here is the probability distribution for Danish two-child families
Step28:  And here are some predicates that will allow us to answer some questions
Step29:  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.
Step30:  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
Step32:  Next, define MM as the joint distribution&mdash;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.
Step33:  First we'll look at the "One is yellow and one is green" part
Step34:  Now we can answer the question
Step37:  So there is a 74% chance that the yellow comes from the 1994 bag.
My favorite Bayes's Theorem Problems
Answering this question is straightforward
Step38:  And the results
Step39:  I'll show a histogram of the squares, with a dotted red line at the average
Step40:  Another way to see the results
Step41:  There is one square far above average
Step42:  And here is a function to sample a random variable k times, show a histogram of the results, and return the mean
Step43:  The two top-scoring players have scoring distributions that are slightly skewed from normal
Step44:  The next two players have bi-modal distributions; some games they score a lot, some games not
Step45:  The fifth "player" represents all the others on the team; the distribution appears to be a cross between a normal and a uniform distribution
Step46:  Now we define the team score to be the sum of the players, and look at the distribution
Step47:  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 the Strength 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 seen how to compute probabilities. Just 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. 
<center><a href="https
Step48:  Now compare the numbers returned with a cutoff of A versus a cutoff of B, and repeat for a large number of trials; this gives us an estimate of the probability that cutoff A is better than cutoff B
Step49:  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 the N top cutoffs (the ones that defeated the most number of opponent cutoffs), and the number of opponents they defeat
Step50:  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
Step51:  We're now ready to tackle the full game. There are four cases to consider, depending on whether A and B gets a first number that is above or below their cutoff choices
Step52:  The exact analysis code is much more complex than the simulation code, relying on a lot of algebra and geometry and case analysis, and thus more likely to have an error. Let's define a few tests to check for obvious errors
Step53:  Let's repeat the calculation with our new, exact Pwin
Step54:  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
Step55:  This says 0.618 is best, better than 0.620. We can get even more accuracy
Step56:  So 0.61803 is best. I don't think we need more accuracy than that.
0.61803 defeats all of the other cutoffs. But in other game-theoretic problems it might be that there is no one strategy that beats all others. Sometimes A beats B and B beats C, but C beats A. To understand what we are dealing with, it is helpful to draw a 3D plot of Pwin(A, B) for values of A and B between 0 and 1
Step57:  What does this show us? The highest win percentage for A, the peak of the surface, occurs when A is around 0.5 and B is 0 or 1. We can confirm that, finding the maximum Pwin(A, B) for many different cutoff values of A and B
Step58:  So A could win 62.5% of the time if only B would chose a cutoff of 0. But, unfortunately for A, a rational player B is not going to do that. We can ask what happens if the game is changed so that player A has to declare a cutoff first, and then player B gets to respond with a cutoff, with full knowledge of A's choice. In other words, what cutoff should A choose to maximize Pwin(A, B), given that B is going to take that knowledge and pick a cutoff that minimizes Pwin(A, B)?
Step59:  And what if we run it the other way around, where B chooses a cutoff first, and then A responds?
<ASSISTANT_TASK:>
Python Code:
from fractions import Fraction
def P(event, space): 
    "The probability of an event, given a sample space of equiprobable outcomes."
    return Fraction(len(event & space), 
                    len(space))
Explanation: <div style="text-align: right">Peter Norvig, 12 Feb 2016</div>
A Concrete Introduction to Probability (using Python)
In this notebook, we cover the basics of probability theory, and show how to implement the theory in Python 3. (You should have a little background in probability and Python.) 
In 1814, Pierre-Simon Laplace wrote (originally in French):
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.
<center><a href="https://en.wikipedia.org/wiki/Pierre-Simon_Laplace">Pierre-Simon Laplace</a></center>
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](https://en.wikipedia.org/wiki/Experiment_(probability_theory%29):
  An occurrence with an uncertain outcome that we can observe.
  <br>For example, rolling a die.
[Outcome](https://en.wikipedia.org/wiki/Outcome_(probability%29):
  The result of an experiment; one particular state of the world. What Laplace calls a "case."
  <br>For example: 4.
Sample Space:
  The set of all possible outcomes for the experiment. 
  <br>For example, {1, 2, 3, 4, 5, 6}.
[Event](https://en.wikipedia.org/wiki/Event_(probability_theory%29):
  A subset of possible outcomes that together have some property we are interested in.
  <br>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).
<br>For example, the probability of an even die roll is 3/6 = 1/2.
This notebook will develop all these concepts; I also have a second part that covers  paradoxes in Probability Theory.
Code for P
P is the traditional name for the Probability function:
End of explanation
D    = {1, 2, 3, 4, 5, 6}
even = {   2,    4,    6}
P(even, D)
Explanation: 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 space D and the event even, and compute the probability:
End of explanation
even = {2, 4, 6, 8, 10, 12}
P(even, D)
Explanation: It is good to confirm what we already knew.
You may ask: Why does the definition of P use len(event &amp; space) rather than len(event)? Because I don't want to count outcomes that were specified in event but aren't actually in the sample space. Consider:
End of explanation
def cross(A, B):
    "The set of ways of concatenating one item from collection A with one from B."
    return {a + b 
            for a in A for b in B}
urn = cross('W', '12345678') | cross('B', '123456') | cross('R', '123456789') 
urn
len(urn)
Explanation: Here, len(event) and len(space) are both 6, so if just divided, then P would be 1, which is not right.
The favorable cases are the intersection of the event and the space, which in Python is (event &amp; space).
Also note that I use Fraction 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 treatise Ars Conjectandi, and ever since then, explanations of probability have relied on urn problems. (You'd think the urns would be empty by now.) 
<center><a href="https://en.wikipedia.org/wiki/Jacob_Bernoulli">Jacob Bernoulli</a></center>
For example, here is a three-part problem adapted  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:
all balls are red
3 are blue, 2 are white, and 1 is red
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 our P 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 a set of balls, where order doesn't matter, not a sequence, 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 the P function, but if I want to do calculations by hand, I will sometimes first count the number of permutations of balls, then get the number of combinations by dividing the number of permutations by c!, where c 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 &times; 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:
End of explanation
import itertools
def combos(items, n):
    "All combinations of n items; each combo as a concatenated str."
    return {' '.join(combo) 
            for combo in itertools.combinations(items, n)}
U6 = combos(urn, 6)
len(U6)
Explanation: Now we can define the sample space, U6, as the set of all 6-ball combinations.  We use itertools.combinations to generate the combinations, and then join each combination into a string:
End of explanation
import random
random.sample(U6, 10)
Explanation: I don't want to print all 100,947 members of the sample space; let's just peek at a random sample of them:
End of explanation
from math import factorial
def choose(n, c):
    "Number of ways to choose c items from a list of n items."
    return factorial(n) // (factorial(n - c) * factorial(c))
choose(23, 6)
Explanation: Is 100,947 really the right number of ways of choosing 6 out of 23 items, or  "23 choose 6", as  mathematicians call 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:
End of explanation
red6 = {s for s in U6 if s.count('R') == 6}
P(red6, U6)
Explanation: Now we're ready to answer the 4 problems: 
Urn Problem 1: what's the probability of selecting 6 red balls?
End of explanation
len(red6)
Explanation: Let's investigate a bit more. How many ways of getting 6 red balls are there?
End of explanation
choose(9, 6)
Explanation: 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:
End of explanation
P(red6, U6) == Fraction(choose(9, 6), 
                        len(U6))
Explanation: So the probabilty of 6 red balls is then just 9 choose 6 divided by the size of the sample space:
End of explanation
b3w2r1 = {s for s in U6 if
          s.count('B') == 3 and s.count('W') == 2 and s.count('R') == 1}
P(b3w2r1, U6)
Explanation: Urn Problem 2: what is the probability of 3 blue, 2 white, and 1 red?
End of explanation
P(b3w2r1, U6) == Fraction(choose(6, 3) * choose(8, 2) * choose(9, 1), 
                          len(U6))
Explanation: 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:
End of explanation
 P(b3w2r1, U6) == Fraction((6 * 5 * 4) * (8 * 7) * 9, 
                           factorial(3) * factorial(2) * len(U6))
Explanation: Here we don't need to divide by any factorials, because choose 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:
End of explanation
w4 = {s for s in U6 if
      s.count('W') == 4}
P(w4, U6)
P(w4, U6) == Fraction(choose(8, 4) * choose(15, 2),
                      len(U6))
P(w4, U6) == Fraction((8 * 7 * 6 * 5) * (15 * 14),
                      factorial(4) * factorial(2) * len(U6))
Explanation: 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:
End of explanation
def even(n): return n % 2 == 0
Explanation: Revised Version of P, with more general events
To calculate the probability of an even die roll, I originally said
even = {2, 4, 6}
But that's inelegant&mdash;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 change  even.  I would prefer to define even once and for all like this:
End of explanation
def P(event, space): 
    The probability of an event, given a sample space of equiprobable outcomes.
    event can be either a set of outcomes, or a predicate (true for outcomes in the event).
    if is_predicate(event):
        event = such_that(event, space)
    return Fraction(len(event & space), len(space))
is_predicate = callable
def such_that(predicate, collection): 
    "The subset of elements in the collection for which the predicate is true."
    return {e for e in collection if predicate(e)}
Explanation: Now in order to make P(even, D) work, I'll have to modify P to accept an event as either
a set of outcomes (as before), or a predicate over outcomes&mdash;a function that returns true for an outcome that is in the event:
End of explanation
such_that(even, D)
P(even, D)
D12 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}
such_that(even, D12)
P(even, D12)
Explanation: Here we see how such_that, the new even predicate, and the new P work:
End of explanation
D3 = {(d1, d2, d3) for d1 in D for d2 in D for d3 in D}
def prime_sum(outcome): return is_prime(sum(outcome))
def is_prime(n): return n > 1 and not any(n % i == 0 for i in range(2, n))
P(prime_sum, D3)
Explanation: Note: such_that is just like the built-in function filter, except such_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 of is_prime that is efficient enough for small n):
End of explanation
suits = 'SHDC'
ranks = 'A23456789TJQK'
deck  = cross(ranks, suits)
len(deck)
Hands = combos(deck, 5)
assert len(Hands) == choose(52, 5)
random.sample(Hands, 5)
Explanation: Card Problems
Consider dealing a hand of five playing cards. We can define deck as a set of 52 cards, and Hands as the sample space of all combinations of 5 cards:
End of explanation
def flush(hand):
    return any(hand.count(suit) == 5 for suit in suits)
P(flush, Hands)
Explanation: Now we can answer questions like the probability of being dealt a flush (5 cards of the same suit):
End of explanation
def four_kind(hand):
    return any(hand.count(rank) == 4 for rank in ranks)
P(four_kind, Hands)
Explanation: Or the probability of four of a kind:
End of explanation
class ProbDist(dict):
    "A Probability Distribution; an {outcome: probability} mapping."
    def __init__(self, mapping=(), **kwargs):
        self.update(mapping, **kwargs)
        # Make probabilities sum to 1.0; assert no negative probabilities
        total = sum(self.values())
        for outcome in self:
            self[outcome] = self[outcome] / total
            assert self[outcome] >= 0
Explanation: 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. An article gives the following counts for two-child families in Denmark, where GB means a family where the first child is a girl and the second a boy:
GG: 121801    GB: 126840
BG: 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 been normalized so that the sum of the frequencies is 1.
We define ProbDist to take the same kinds of arguments that dict does: either a mapping or an iterable of (key, val) pairs, and/or optional keyword arguments.
End of explanation
def P(event, space): 
    The probability of an event, given a sample space of equiprobable outcomes. 
    event: a collection of outcomes, or a predicate that is true of outcomes in the event. 
    space: a set of outcomes or a probability distribution of {outcome: frequency} pairs.
    if is_predicate(event):
        event = such_that(event, space)
    if isinstance(space, ProbDist):
        return sum(space[o] for o in space if o in event)
    else:
        return Fraction(len(event & space), len(space))
    
def such_that(predicate, space): 
    The outcomes in the ssample pace for which the predicate is true.
    If space is a set, return a subset {outcome,...};
    if space is a ProbDist, return a ProbDist {outcome: frequency,...};
    in both cases only with outcomes where predicate(element) is true.
    if isinstance(space, ProbDist):
        return ProbDist({o:space[o] for o in space if predicate(o)})
    else:
        return {o for o in space if predicate(o)}
Explanation: We also need to modify the functions P and such_that to accept either a sample space  or a probability distribution as the second argument.
End of explanation
DK = ProbDist(GG=121801, GB=126840,
              BG=127123, BB=135138)
DK
Explanation: Here is the probability distribution for Danish two-child families:
End of explanation
def first_girl(outcome):  return outcome[0] == 'G'
def first_boy(outcome):   return outcome[0] == 'B'
def second_girl(outcome): return outcome[1] == 'G'
def second_boy(outcome):  return outcome[1] == 'B'
def two_girls(outcome):   return outcome    == 'GG'
P(first_girl, DK)
P(second_girl, DK)
Explanation: And here are some predicates that will allow us to answer some questions:
End of explanation
P(second_girl, such_that(first_girl, DK)), P(second_girl, such_that(first_boy, DK))
P(second_boy, such_that(first_girl, DK)), P(second_boy, such_that(first_boy, DK))
Explanation: 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.
End of explanation
bag94 = ProbDist(brown=30, yellow=20, red=20, green=10, orange=10, tan=10)
bag96 = ProbDist(blue=24, green=20, orange=16, yellow=14, red=13, brown=13)
Explanation: 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 author Allen 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 and bag96:
End of explanation
def joint(A, B, sep=''):
    The joint distribution of two independent probability distributions. 
    Result is all entries of the form {a+sep+b: P(a)*P(b)}
    return ProbDist({a + sep + b: A[a] * B[b]
                    for a in A
                    for b in B})
MM = joint(bag94, bag96, ' ')
MM
Explanation: Next, define MM as the joint distribution&mdash;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.
End of explanation
def yellow_and_green(outcome): return 'yellow' in outcome and 'green' in outcome
such_that(yellow_and_green, MM)
Explanation: First we'll look at the "One is yellow and one is green" part:
End of explanation
def yellow94(outcome): return outcome.startswith('yellow')
P(yellow94, such_that(yellow_and_green, MM))
Explanation: 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?
End of explanation
from collections import Counter, deque
import random
# The board: a list of the names of the 40 squares
# As specified by https://projecteuler.net/problem=84
board = GO   A1 CC1 A2  T1 R1 B1  CH1 B2 B3
           JAIL C1 U1  C2  C3 R2 D1  CC2 D2 D3 
           FP   E1 CH2 E2  E3 R3 F1  F2  U2 F3 
           G2J  G1 G2  CC3 G3 R4 CH3 H1  T2 H2.split()
def monopoly(steps):
    Simulate given number of steps of Monopoly game, 
    yielding the number of the current square after each step.
    global here
    here = 0    # The square number where we currently are
    CC_deck = Deck('GO JAIL' + 14 * ' ?')
    CH_deck = Deck('GO JAIL C1 E3 H2 R1 R R U -3' + 6 * ' ?')
    doubles = 0
    jail = board.index('JAIL')
    for _ in range(steps):
        d1, d2 = random.randint(1, 6), random.randint(1, 6)
        goto(here + d1 + d2)
        doubles = (doubles + 1) if (d1 == d2) else 0
        if doubles == 3 or board[here] == 'G2J':  
            goto(jail)
        elif board[here].startswith('CC'):
            do_card(CC_deck)
        elif board[here].startswith('CH'):
            do_card(CH_deck)
        yield here         
def goto(square):
    "Update 'here' to be square."
    global here
    here = square % len(board)
                   
def Deck(names):
    "Make a shuffled deck of cards, given space-delimited names."
    cards = names.split()
    random.shuffle(cards)
    return deque(cards)         
def do_card(deck):
    "Take the top card from deck and do what it says."
    global here
    card = deck[0]             # The top card
    deck.rotate(-1)            # Move top card to bottom of deck
    if card == 'R' or card == 'U': 
        while not board[here].startswith(card):
            goto(here + 1)     # Advance to next railroad or utility
    elif card == '-3':
        goto(here - 3)         # Go back 3 spaces
    elif card != '?':
        goto(board.index(card))# Go to destination named on card
Explanation: So there is a 74% chance that the yellow comes from the 1994 bag.
My favorite Bayes's Theorem Problems
Answering this question is straightforward: just like all the other probability problems, we simply create a sample space, and use P to pick out the probability of the event in question, given what we know about the outcome.
<center><a href="https://en.wikipedia.org/wiki/Thomas_Bayes">Rev. Thomas Bayes</a></center>
It is curious that we were able to solve this problem with the same methodology as the others: this problem comes from a section titled My 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, we could solve it using Bayes Theorem, and if you only have pencil and paper, and not a computer, that's a great way to go. Why is Bayes Theorem recommended? Because we are asked about the probability of an event given the evidence, where that probability 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 and B, both with equal probability:
A: first M&amp;M from 94 bag, second from 96 bag
B: first M&amp;M from 96 bag, second from 94 bag
P(A) = P(B) = 0.5
Then we get some evidence:
E: first M&amp;M yellow, second green
We want to know the probability of hypothesis A, 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.04
P(E | B) = 0.10 * 0.14 
         = 0.014
P(A)     = 0.5
P(B)     = 0.5
P(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: you can use Bayes Theorem to calculate answers like this algebraically, or you can use sample spaces to calculate answers exhaustively.
There is one important question that Allen Downey does not address: would you  eat twenty-year-old M&Ms?
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 to simulate the situation, rather than one to enumerate the complete sample space. Random sampling from the simulation
can give an accurate estimate of the probability.
Simulating Monopoly
Consider problem 84 from the excellent Project 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 can 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:
End of explanation
results = list(monopoly(400000))
Explanation: And the results:
End of explanation
%matplotlib inline 
import matplotlib.pyplot as plt
plt.hist(results, bins=40)
avg = len(results) / 40
plt.plot([0, 39], [avg, avg], 'r--');
Explanation: I'll show a histogram of the squares, with a dotted red line at the average:
End of explanation
ProbDist(Counter(board[i] for i in results))
Explanation: Another way to see the results:
End of explanation
from random import gauss, triangular, choice, vonmisesvariate, uniform
def SC(): return posint(gauss(15.1, 3) + 3 * triangular(1, 4, 13)) # 30.1
def KT(): return posint(gauss(10.2, 3) + 3 * triangular(1, 3.5, 9)) # 22.1
def DG(): return posint(vonmisesvariate(30, 2) * 3.08) # 14.0
def HB(): return posint(gauss(6.7, 1.5) if choice((True, False)) else gauss(16.7, 2.5)) # 11.7
def OT(): return posint(triangular(5, 17, 25) + uniform(0, 30) + gauss(6, 3)) # 37.0
def posint(x): "Positive integer"; return max(0, int(round(x)))
Explanation: 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, and CH3, 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 an outcome 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 components random 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.)
The Central 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 a normal distribution (also called a Gaussian distribution or a bell-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:
End of explanation
from statistics import mean
def repeated_hist(rv, bins=10, k=100000):
    "Repeat rv() k times and make a histogram of the results."
    samples = [rv() for _ in range(k)]
    plt.hist(samples, bins=bins)
    return mean(samples)
Explanation: And here is a function to sample a random variable k times, show a histogram of the results, and return the mean:
End of explanation
repeated_hist(SC, bins=range(60))
repeated_hist(KT, bins=range(60))
Explanation: The two top-scoring players have scoring distributions that are slightly skewed from normal:
End of explanation
repeated_hist(DG, bins=range(60))
repeated_hist(HB, bins=range(60))
Explanation: The next two players have bi-modal distributions; some games they score a lot, some games not:
End of explanation
repeated_hist(OT, bins=range(60))
Explanation: The fifth "player" represents all the others on the team; the distribution appears to be a cross between a normal and a uniform distribution:
End of explanation
def GSW(): return SC() + KT() + DG() + HB() + OT()
repeated_hist(GSW, bins=range(70, 160, 2))
Explanation: Now we define the team score to be the sum of the players, and look at the distribution:
End of explanation
def number(cutoff):
    "Play the game with given cutoff, returning the first or second random number."
    first = random.random()
    return first if first > cutoff else random.random()
number(.5)
Explanation: 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 the Strength 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 seen how to compute probabilities. Just 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. 
<center><a href="https://en.wikipedia.org/wiki/Count_von_Count">The Count</a></center>
<hr>
Appendix: Continuous Sample Spaces
Everything up to here has been about discrete, finite sample spaces, where we can enumerate all the possible outcomes. 
ButI was asked about continuous 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 of counting 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 posed this problem in the 538 Riddler blog:
Two players go on a hot new game show called Higher 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 player A chooses a cutoff of A = 0.6, that means that A would accept any first number greater than 0.6, and reject any number below that cutoff. The question is: What cutoff, A, should player A 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 that random.random() returns a float sampled uniformly from the interval [0..1]):
End of explanation
def Pwin(A, B, trials=30000):
    "The probability that cutoff A wins against cutoff B."
    Awins = sum(number(A) > number(B) 
                for _ in range(trials))
    return Awins / trials
Pwin(.5, .6)
Explanation: Now compare the numbers returned with a cutoff of A versus a cutoff of B, and repeat for a large number of trials; this gives us an estimate of the probability that cutoff A is better than cutoff B:
End of explanation
def top(N, cutoffs):
    "Return the N best cutoffs and the number of opponent cutoffs they beat."
    winners = Counter(A if Pwin(A, B) > 0.5 else B
                      for (A, B) in itertools.combinations(cutoffs, 2))
    return winners.most_common(N)
from numpy import arange
%time top(5, arange(0.50, 0.99, 0.01))
Explanation: 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 the N top cutoffs (the ones that defeated the most number of opponent cutoffs), and the number of opponents they defeat:
End of explanation
def Phigher(A, B):
    "Probability that a sample from [A..1] is higher than one from [B..1]."
    if A <= B:
        return (1 - B) / (2 * (1 - A))
    else:
        return 1 - Phigher(B, A)
Phigher(0.5, 0.6)
Explanation: 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 making Pwin(A, B) an exact calculation. But before we get to Pwin(A, B), let's solve a simpler problem: assume that both players A and B have chosen a cutoff, and have each received a number above the cutoff. What is the probability that A gets the higher number? We'll call this Phigher(A, B). We can think of this as a two-dimensional sample space of points in the (a, b) plane, where a ranges from the cutoff A to 1 and b ranges from the cutoff B to 1. Here is a diagram of that two-dimensional sample space, with the cutoffs A=0.5 and B=0.6:
<img src="http://norvig.com/ipython/probability2da.jpg" width=413>
The total area of the sample space is 0.5 &times; 0.4 = 0.20, and in general it is (1 - A) &middot; (1 - B). What about the favorable cases, where A beats B? That corresponds to the shaded triangle below:
<img src="http://norvig.com/ipython/probability2d.jpg" width=413>
The area of a triangle is 1/2 the base times the height, or in this case, 0.4<sup>2</sup> / 2 = 0.08, and in general, (1 - B)<sup>2</sup> / 2. So in general we have:
Phigher(A, B) = favorable / total
favorable     = ((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 cutoff A &le; B; when A > B, we need to reverse things. That gives us the code:
End of explanation
def Pwin(A, B):
    "With what probability does cutoff A win against cutoff B?"
    return (  A * B         * Phigher(0, 0)  # both below cutoff
            + (1-A) * (1-B) * Phigher(A, B)  # both above
            + (1-A) * B     * Phigher(A, 0)  # A above, B below
            + A * (1-B)     * Phigher(0, B)) # A below, B above
Explanation: We're now ready to tackle the full game. There are four cases to consider, depending on whether A and B gets a first number that is above or below their cutoff choices:
| first a | first b | P(a, b) | P(A wins &vert; a, b) |  Comment |
|:-----:|:-----:| ----------- | ------------- | ------------ |
| a > A | b > B | (1 - A) &middot; (1 - B) | Phigher(A, B) | Both above cutoff; both keep first numbers |
| a < A | b < B | A &middot; B | Phigher(0, 0)  | Both below cutoff, both get new numbers from [0..1] |
| a > A | b < B | (1 - A) &middot;  B | Phigher(A, 0) | A keeps number; B gets new number from [0..1] |
| a < A | b > B | A &middot; (1 - B) | Phigher(0, B) | A gets new number from [0..1]; B keeps number |
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) &middot; (1 - B), and if this does occur, then the probability of A winning is Phigher(A, B).
We're ready to replace the old simulation-based Pwin with a new calculation-based version:
End of explanation
def test():
    assert Phigher(0.5, 0.5) == Phigher(0.7, 0.7) == Phigher(0, 0) == 0.5
    assert Pwin(0.5, 0.5) == Pwin(0.7, 0.7) == 0.5
    assert Phigher(.6, .5) == 0.6
    assert Phigher(.5, .6) == 0.4
    return 'ok'
test()
Explanation: The exact analysis code is much more complex than the simulation code, relying on a lot of algebra and geometry and case analysis, and thus more likely to have an error. Let's define a few tests to check for obvious errors:
End of explanation
top(5, arange(0.50, 0.99, 0.01))
Explanation: Let's repeat the calculation with our new, exact Pwin:
End of explanation
top(10, arange(0.500, 0.700, 0.001))
Explanation: 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:
End of explanation
top(5, arange(0.61700, 0.61900, 0.00001))
Explanation: This says 0.618 is best, better than 0.620. We can get even more accuracy:
End of explanation
import numpy as np
from mpl_toolkits.mplot3d.axes3d import Axes3D
def map2(fn, A, B):
    "Map fn to corresponding elements of 2D arrays A and B."
    return [list(map(fn, Arow, Brow))
            for (Arow, Brow) in zip(A, B)]
cutoffs = arange(0.00, 1.00, 0.02)
A, B = np.meshgrid(cutoffs, cutoffs)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.set_xlabel('A')
ax.set_ylabel('B')
ax.set_zlabel('Pwin(A, B)')
ax.plot_surface(A, B, map2(Pwin, A, B));
Explanation: So 0.61803 is best. I don't think we need more accuracy than that.
0.61803 defeats all of the other cutoffs. But in other game-theoretic problems it might be that there is no one strategy that beats all others. Sometimes A beats B and B beats C, but C beats A. To understand what we are dealing with, it is helpful to draw a 3D plot of Pwin(A, B) for values of A and B between 0 and 1:
End of explanation
cutoffs = (set(arange(0.00,    1.00,    0.01)) | 
           set(arange(0.500,   0.700,   0.001)) | 
           set(arange(0.61700, 0.61900, 0.00001)))
max([Pwin(A, B), A, B]
    for A in cutoffs for B in cutoffs)
Explanation: What does this show us? The highest win percentage for A, the peak of the surface, occurs when A is around 0.5 and B is 0 or 1. We can confirm that, finding the maximum Pwin(A, B) for many different cutoff values of A and B:
End of explanation
max(min([Pwin(A, B), A, B] for B in cutoffs)
    for A in cutoffs)
Explanation: So A could win 62.5% of the time if only B would chose a cutoff of 0. But, unfortunately for A, a rational player B is not going to do that. We can ask what happens if the game is changed so that player A has to declare a cutoff first, and then player B gets to respond with a cutoff, with full knowledge of A's choice. In other words, what cutoff should A choose to maximize Pwin(A, B), given that B is going to take that knowledge and pick a cutoff that minimizes Pwin(A, B)?
End of explanation
min(max([Pwin(A, B), A, B] for A in cutoffs)
    for B in cutoffs)
Explanation: And what if we run it the other way around, where B chooses a cutoff first, and then A responds?
End of explanation
<END_TASK>