There are several tools available for simulating outcomes from a probability space (or values of random variables or random 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().
.sim() simulates outcomes of a probability model, or values of random variables or process.get() returns the results of a particular realization of the simulation..apply() applies a function to each outcome or value.tabulate() creates a table summary of simulated outcomes of a probability model or simulated values of random variables.filter() and its relatives create subsets of the values which satsify some criteria..count() and its relatives count the values which satsify some criteria.tabulate, filter, and count to manipulate simulation results.from symbulate import *
%matplotlib inline
The .draw() extension simulates one outcome from a probability space. Many draws can be simulated .sim(); the argument is the number of outcomes or values to simulate.
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 values are simulated. Store simulated values as a variable in order to perform multiple operations in different lines of code on the same set of 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
The outcome of a particular repetition of the simulation can be accessed with .get(). Recall that Python starts an index at 0, so .get(0) returns the first simulated value, .get(1) the second, etc.
die = list(range(1, 6+1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size = 2)
sims = roll.sim(4)
sims
sims.get(0)
sims.get(2)
Use .apply() to apply a function to each outcome of a simulation.
Example. Roll two fair six-sided dice and compute their sum.
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.
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.)
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
def is_match(x):
for i in range(n):
if x[i] == labels[i]:
return 'At least one match'
return 'No match'
sims.apply(is_match)
The outcomes of a simulation can be quickly tabulated using .tabulate(). Tabulate counts the number of times each outcome occurs among the simulated values. Use .tabulate(normalize=True) to find the proportion (relative frequency) of times each outcome occurs.
Example. Roll two fair six-sided.
die = list(range(1, 6+1, 1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size=2)
rolls = roll.sim(10000)
rolls.tabulate()
Now sum the dice, and approximate the distribution of the sum using 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[10] + roll_sum[11] + roll_sum[12]
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, 6+1, 1)) # this is just a list of the number 1 through 6
rolls = BoxModel(die).sim(5)
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])
You can get the subset of simulations equal to a particular value using .filter_eq().
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 argumentdie = list(range(1, 1+6)) # this is just a list of the number 1 through 6
sims = BoxModel(die, size=2).sim(1000).apply(sum)
len(sims.filter_geq(10)) / 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(10).
def greater_than_or_equal_to_10(x):
return x >= 10
len(sims.filter(greater_than_or_equal_to_10)) / 1000
You can count the number of simulated outcomes equal to a particular value using count_eq().
BoxModel(['H','T']).sim(10000).count_eq('H')
die = list(range(1, 6+1, 1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size = 2)
rolls = roll.sim(10000)
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 argumentrolls.apply(sum).count_geq(10) / 10000
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 keep, and pass the function into .count(). For example, the following code is equivalent to using .count_geq(10).
def greater_than_or_equal_to_10(x):
return x >= 10
rolls.apply(sum).count(greater_than_or_equal_to_10) / 10000
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 >= 10) / 10000
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 >= 10)
The tabulate, filter, and count 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 six-sided dice is at least 10.
die = list(range(1, 6+1, 1)) # this is just a list of the number 1 through 6
roll = BoxModel(die, size = 2)
rolls = roll.sim(10000)
rollsums = rolls.apply(sum)
(rollsums.tabulate()[10] + rollsums.tabulate()[11] + rollsums.tabulate()[12]) /10000
(rollsums.tabulate(normalize=True)[10] + rollsums.tabulate(normalize=True)[11] + rollsums.tabulate(normalize=True)[12])
len(rollsums.filter_geq(10)) / 10000
rollsums.count_geq(10) / 10000
sum(rollsums >= 10) / 10000
mean(rollsums >= 10)
rollsums = RV(roll, sum)
(rollsums >= 10).sim(10000).tabulate(normalize=True)