There are several tools available for simulating realizations of probabilistic objects — outcomes from a probability space, events, values of random variables or random vectors, stochastic processes — and analyzing the results. In general, the methods below can be chained, one after the other, e.g. P.sim(1000).apply(sum).tabulate().
- Simulate:
.sim()simulates outcomes of a probability model, events, or values of random variables or processes - Apply:
.apply()applies a function to each outcome or value - Tabulate:
.tabulate()creates a table summary of simulated outcomes of a probability model or simulated values of random variables - Count:
.count()and its relatives count the values which satsify some criteria. - Filter:
.filter()and its relatives create subsets of the values which satsify some criteria. - Get:
.get()returns the results of a particular realization of the simulation. - Summary of commands: Examples using
tabulate,filter, andcountto manipulate simulation results.
from symbulate import *
%matplotlib inline
Simulate
The .draw() extension simulates one realization of the simulation. Many realizations can be simulated using .sim(); the single argument is the number of repetitions of the simulation to run. (Note that draw is most useful when defining probability spaces, while sim is most useful when actually running simulations.)
Example. Simulate 100 repetitions of rolling two fair six-sided dice; each repetition involves a pair of values.
die = list(range(1, 6+1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size = 2)
roll.sim(100)
Caution: Note that every time .sim() is called new realizations are simulated. Store simulated realizations as variables in order to perform multiple operations in different lines of code on the same set of simulated values.
Example. Ten percent of all e-mail is spam. Thirty percent of spam e-mails contain the word "money", while 2% of non-spam e-mails contain the word "money". Simulate the email status (spam or not) and wording (money or not) for 1000 emails.
def spam_sim():
email_type = BoxModel(["spam", "not spam"], probs=[.1, .9]).draw()
if email_type == "spam":
has_money = BoxModel(["money", "no money"], probs=[.3, .7]).draw()
else:
has_money = BoxModel(["money", "no money"], probs=[.02, .98]).draw()
return email_type, has_money
P = ProbabilitySpace(spam_sim)
sims = P.sim(1000)
sims
die = list(range(1, 6+1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size = 2)
roll.sim(1000).apply(sum)
User defined functions can also be applied. Standard Python commands can be used to define functions.
Example. Ten cards labeled 1, 2, $\ldots$, 10 are shuffled and dealt one at a time. Find the the probability that the number on the card matches its position in the deal for at least one of the cards. (For example, a match occurs if card 3 is the third card dealt.)
First we define a probability space whose outcomes correspond to ordered deals of the cards. Each realization of the simulation is a permutation of the numbers 0, 1, $\ldots$ 9 (using Python indexing). Simulate 10000 realizations and store as sims.
n = 10
labels = list(range(n)) # remember, Python starts the index at 0, so the cards are labebeled 0, ..., 9
P = BoxModel(labels, size = n, replace = False)
sims = P.sim(10000)
sims
Next we define a function is_match which takes as an input a permutation of 0, $\ldots$, 9 and checks if there is at least once match. (Defining custom functions like this is the main way that Python code is used with Symbulate.)
def is_match(x):
for i in range(n):
if x[i] == labels[i]:
return 'At least one match'
return 'No match'
is_match((1, 2, 0, 3, 5, 6, 7, 8, 9, 4)) # for the outcome, card 3 matches position 3
Now we apply the is_match function to the simulated realizations, stored in sims.
sims.apply(is_match)
Finally, there are various ways to estimate the probability of at least one match, which will be described below. For example, the "match/no match" results could be tabulated to determine the relative frequence of a match.
sims.apply(is_match).tabulate(normalize = True)
Tabulate
The results of a simulation can be quickly tabulated using .tabulate(). Tabulate counts the number of times each particular outcome occurs among the simulated realizations. Use .tabulate(normalize=True) to find the proportion (relative frequency) of times each outcome occurs.
Example. Roll two fair four-sided. Each realization is an ordered pair of rolls. There are 16 possible ordered pairs - (1, 1), (1, 2), ..., (4, 3), (4, 4) - all equally likely.
die = list(range(1, 4+1, 1)) # this is just a list of the numbers 1 through 4
roll = BoxModel(die, size=2)
rolls = roll.sim(16000)
rolls.tabulate()
Now sum the dice, and approximate the probability distribution of the sum using tabulate with the normalize option.
rolls.apply(sum).tabulate(normalize=True)
Individual entries of the table can be referenced using .tabulate()[label] where label represents the value of interest.
rolls.tabulate()[(2,4)]
roll_sum = rolls.apply(sum).tabulate(normalize=True)
roll_sum[6] + roll_sum[7] + roll_sum[8]
By default, tabulate only tabulates those outcomes which are among the simulated values, rather than all possible outcomes. An argument can be passed to .tabulate() to tabulate all outcomes in a given list.
die = list(range(1, 4+1, 1)) # this is just a list of the number 1 through 4
rolls = BoxModel(die).sim(2)
rolls.tabulate(die)
# Compare with
rolls.tabulate()
By default, the outcomes in the table produced by .tabulate() are in alphanumeric order. A list can be passed to .tabulate() to achieve a specified order.
BoxModel(['a', 'b', 1, 2, 3]).sim(10).tabulate([3, 'a', 2, 'b', 1])
BoxModel(['H','T']).sim(10000).count_eq('H')
die = list(range(1, 4+1, 1)) # this is just a list of the number 1 through 4
roll = BoxModel(die, size = 2)
rolls = roll.sim(16000)
rolls.count_eq((2,4))
In addition to .count_eq(), the following count functions can be used when the values are numerical.
.count_neq()counts the values not equal to the argument.count_lt()counts the values less than the argument.count_leq()counts the values less than or equal to the argument.count_gt()counts the values greater than the argument.count_geq()counts the values greater than or equal to the argument
rolls.apply(sum).count_geq(6) / 16000
You can also count the number of outcomes which satisfy some criteria specified by a user defined function. Define a function that returns True for the outcomes you want to count, and pass the function into .count(). For example, the following code is equivalent to using .count_geq(6).
def greater_than_or_equal_to_6(x):
return x >= 6
rolls.apply(sum).count(greater_than_or_equal_to_6) / 16000
Custom functions can also be used to count based on multiple criteria. For example, the following counts the pairs of rolls in which the first roll is equal to 2 and the second roll is at most 3.
def custom_count(x): # x represents a pair of values (x[0], x[1])
if (x[0] == 2) & (x[1] <= 3):
return True
rolls.count(custom_count)
Counting can also be accomplished by creating logical (True = 1, False = 0) values according to whether an outcome satisfies some criteria and then summing.
rollsums = rolls.apply(sum)
sum(rollsums >= 6) / 16000
Since a mean (average) is the sum of values divided by the number of values, changing sum to mean in the above method returns the relative frequency directly (without having to divide by the number of values).
mean(rollsums >= 6)
Heads = BoxModel(['H','T']).sim(10000).filter_eq('H')
Heads
Using len (length) with the filter functions is one way to count the simulated occurrences of outcomes which satisfy some criteria.
len(Heads)
In addition to .filter_eq(), the following filter functions can be used when the values are numerical.
.filter_neq()subsets the values not equal to the argument.filter_lt()subsets the values less than the argument.filter_leq()subsets the values less than or equal to the argument.filter_gt()subsets the values greater than the argument.filter_geq()subsets the values greater than or equal to the argument
die = list(range(1, 1+4)) # this is just a list of the number 1 through 4
sims = BoxModel(die, size=2).sim(1000).apply(sum)
len(sims.filter_geq(6)) / 1000
You can also define your own custom filter function. Define a function that returns True for the outcomes you want to keep, and pass the function into .filter(). For example, the following code is equivalent to using .filter_geq(6).
def greater_than_or_equal_to_6(x):
return x >= 6
len(sims.filter(greater_than_or_equal_to_6)) / 1000
die = list(range(1, 4+1)) # this is just a list of the number 1 through 4
roll = BoxModel(die, size = 2)
sims = roll.sim(4)
sims
sims.get(0)
sims.get(2)
Summary of commands
The tabulate, count, and filter functions in Symbulate allow for several ways of manipulating the outcomes of a simulation. As an example, here is a recap of some of the ways Symbulate can be used to estimate the probability that the sum of two fair four-sided dice is at least 6.
die = list(range(1, 4+1, 1)) # this is just a list of the number 1 through 4
roll = BoxModel(die, size = 2)
rolls = roll.sim(10000)
rollsums = rolls.apply(sum)
rollsums.count_geq(6) / 10000
(rollsums.tabulate()[6] + rollsums.tabulate()[7] + rollsums.tabulate()[8]) / 10000
(rollsums.tabulate(normalize=True)[6] + rollsums.tabulate(normalize=True)[7] + rollsums.tabulate(normalize=True)[8])
len(rollsums.filter_geq(6)) / 10000
sum(rollsums >= 6) / 10000
mean(rollsums >= 6)
We will also see that Random variables (RV) and events play an extremely important roll in defining and simulating probabilistics objects. The following code (1) defines a probability space corresponding to two rolls of a fair four-sided die, (2) define a RV X equal to the sum of the rolls, (3) simulates realizations of X and uses a count function.
die = list(range(1, 4+1, 1)) # this is just a list of the number 1 through 4
roll = BoxModel(die, size = 2)
X = RV(roll, sum)
X.sim(10000).count_geq(6) / 10000
We can also simulate and summarize realizations of the event (X >= 6).
(X >= 6).sim(10000).tabulate(normalize=True)