Hands-on Tutorials
Monte Carlo Methods and Simulations explained in real-life: modeling insomnia
Monte Carlo Methods is a group of algorithms that simulate the behavior of a complex system, or probabilistic phenomena, using inferential statistics. They simulate physical processes that are typically time-consuming, or too expensive to setup and run for a large number times.
Since it is a tool to model probabilistic real-world processes, Monte Carlo Methods are widely used in areas ranging from particle Physics and Biochemistry to Engineering. So, if you can model it, you can use Monte Carlo Methods and run simulations!
Why you want to run simulations
Monte Carlo simulations are great methodology when you want to:
- Simulate processes that are time consuming, i.e., they take a lot of time to setup the right conditions or take a long time to run before you can observe meaningful results.
- Run a large number of experiments in a short time-frame. With enough computing power, you can run thousands, even millions, of experiments in a relative short time and in parallel.
- Quickly adjust and repeat experiments. Simulations can run much faster than performing the experiment by hand or waiting for the optimal conditions to observe the phenomena. So you can run the simulation a few times, quickly make adjustments, fix bugs and get the experiment back up again. For instance, in a computer simulation, changing the concentration of a molecule means updating a number in your code, while in real life it would mean you’d have to create the compound from scratch, be precise in your measurements and get all the lab equipment running again. Sometimes that’s not an easy feat!
At the end of the day, simulations help find the optimal trade-off between time to run your experiments, having faster cycles of iteration and achieving a volume of experiments that could be much difficult to manage and maintain if they were not computer simulations.
One of the earliest Monte Carlo simulations was conducted by mathematician and naturalist Buffon, in 1777. He tossed a coin 2,048 times and recorded the results, to study the distribution of the possible outcomes[1].
But Buffon’s coin toss experiment is easy to setup, run and obtain results. Other experiments may have constraints that are difficult to overcome.
For instance, in Biochemistry, Monte Carlo simulations are used to determine the parameters of metabolic processes. These processes are usually complex and simulations are used to create parallel virtual worlds, creating multiple what if scenarios. With well known parameters, they can also use simulations to create pseudo-real samples from the results of a biochemical process for further studies[2].
Not to mention the project that actually gave a name to this methodology.
Stanislaw Ulam was working on the Manhattan project when he proposed a way to use random experiments to solve deterministic problems, such as neutron diffusion in fissionable material. And because they were working on a secret project, they needed a code name. Another colleague, Nicholas Metropolis, suggested the Monte Carlo, which was where Ulam’s uncle used to borrow money to gamble.
But instead of doing everything by hand like everyone else before them, John Von Neumann helped program the ENIAC, the first programmable computer, to run those simulations.
Different methods, the same framework
There are many Monte Carlo methods, ranging different levels of complexity.
But simpler methods include the one used by Buffon, where you have Bernoulli trials and record the probability associated with success and failure, success being heads and failure tails, or vice-versa.
Another one is to calculate the expected value of a random variable.
In this case, the output of the random variable is defined by a function[3], and this function can be anything.
If your random variable represents the toss of a fair coin, that function returns heads or tails with a probability of 1/2. But you have control over the random variable, and it can be anything that makes sense in the context of your problem.
To calculate the expected value of the random variable, you repeat the following process multiple times:
- Pick a value from the set of possible values the random variable accepts,
- Plug the value into the function that generates the output of the random variable,
- Record the output.
Take the average of all outputs, and you’ve just simulated the expected value of that random variable!
Why does this work?
This only works if you run a big enough number of simulations, and so, the Law of Large Numbers can be applied.
The Law of Large Numbers states something like this:
The results obtained from performing a large number of trials should be close to the expected value. And it will become closer to the true expected value, the more trials you perform.
If you check this criteria, you can trust the result of your simulations.
Monte Carlo simulations: the framework
Regardless of the method you’re using, Monte Carlo simulations share a framework.
- Model the process or phenomenon you want to simulate.
- Define the inputs and outputs.
- Run the experiment, i.e., run the inputs through the model, and aggregate results.
1. Model the process
To run Monte Carlo Simulations you need a model, a way to represent the conditions of your experiment and how the different elements in the experiment interact (or not) to achieve a result or output.
In a simplistic way, a model transforms inputs into an output.
The model represents a probabilistic system, something that takes inputs and returns an output with a given probability. Everything depends on the model. If you don’t have a good model, you’ll have consistently poor results, Garbage in, Garbage Out.
2. Define inputs
What’s particular about Monte Carlo simulations, is that the inputs are randomly generated. Continuous inputs are generated based on a probability density function, and discrete inputs based on a probability mass function.
Continuous inputs
Probability Density Functions (PDFs) describe the probability of each output of a continuous random variable. They give you a relative likelihood that the random variable will be equal to any given value. The interesting part is that it’s the relative probability, not an absolute one.
Since we are dealing with continuous random variables, it’s hard to isolate the probability of a single value.
It’s difficult to say, for example, what is the probability that all cats in my neighborhood have a tail with exactly 12 inches?
Tail length is continuous, so we tend to look at the probabilities of continuous random variables as intervals.
You might ask instead, What’s the probability that all cats in my neighborhood have a tail with length between 10 and 12 inches?
Visually, you’d be interested in the following areas under the PDF curve.
Discrete inputs
For inputs that can be described using discrete random variables, you can only get results from a preset list of outcomes. So you use a Probability Mass Function (PMF) to get the probability of a particular outcome. It’s the equivalent of a PDF, but for discrete random variables.
When it’s your turn to play Risk and you choose to attack an enemy, you roll the dice. In 6-sided fair dice, you know you can only get a value between 1 and 6 for each die. You can’t get 2.3 on a die roll.
With discrete random variables, you can ask questions about absolute probabilities, like what’s the probability that I get a 4 in the next roll? Or, if your opponent rolls the die, and gets a 3 and it’s your turn to attack, you might wonder what’s the probability of getting a value higher than 3?
3. Run the experiment and aggregate results
To run the experiment you plug the inputs in your model, or function, and record the results.
But what do you do with the result of 100,000 simulations?
Even though you run multiple simulations, the goal is to have an aggregated output that explains the phenomenon you’ve modeled.
This aggregated output, can be the mean value for what you’ve modeled, like when you’re modeling temperatures, or the probability of a success/failure scenario, such as the likelihood of a team winning the Super Bowl.
Running these simulations is a way to get all the pieces in place to use inferential statistics. The idea is to get lots outputs from the simulations, very much like you take random samples to study a population, and then use inferential statistics to draw conclusions about the actual phenomenon.
But again, the key is to run a big number of simulations. If you run your simulation just a handful of times you might get just one-off, extreme outcomes.
Try tossing a coin 5 times in a row. Maybe you’ll get all heads, or all tails. How sure are you that the probability of getting heads/tails is 1? Based on your 5 tosses it is, but we know that if they are truly independent tosses the probability would be close to 0.5. So there’s a lot of variance in the result from those 5 tosses.
When you’re studying a population, the way to reduce variance it to increase the sample size. With Monte Carlo methods, you simply run more simulations.
Let’s put this all together!
Simulating (work from home) insomnia 🖥☕️
You’ve noticed that since starting to work from home, you’ve had a harder time falling asleep at night, and even some episodes of insomnia. Although there may be several other factors, you quickly realize that you’ve just been having too much coffee altogether.
Usually the day starts with a fresh brew that used to sit there half empty until the next day. But recently that is all gone by the end of the day! You start getting a sense that, on most days, you go over 5 cups of coffee. And for you, that’s a recipe for insomnia.
To grab a hold of the situation and control the insomnia episodes, you start thinking about what may be leading you to drink more coffee:
- Working from couch On days when you spend most of time working from the couch instead of from your desk, you’re slouching more so you feel kind of sleepy all the time. On those days, you drink a cup of coffee every couple of hours.
- Pajama pants some days you don’t feel like changing from your pajama pants. Besides making you feel more comfy, you actually feel like you haven’t really woken up. So on those days, you end up drinking an extra cup in the morning, just in case.
There are some exceptions, though. There are days when you’re so tired, you’ll fall asleep as soon as you get to bed.
With these clues about what makes you drink more coffee, you want to know how likely are you to get a good night’s sleep.
To do that, you can model your insomnia patterns with Monte Carlo methods, and run some simulations.
To prepare the inputs of the simulations you took note, during the course of a few weeks, how many times you worked from the couch, wore pajama pants and the days you felt so tired that you quickly fell asleep. Then you turned that into probabilities:
- Probability of working from the couch: 0.7,
- Probability of wearing pajama pants: 0.35,
- Probability of being so tired, you quickly fall asleep: 0.05.
You can think of each of these as the probability of success/failure event, meaning, you either wore pajama pants or not. This helps in deciding that the best distribution to model these probabilities is the Bernoulli distribution.
There are also a couple of constants in your model, which you can also tweak:
- Number of working hours in a day,
- Number of cups that you can drink without getting insomnia.
Here’s the model in Python.
import randomdef insomnia_model(working_hours, coffee_threshold, p_couch, p_pajamas, p_tired):''' The model, i.e., how the different inputs interact to generate the outputModel constants (not picked at random):param working_hours: number of working hours in a day
:param coffee_threshold: amount of coffee you need that will for sure result in insomniaModel inputs (picked at random)
:param p_couch: probability of working from the couch on any given day
:param p_pajamas: probability of wearing pajama pants on any given day
:param p_tired: probability of ending the day so tired that you'll fall asleep regardless of how many cups of coffee you have:return: The result of the model is either insomnia or sleep.'''# taking 1 sample from the bernoulli distribution with a given probability of success.
# a binomial distribution with just 1 trial is equivalent to the bernoulli distribution. that's done by setting size = 1 if bool(np.random.binomial(1, p_tired, 1)[0]):
# you're too tired today, you'll fall asleep regardless
return False # starting your day, picking the outfit and deciding if you're going to work from the couch
day_start = True
pajamas = bool(np.random.binomial(1, p_pajamas, 1)[0])
work_from_couch = bool(np.random.binomial(1, p_couch, 1)[0]) while working_hours > 0 and coffee_threshold > 0:
if pajamas and day_start:
# starting the day with an extra cup of coffee
coffee_threshold -= 2
day_start = False
if work_from_couch and working_hours % 2 == 0:
# if you're working from the couch, you get a cup of coffee every couple of hours
coffee_threshold -= 1
working_hours -= 1return coffee_threshold == 0
Then you just need to run the model multiple times, aggregate and plot the results. Here, by default, you’ll be running 10,000 simulations.
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as pltdef run_insomnia_simulations(runs=10000, working_hours=8, coffee_threshold=6, p_couch=0.7, p_pajamas=0.35, p_tired=0.05):''' Running 10,000 simulations (default) and plotting the resultsOnly highlighting here the parameters not mentioned in the model function:param runs: number of times the model will run
:return:'''
results = []
sim_count = 0 while sim_count < runs:
results.append(insomnia_model(working_hours, coffee_threshold, p_couch, p_pajamas, p_tired))
sim_count += 1 # making the output easier to interpret
results = list(map(lambda x: 'insomnia' if x else 'sleep', results))
results_pd = pd.DataFrame(data=results)
results_pd.columns = ['type'] p_insomnia = round((results_pd[results_pd['type'] == 'insomnia'].count().values[0])/float(runs), 2)
p_sleep = round((results_pd[results_pd['type'] == 'sleep'].count().values[0]) / float(runs), 2) print('Probability of having insomnia = ' + str(p_insomnia) + ' || Probability of a good night\'s sleep = ' + str(p_sleep)) # plotting the output of all model runs
fig, ax = plt.subplots(figsize=(12, 7), constrained_layout=True)
# removing to and right border
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False) # setting the color palette
sns.set_palette("Set2")
sns.histplot(data=results_pd, x='type', hue='type', legend=False, linewidth=0)
ax.xaxis.label.set_visible(False)
plt.show()
Running the simulations with the default values.
# running simulations with all the default values
run_insomnia_simulations()
With these probabilities, you’re actually getting a good night’s sleep 77% of the days!
Tweaking the simulation
This is a simulation, so you can refine your model any time you want 😀 All you need to do is change the inputs, in this case, the probabilities of the events that lead you to drink more coffee.
Maybe a more realist model is one where you’ll always be working from the couch and roughly 60% of the days you’re wearing your pajama pants.
To see that in action, you just need to change the probability of working from the couch to be equal to 1 and probability of wearing pajama pants to 0.6.
run_insomnia_simulations(p_couch=1, p_pajamas=0.6)
With this model, you’re likely to have insomnia. In fact, 57% of the days you’re going to have insomnia if you always work from the couch and wear pajama pants 3 out of 5 days of the week.
Creating this model opens up a lot of interesting investigations! For instance, you can continue to play around with the inputs and try to find the sweet-spot where you have as few insomnia episodes as possible!
Hope you enjoyed this article, and got a better idea about how Monte Carlo methods and simulations can help you understand different phenomena.
Thanks for reading!
References
[1] Sawilowsky, Shlomo S. (2003). “You think you’ve got trivials?”. Journal of Modern Applied Statistical Methods. 2 (1): 218–225.
[2] Tenekedjiev KI, Nikolova ND, Kolev K. Applications of Monte Carlo Simulation in Modelling of Biochemical Processes. In: Mode CJ, editor. Applications of Monte Carlo Methods in Biology, Medicine and Other Fields of Science [Internet]. Rijeka (HR): InTech; 2011 Feb 28. Chapter 4. PMID: 28045483.
[3] Claus Beisbart & John D. Norton (2012): Why Monte Carlo Simulations Are Inferences and Not Experiments, International Studies in the Philosophy of Science, 26:4, 403–422
