In January, I enrolled in a master’s program in mathematics at a moderately sized state school. Because the mathematics department here is pretty small (maybe 20 tenure track professors, and another 10 or 15 other faculty), and because a lot of departmental resources go into the “service classes,” the offerings for graduate students can sometimes be a little sparse. Because of this, I ended up signing up for a course in stochastic processes.^{[1]}
Going into it, I didn’t have high expectations. I didn’t know the professor from my time spent here as an undergraduate, and, while I have my bachelor’s degree in probability and statistics, I was rather hoping to move away from those topics into something a bit more abstract. I took the course because it was about the only thing offered that I could take in order to fill out my schedule.
It is a happy accident, then, that this course has ended up being one of the more interesting courses in my college career. The study of these random processes pulls a lot of theory from seemingly different branches of mathematics: over the course of the semester, I have used material from freshman level calculus, upper division probability, differential equations, and linear algebra (just to name a few).
And, now that the semester is winding down, we have been given one more tool for coming to grips with these non-deterministic systems: simulation, which draws on material that I was first introduced to in a computer science class many years ago. There are times when we just need a numerical answer. We don’t really care if it is exact, only that it is “good enough.”^{[2]} Other times, an analytical answer isn’t even possible (this happens quite a bit when differential equations get involved—differential equations are often analytically intractable). In these situations, it is often faster to simply have a computer model the process a hundred thousand times, and see what happens.
To give a fairly trivial example, suppose you have a coin. This coin has a certain probability of landing heads when it is flipped, and a certain probability of landing tails. Formally, we might say that the chance of landing heads is \(p\), and the chance of landing tails is \((1-p)\)—if we have a fair coin, then \(p=0.5\), as there is a 50% chance of tossing heads, and a 50% chance of tossing tails. Perhaps we want to know how many heads we will get if we flip the coin some number of times.
It turns out that there is an analytical answer to that question. There is a probability distribution called the Binomial Distribution which models exactly this situation. If we recall that fact, or if we are willing to work out the probability from first principles, we might find out that the probability of tossing exactly n heads is \({N\choose n}p^n(1-p)^{N-n}\), where N is the number of times we toss the coin. Once we figure that out, we have yet more calculations to do in order to find out the number of heads that we expect to get—we have to calculate the number above for every possible number of heads, then average those numbers over how likely the outcomes are. After some more work, we might discover that we should, on average, expect to see Np heads after N tosses.^{[3]}
That general formula is nice an all, but maybe we don’t have the mathematical chops to work it out, or maybe we just don’t want to. We have powerful computers now, and we could run a huge number of simulations in the time it takes to work through the mathematics. So here is another approach: pretty much any computer is capable of generating random^{[4]} numbers. We can use those random numbers to have a computer flip coins over and over again, and track the results to determine the final outcome.
In the situation above, let’s say we are going to flip a weighted coin 5 times, and we want to know how many times it comes up heads. The coin is weighted such that it will come up heads one quarter of the time, and tails the remaining 75% of the time.
Our first task is to determine how to tell the computer when a head has been tossed, and when a tail has been tossed. Using a spreadsheet program (such as Excel), we can generate a random number between 0 and 1.^{[5]} If that number is less than 0.25, then a head has been tossed. Otherwise, a tail has been tossed. Because any number between 0 and 1 is equally likely, about 25% of the random numbers will be less than 0.25, and the remainder will be larger. Using Excel, we might put the following into a cell in order to simulate a single toss of a coin:
=IF(RAND()<0.25,1,0)
This little snippet of code generates a random number (using the RAND() function), and checks to see if that number is less than 0.25. If it is, it puts a 1 in the cell. Otherwise, it puts a 0 in the cell. We are then using 1s to represent heads, and 0s to represent tails.
Of course, we are trying to simulate several coin tosses, not just one. So we copy this same snippet of code into several cells. We want to toss a coin five times, so we can copy it into five cells, which might look like this:
A | B | C | D | E | F | G | |
1 | 0 | 1 | 1 | 0 | 0 |
Then in cell F1, we enter the code
=SUM(A1:E1)
which adds up the cells A1 through E1, and tells us the number of heads that have been tossed. At this point, we have simulated one experiment. That is, we have used the computer to flip five coins and count the number of heads.
Now comes the fun part. If we were actually flipping coins, it might take all day to run the experiment 1,000 times. However, modern computers are pretty fast, and can be used to run the experiment many, many times very quickly. All we have to do is select the row above and the 999 rows below, and use the “Fill Down” command to fill all of those rows with coin tossing experiments. It takes less than a second for Excel to compute all of these random values. The first several lines might look like this:
A | B | C | D | E | F | G | |
1 | 0 | 0 | 0 | 0 | 0 | 0 | |
2 | 1 | 0 | 0 | 1 | 1 | 3 | |
3 | 0 | 0 | 0 | 0 | 0 | 0 | |
4 | 1 | 1 | 0 | 0 | 0 | 2 | |
5 | 0 | 0 | 0 | 1 | 0 | 1 |
All of these data are nice, but they are not overwhelmingly helpful without some analysis. We are interested in the average number of heads that we expect to see, so we are going to have to compute an average. Fortunately, Excel makes this really easy: into cell G1, we can enter
=AVERAGE(F:F)
This function computes the average of all of the entries in column F. And of course, if the AVERAGE() function is too arcane, we can compute the average in a more traditional way using the command
=SUM(F:F)/1000
which adds up everything in column F then divides by the number of experiments that we ran in order to get an average. The final result might look something like this:
A | B | C | D | E | F | G | |
1 | 0 | 0 | 0 | 0 | 0 | 0 | 1.242 |
2 | 1 | 0 | 0 | 1 | 1 | 3 | |
3 | 0 | 0 | 0 | 0 | 0 | 0 | |
4 | 1 | 1 | 0 | 0 | 0 | 2 | |
5 | 0 | 0 | 0 | 1 | 0 | 1 |
This tells us that out of a thousand runs of the experiment, we averaged 1.242 heads per experiment. In other words, when we flip our weighted coin 5 times, we expect to see an average of between one and two heads. Using the formula from above, we would expect 1.25 tosses out of 5 to be heads, so our simulation comes pretty close.
For those that are curious to play around with this kind of experiment, I have created an Excel spreadsheet which can be used to play around. Of note, I’ve set it up so that the probability of tossing a head can be changed on the fly.
In this particular example, none of the math to figure out the analytical formula is terribly difficult, and any moderately advanced college math student could probably figure it out without resorting to simulation. However, even in this case, simulation gives a very nice result very quickly—it took me maybe 30 seconds to set it up and get an answer. Moreover, as noted above, there are many situations in the real world where analysis simply fails to produce anything useful. For those situations, simulation can be extraordinarily helpful.
Notes:
- [1]
- “Stochastic” is the fancy mathematics term that basically means “random.” A stochastic process is a process that unfolds through time in a way that is random. When studying stochastic processes, the goal is generally to use what we know of probability and statistics to predict how the system will behave, to within certain probabilistic limits. ↩
- [2]
- Where the exact meaning of “good enough” depends upon the particular application under consideration. ↩
- [3]
- As an aside, this should make some kind of intuitive sense. If we toss a fair coin (\(p=0.5\)) a total of 100 times (\(N=100\)), then we would expect to see 50 heads (\(Np = (100)(0.5) = 50\)). ↩
- [4]
- Well, pseudo-random, but the distinction is fairly trivial in this case. ↩
- [5]
- Of course, spreadsheets are not the only tools that we can use. There are mathematical packages like Maple or Mathematica for which simulations can be written, and if we are really into hard-core computer programming, we might write some code in a programming language like python or C++. ↩
“Stochastic happens!”
Are we friends again?
I suppose. :P