Computing Probabilities
With the mathematical rules from probability theory one may compute the probability that a certain event happens, say the probability that you get one black ball when drawing three balls from a hat with four black balls, six white balls, and three green balls. Unfortunately, theoretical calculations of probabilities may soon become hard or impossible if the problem is slightly changed. There is a simple "numerical way" of computing probabilities that is generally applicable to problems with uncertainty. The principal ideas of this approximate technique is explained below, followed by with three examples of increasing complexity.
8.3.1 Principles of Monte Carlo Simulation
Assume that we perform N experiments where the outcome of each experiment is random. Suppose that some event takes place M times in these N experiments. An estimate of the probability of the event is then M/N. The estimate becomes more accurate as N is increased, and the exact probability is assumed to be reached in the limit as N — to. (Note that in this limit, M — to too, so for rare events, where M may be small in a program, one must increase N such that M is sufficiently large for M/N to become a good approximation to the probability.)
Programs that run a large number of experiments and record the outcome of events are often called simulation programs5. The mathematical technique of letting the computer perform lots of experiments
5 This term is also applied for programs that solve equations arising in mathematical models in general, but it is particularly common to use the term when random numbers are used to estimate probabilities.
based on drawing random numbers is commonly called Monte Carlo simulation. This technique has proven to be extremely useful throughout science and industry in problems where there is uncertain or random behavior is involved6. For example, in finance the stock market has a random variation that must be taken into account when trying to optimize investments. In offshore engineering, environmental loads from wind, currents, and waves show random behavior. In nuclear and particle physics, random behavior is fundamental according to quantum mechanics and statistical physics. Many probabilistic problems can be calculated exactly by mathematics from probability theory, but very often Monte Carlo simulation is the only way to solve statistical problems. Chapters 8.3.2-8.3.4 applies examples to explain the essence of Monte Carlo simulation in problems with inherent uncertainty. However, also deterministic problems, such as integration of functions, can be computed by Monte Carlo simulation (see Chapter 8.5).
8.3.2 Example: Throwing Dice
What is the probability of getting at least six eyes twice when rolling four dice? The experiment is to roll four dice, and the event we are looking for appears when we get two or more dice with six eyes. A program roll_dice1.py simulating N such experiments may look like this:
import random as random_number import sys
N = int(sys.argv[1]) # no of experiments M = 0 # no of successful events for i in range(N):
six =0 # count the no of dice with a six r1 = random_number.randint(1, 6) if r1 == 6: six += 1
r2 = random_number.randint(1, 6) if r2 == 6: six += 1
r3 = random_number.randint(1, 6) if r3 == 6: six += 1
six += 1 # successful event? if six >= 2: M += 1 p = float(M)/N print 'probability:', p
Generalization. We can easily parameterize how many dice (ndice) we roll in each experiment and how many dice with six eyes we want to see
6 "As far as the laws of mathematics refer to reality, they are not certain, as far as they are certain, they do not refer to reality." -Albert Einstein, physicist, 1879-1955.
(nsix). Thereby, we get a shorter and more general code. The increased generality usually makes it easier to apply or adapt to new problems. The resulting program is found in roll_dice2.py and is listed below:
import random as random_number import sys
N = int(sys.argv[1]) # no of experiments ndice = int(sys.argv[2]) # no of dice nsix = int(sys.argv[3]) # wanted no of dice with six eyes M = 0 # no of successful events for i in range(N):
six =0 # how many dice with six eyes?
for j in range(ndice): # roll die no. j: r = random_number.randint(1, 6) if r == 6: six += 1 # successful event? if six >= nsix: M += 1 p = float(M)/N print 'probability:', p
With this program we may easily change the problem setting and ask for the probability that we get six eyes q times when we roll q dice. The theoretical probability can be calculated to be 6-4 w 0.00077, and a program performing 105 experiments estimates this probability to 0.0008. For such small probabilities the number of successful events M is small, and M/N will not be a good approximation to the probability unless M is reasonably large, which requires a very large N. The roll_dice2. py program runs quite slowly for one million experiments, so it is a good idea to try to vectorize the code to speed up the experiments. Unfortunately, this may constitute a challenge for newcomers to programming, as shown below.
Vectorization. In a vectorized version of the roll_dice2.py program, we generate a two-dimensional array of random numbers where the first dimension reflects the experiments and the second dimension reflects the trials in each experiment:
from numpy import random, sum eyes = random.randint(1, 7, (N, ndice))
The next step is to count the number of successes in each experiment. For this purpose, we must avoid explicit loops if we want the program to run fast. In the present example, we can compare all rolls with 6, resulting in an array compare (dimension as eyes) with ones for rolls with 6 and 0 otherwise. Summing up the rows in compare, we are interested in the rows where the sum is equal to or greater than nsix. The number of such rows equals the number of successful events, which we must divide by the total number of experiments to get the probability7:
nthrows_with_6 = sum(compare, axis=1) # sum over columns (index 1) nsuccesses = nthrows_with_6 >= nsix M = sum(nsuccesses) p = float(M)/N
The complete program is found in the file roll_dice2_vec.py. Getting rid of the two loops, as we obtained in the vectorized version, speeds up the probability estimation with a factor of 40. However, the vector-ization is highly non-trivial, and the technique depends on details of how we define success of an event in an experiment.
8.3.3 Example: Drawing Balls from a Hat
Suppose there are 12 balls in a hat: four black, four red, and four blue. We want to make a program that draws three balls at random from the hat. It is natural to represent the collection of balls as a list. Each list element can be an integer 1, 2, or 3, since we have three different types of balls, but it would be easier to work with the program if the balls could have a color instead of an integer number. This is easily accomplished by defining color names:
colors = 'black', 'red', 'blue' # (tuple of strings) hat = []
for color in colors: for i in range(4):
hat.append(color)
Drawing a ball at random is performed by import random as random_number color = random_number.choice(hat) print color
Drawing n balls without replacing the drawn balls requires us to remove an element from the hat when it is drawn. There are three ways to implement the procedure: (i) we perform a hat.remove(color), (ii) we draw a random index with randint from the set of legal indices in the hat list, and then we do a del hat[index] to remove the element, or (iii) we can compress the code in (ii) to hat.pop(index).
def draw_ball(hat):
color = random_number.choice(hat) hat.remove(color)
7 This code is considered advanced so don't be surprised if you dont't understand much of it. A first step toward understanding is to type in the code and write out the individual arrays for (say) N = 2. The use of numpy's sum function is essential for efficiency.
return color, hat def draw_ball(hat):
index = random_number.randint(0, len(hat)-1) color = hat[index] del hat[index] return color, hat def draw_ball(hat):
index = random_number.randint(0, len(hat)-1) color = hat.pop(index) return color, hat
# draw n balls from the hat:
color, hat = draw_ball(hat) balls.append(color) print 'Got the balls', balls
We can extend the experiment above and ask the question: What is the probability of drawing two or more black balls from a hat with 12 balls, four black, four red, and four blue? To this end, we perform N experiments, count how many times M we get two or more black balls, and estimate the probability as M/N. Each experiment consists of making the hat list, drawing a number of balls, and counting how many black balls we got. The latter task is easy with the count method in list objects: hat.count('black') counts how many elements with value 'black' we have in the list hat. A complete program for this task is listed below. The program appears in the file balls_in_hat.py.
import random as random_number def draw_ball(hat):
"""Draw a ball using list index."""
index = random_number.randint(0, len(hat)-1)
color = hat.pop(index)
return color, hat def draw_ball(hat):
"""Draw a ball using list index."""
index = random_number.randint(0, len(hat)-1)
color = hat[index]
del hat[index]
return color, hat def draw_ball(hat):
"""Draw a ball using list element.""" color = random_number.choice(hat) hat.remove(color) return color, hat def new_hat():
colors = 'black', 'red', 'blue' # (tuple of strings) hat = []
for color in colors: for i in range(4):
hat.append(color) return hat n = int(raw_input('How many balls are to be drawn? ')) N = int(raw_input('How many experiments? '))
# run experiments: M = 0 # no of successes for e in range(N): hat = new_hat()
balls = [] # the n balls we draw for i in range(n):
color, hat = draw_ball(hat) balls.append(color) if balls.count('black') >=2: # at least two black balls? M += 1
print 'Probability:', float(M)/N
Running the program with n = 5 (drawing 5 balls each time) and N = 4000 gives a probability of 0.57. Drawing only 2 balls at a time reduces the probability to about 0.09.
One can with the aid of probability theory derive theoretical expressions for such probabilities, but it is much simpler to let the computer perform a large number of experiments to estimate an approximate probability.
A class version of the code in this section is better than the code presented, because we avoid shuffling the hat variable in and out of functions. Exercise 8.17 asks you to design and implement a class Hat.
8.3.4 Example: Policies for Limiting Population Growth
China has for many years officially allowed only one child per couple. However, the success of the policy has been somewhat limited. One challenge is the current overrepresentation of males in the population (families have favored sons to live up). An alternative policy is to allow each couple to continue getting children until they get a son. We can simulate both policies and see how a population will develop under the "one child" and the "one son" policies. Since we expect to work with a large population over several generations, we aim at vectorized code at once.
Suppose we have a collection of n individuals, called parents, consisting of males and females randomly drawn such that a certain portion (male_portion) constitutes males. The parents array holds integer values, 1 for male and 2 for females. We can introduce constants, MALE=1 and FEMALE=2, to make the code easier to read. Our task is to see how the parents array develop from one generation to the next under the two policies. Let us first show how to draw the random integer array parents where there is a probability male_portion of getting the value MALE:
r = random.random(n)
parents[r < male_portion] = MALE
parents[r >= male_portion] = FEMALE
The number of potential couples is the minimum of males and females. However, only a fraction (fertility) of the couples will actually get a child. Under the perfect "one child" policy, these couples can have one child each:
males = len(parents[parents==MALE]) females = len(parents) - males couples = min(males, females)
n = int(fertility*couples) # couples that get a child
# the next generation, one child per couple: r = random.random(n)
children = zeros(n, int) children[r < male_portion] = MALE children[r >= male_portion] = FEMALE
The code for generating a new population will be needed in every generation. Therefore, it is natural to collect the last statements statements in a separate function such that we can repeat the statements when needed.
def get_children(n, male_portion, fertility): n = int(fertility*n) r = random.random(n) children = zeros(n, int) children[r < male_portion] = MALE children[r >= male_portion] = FEMALE return children
Under the "one son" policy, the families can continue getting a new child until they get the first son:
children = get_children(couples, male_portion, fertility)
# continue with getting a new child for each daughter: daughters = children[children == FEMALE]
while len(daughters) > 0:
new_children = get_children(len(daughters), male_portion, fertility) children = concatenate((children, new_children)) daughters = new_children[new_children == FEMALE]
The program birth_policy.py organizes the code segments above for the two policies into a function advance_generation, which we can call repeatedly to see the evolution of the population.
def advance_generation(parents, policy='one child', male_portion=0.5, fertility=1.0): males = len(parents[parents==MALE]) females = len(parents) - males couples = min(males, females)
children = get_children(couples, male_portion, fertility) elif policy == 'one son':
# first try at getting a child:
children = get_children(couples, male_portion, fertility)
# continue with getting a new child for each daughter:
daughters = children[children == FEMALE] while len(daughters) > 0:
new_children = get_children(len(daughters), male_portion, fertility) children = concatenate((children, new_children)) daughters = new_children[new_children == FEMALE] return children
The simulation is then a matter of repeated calls to advance_generation:
N = 1000000 # population size male_portion = 0.51 fertility =0.92
# start with a "perfect" generation of parents:
parents = get_children(N, male_portion=0.5, fertility=1.0)
parents = advance_generation(parents, 'one son', male_portion, fertility) print '/3d: /d' / (i+1, len(parents))
Under ideal conditions with unit fertility and a male_portion of 0.5, the program predicts that the "one child" policy halves the population from one generation to the next, while the "one son" policy, where we expect each couple to get one daughter and one son on average, keeps the population constant. Increasing male_portion slightly and decreasing fertility, which corresponds more to reality, will in both cases lead to a reduction of the population. You can try the program out with various values of these input parameters.
An obvious extension is to incorporate the effect that a portion of the population does not follow the policy and get c children on average. The program birth_policy.py can account for the effect, which is quite dramatic: If 1% of the population does not follow the "one son" policy and get 4 children on average, the population grows with 50% over 10 generations (male_portion and fertility kept at the ideal values 0.5 and 1, respectively).
Normally, simple models like the difference equations (5.9) and (5.12), or the differential equations (B.11) or (B.23), are used to model population growth. However, these models track the number of individuals through time with a very simple growth factor from one generation to the next. The model above tracks each individual in the population and applies rules involving random actions to each individual. Such a detailed and much more computer-time consuming model can be used to see the effect of different policies. Using the results of this detailed model, we can (sometimes) estimate growth factors for simpler models so that these mimic the overall effect on the population size. Exercise 8.24 asks you to investigate if a certain realization of the "one son" policy leads to simple exponential growth.
Post a comment