# Eight to Late

Sensemaking and Analytics for Organizations

## Monte Carlo Simulation of Projects – an (even simpler) explainer

In this article I’ll explain how Monte Carlo simulation works using an example of a project that consists of two tasks that must be carried out sequentially as shown in the figure:

Task 1 takes 3 to 7 days

Task 2 takes 2 to 5 days

The two tasks do not have any dependencies other than that they need to be completed in sequence.

(Note: in case you’re wondering about “even simpler” bit in the title – the current piece is, I think, even easier to follow than this one I wrote up some years ago).

Assume the project has been carried a number of times in the past – say 20 times – and we have the data shown below for the two tasks. For each task, we have the frequency of completion by day. So, Task 1 was completed twice on day 3 , four times on day 4 and so on. Similarly, Task 2 was completed twice on the 2nd day after the task started and 10 times the 3rd day after the task started and so on.

Consider Task 1. Since it was completed 2 times on day 3 and 4 times on day 4, it is reasonable to assume that it twice as likely that it will finish on day 4 than on day 3. In other words, the number of times a task is completed on a particular day is proportional to the probability of finishing on that day.

One can therefore approximate the probability of finishing on a particular day by dividing the number of completions on that day by the total number of times the task was performed. So, for example, the probability of finishing task 1 on day 3 is 2/20 or 0.1 and the probability of finishing it on day 4 is 0.2.

It is straightforward to calculate the probability for each of the completion days. The tables displayed below show the calculated probabilities. The tables also show the cumulative probability – this is sum of all probabilities of completion prior to (and including) current completion day. This gives the probability of finishing by the particular day – that is, on that day or any day before it.  This, rather than the probability, is typically what you want to know.

The cumulative probability has two useful properties

1. It is an increasing function (that is, it increases as the completion day increases)
2. It lies between 0 and 1

What this means is that if we pick any number between 0 and 1, we will be able to find the “completion day” corresponding to that number. Let’s try this for task one:

Say we pick 0.35. Since 0.35 lies between 0.3 and 0.75, it corresponds to a completion between day 4 and day 5. That is, the task will be completed by day 5. Indeed, any number picked between 0.3 and 0.75 will correspond to a completion by day 5.

Say we pick 0.79. Since 0.79 lies between 0.75 and 0.95, it corresponds to a completion between day 5 and day 6. That is, the task will be completed by day 6.

….and so on.  It is easy to see that any random number between 0 and 1 corresponds to a specific completion day depending on which cumulative probability interval it lies in.

Let’s pick a thousand random numbers between 0 and 1 and find the corresponding completion days for each.  It should be clear from what I have said so far that these correspond to 1000 simulations of task 1, consistent with the historical data that we have on the task.

We will do the simulations in Excel. You may want to download the workbook that accompanies this post and follow along.

Enter the completion days and the cumulative probabilities corresponding to them in rows 1 through 8 of columns A and B  as shown below.

Then enter the Excel RAND() function in cell A10 as shown in the figure below. This generates a random number between 0 and 1 (note that the random number you generate will be different from mine).

Next, fill down to cell A1009 to generate 1000 random numbers between 0 and 1 – see figure below ( again your random numbers will be different from mine)

Now in cell B10, input the formula shown below:

This nested IF() function checks which cumulative probability interval the random number lies in and returns the corresponding completion day. This is the completed by day corresponding to the inputted probability.

Fill this down to cell B1009. Your first few rows will look something like shown in the figure below:

You have now simulated Task 1 thousand times.

Next, enter the data for task 2 in columns D and E (from rows 1 through 7) and follow a similar procedure to simulate Task 2 thousand times. When you’re done, you will have something like what’s shown below (again, your random numbers and hence your completed by days will differ from mine):

Each line from row 10 to 1009 corresponds to a simulation of the project. So, this is equivalent to running the project 1000 times.

We can get completion times for each simulation by summing columns B and E, which will give us 1000 project completion times. Let’s do this in column G.

Using the MIN() and MAX() functions over the range G10:G1009, we see  that the earliest and latest days for project completion are day 5 and day 12 respectively.

Using the simulation results, we can now get approximate cumulative probabilities for each of the possible completion days (i.e days 5 through 12).

Pause for a minute and have a think about how you would do this.

–x–

OK, so here’s how you would do it for day 5

Count the number of 5s in the range G10:G1009 using the COUNTIF() function.   To estimate the probability of completion on day 5, divide this number by the total number of simulations.

To get the cumulative probability you would need to add in the probabilities for all prior completion days.  However, since day 5 is the earliest possible completion day, there is no prior day.

Let’s do day 6

Count the number of 6s in the range G10:G1009 using the COUNTIF() function.   To estimate the probability of completion on day 6, divide this number by the total number of simulations.

To get the cumulative probability you would need to add the estimated  probability of completion for day 5 to the estimated  probability of completion for day 6.

…and so on.

The resulting table, show below, is excerpted from  columns J and K of the Excel workbook linked to above. Your numbers will differ (but hopefully by not too much) from the ones shown in the table.

Now that we have done all this work, we can make statements like:

1. It is highly unlikely that we will finish before day 7.
2. There’s an 80% chance that we will finish by day 9.
3. There’s a 95% chance we’ll finish by day 10.

…and so on.

And that’s how Monte Carlo simulations work in the context of project estimation

Before we close, a word or two about data. The method we have used here assumes that you have detailed historical completion data for the tasks. However, you probably know from experience that it is rarely the case that you have this.

What do you do then?

Well, one can develop probability distributions based on subjective probabilities. Here’s how: ask the task performer for a best guess earliest, most likely and latest completion time. Based on these, one can construct triangular probability distributions that can be used in simulations. It would take me far too long to explain the procedure here so I’ll point you to an article instead.

And that’s it for this explainer. I hope it has given you a sense for how Monte Carlo simulations work.

Written by K

January 4, 2022 at 5:29 pm

## A gentle introduction to Monte Carlo simulation for project managers

This article covers the why, what and how of Monte Carlo simulation using a canonical example from project management –  estimating the duration of a small project. Before starting, however, I’d like say a few words about the tool I’m going to use.

In keeping with the format of the tutorials on this blog, I’ve assumed very little prior knowledge about probability, let alone Monte Carlo simulation. Consequently, the article is verbose and the tone somewhat didactic.

### Introduction

Estimation is key part of a project manager’s role. The most frequent (and consequential) estimates they are asked deliver relate to time and cost.  Often these are calculated and presented as point estimates: i.e. single numbers – as in, this task will take 3 days. Or, a little better, as two-point ranges – as in, this task will take between 2 and 5 days.  Better still, many use a PERT-like approach wherein estimates are based on 3 points: best, most likely and worst case scenarios – as in, this task will take between 2 and 5 days, but it’s most likely that we’ll finish on day 3.  We’ll use three-point estimates as a starting point for Monte Carlo simulation, but first, some relevant background.

It is a truism, well borne out by experience, that it is easier to estimate small, simple tasks than large, complex ones. Indeed, this is why one of the early to-dos in a project is the construction of a work breakdown structure. However, a problem arises when one combines the estimates for individual elements into an overall estimate for a project or a phase thereof. It is that a straightforward addition of individual estimates or bounds will almost always lead to a grossly incorrect estimation of overall time or cost. The reason for this is simple: estimates are necessarily based on probabilities and probabilities do not combine additively. Monte Carlo simulation provides a principled and intuitive way to obtain probabilistic estimates at the level of an entire project based on estimates of the individual tasks that comprise it.

### The problem

The best way to explain Monte Carlo is through a simple worked example. So, let’s consider a 4 task project shown in Figure 1. In the project, the second task is dependent on the first, and third and fourth are dependent on the second but not on each other. The upshot of this is that the first two tasks have to be performed sequentially and the last two can be done at the same time, but can only be started after the second task is completed.

To summarise: the first two tasks must be done in series and the last two can be done in parallel.

Figure 1; A project with 4 tasks.

Figure 1 also shows the three point estimates for each task – that is the minimum, maximum and most likely completion times. For completeness I’ve listed them below:

• Task 1 – Min: 2 days; Most Likely: 4 days; Max: 8 days
• Task 2 – Min: 3 days; Most Likely: 5 days; Max: 10 days
• Task 3 – Min: 3 days; Most Likely: 6 days; Max: 9 days
• Task 4 – Min: 2 days; Most Likely: 4 days; Max: 7 days

OK, so that’s the situation as it is given to us. The first step to  developing  an estimate is to formulate the problem in a way that it can be tackled using Monte Carlo simulation. This bring us to the important topic of the shape of uncertainty aka probability distributions.

### The shape of uncertainty

Consider the data for Task 1. You have been told that it most often finishes on day 4.  However, if things go well, it could take as little as 2 days; but if things go badly it could take as long as 8 days.  Therefore, your range of possible finish times (outcomes) is between 2 to 8 days.

Clearly, each of these outcomes is not equally likely.  The most likely outcome is that you will finish the task in 4 days (from what your team member has told you). Moreover, the likelihood of finishing in less than 2 days or more than 8 days is zero. If we plot the likelihood of completion against completion time, it would look something like Figure 2.

Figure 2: Likelihood of finishing on day 2, day 4 and day 8.

Figure 2 begs a couple of questions:

1. What are the relative likelihoods of completion for all intermediate times – i.e. those between 2 to 4 days and 4 to 8 days?
2. How can one quantify the likelihood of intermediate times? In other words, how can one get a numerical value of the likelihood for all times between 2 to 8 days?  Note that we know from the earlier discussion that this must be zero for any time less than 2 or greater than 8 days.

The two questions are actually related. As we shall soon see, once we know the relative likelihood of completion at all times (compared to the maximum), we can work out its numerical value.

Since we don’t know anything about intermediate times (I’m assuming there is no other historical data available), the simplest thing to do is to assume that the likelihood increases linearly (as a straight line) from 2 to 4 days and decreases in the same way from 4 to 8 days as shown in Figure 3. This gives us the well-known triangular distribution.

Jargon Buster: The term distribution is simply a fancy word for a plot of likelihood vs. time.

Figure 3: Triangular distribution fitted to points in Figure 1

Of course, this isn’t the only possibility; there are an infinite number of others. Figure 4 is another (admittedly weird) example.

Figure 4: Another distribution that fits the points in Figure 2.

Further, it is quite possible that the upper limit (8 days) is not a hard one. It may be that in exceptional cases the task could take much longer (for example, if your team member calls in sick for two weeks) or even not be completed at all (for example, if she then leaves for that mythical greener pasture).  Catering for the latter possibility, the shape of the likelihood might resemble Figure 5.

Figure 5: A distribution that allows for a very long (potentially) infinite completion time

The main takeaway from the above is that uncertainties should be expressed as shapes rather than numbers, a notion popularised by Sam Savage in his book, The Flaw of Averages.

[Aside:  you may have noticed that all the distributions shown above are skewed to the right – that  is they have a long tail. This is a general feature of distributions that describe time (or cost) of project tasks. It would take me too far afield to discuss why this is so, but if you’re interested you may want to check out my post on the inherent uncertainty of project task estimates.

### From likelihood to probability

Thus far, I have used the word “likelihood” without bothering to define it.  It’s time to make the notion more precise.  I’ll begin by asking the question: what common sense properties do we expect a quantitative measure of likelihood to have?

Consider the following:

1. If an event is impossible, its likelihood should be zero.
2. The sum of likelihoods of all possible events should equal complete certainty. That is, it should be a constant. As this constant can be anything, let us define it to be 1.

In terms of the example above, if we denote time by $t$ and the likelihood by $P(t)$  then:

$P(t) = 0$ for $t< 2$ and  $t> 8$

And

$\sum_{t}P(t) = 1$ where $2\leq t< 8$

Where $\sum_{t}$ denotes the sum of all non-zero likelihoods – i.e. those that lie between 2 and 8 days. In simple terms this is the area enclosed by the likelihood curves and the x axis in figures 2 to 5.  (Technical Note:  Since $t$ is a continuous variable, this should be denoted by an integral rather than a simple sum, but this is a technicality that need not concern us here)

$P(t)$ is , in fact, what mathematicians call probability– which explains why I have used the symbol $P$ rather than $L$. Now that I’ve explained what it  is, I’ll use the word “probability” instead of ” likelihood” in the remainder of this article.

With these assumptions in hand, we can now obtain numerical values for the probability of completion for all times between 2 and 8 days. This can be figured out by noting that the area under the probability curve (the triangle in figure 3 and the weird shape in figure 4) must equal 1, and we’ll do this next.  Indeed, for the problem at hand, we’ll assume that all four task durations can be fitted to triangular distributions. This is primarily to keep things  simple. However, I should emphasise that you can use any shape so long as you can express it mathematically, and I’ll say more about this towards the end of this article.

### The triangular distribution

Let’s look at the estimate for Task 1. We have three numbers corresponding to a minimummost likely and maximum time.  To keep the discussion general, we’ll call these $t_{min}$, $t_{ml}$ and $t_{max}$ respectively, (we’ll get back to our estimator’s specific numbers later).

Now, what about the probabilities associated with each of these times?

Since $t_{min}$ and $t_{max}$ correspond to the minimum and maximum times,  the probability associated with these is zero. Why?  Because if it wasn’t zero, then there would be a non-zero probability of completion for a time less than $t_{min}$ or greater than $t_{max}$ – which isn’t possible [Note: this is a consequence of the assumption that the probability varies continuously –  so if it takes on non-zero value, $p_{0}$,  at $t_{min}$ then it must take on a value slightly less than $p_{0}$ – but greater than 0 –  at $t$ slightly smaller than $t_{min}$ ] .   As far as  the most likely time,  $t_{ml}$,  is concerned:  by definition, the probability attains its highest value at time $t_{ml}$.    So, assuming the probability can be described by a triangular function, the distribution must have the form shown in Figure 6 below.

Figure 6: Triangular distribution redux.

For the simulation, we need to know the equation describing the above distribution.  Although Wikipedia will tell us the answer in a mouse-click, it is instructive to figure it out for ourselves. First, note that the area under the triangle must be equal to  1 because the task must finish at some time between $t_{min}$ and $t_{max}$.   As a consequence we have:

$\frac{1}{2}\times{base}\times{altitude}=\frac{1}{2}\times{(t_{max}-t_{min})}\times{p(t_{ml})}=1\ldots\ldots{(1)}$

where $p(t_{ml})$ is the probability corresponding to time $t_{ml}$.  With a bit of rearranging we get,

$p(t_{ml})=\frac{2}{(t_{max}-t_{min})}\ldots\ldots(2)$

To derive the probability for any time $t$ lying between $t_{min}$ and $t_{ml}$, we note that:

$\frac{(t-t_{min})}{p(t)}=\frac{(t_{ml}-t_{min})}{p(t_{ml})}\ldots\ldots(3)$

This is a consequence of the fact that the ratios on either side of equation (3)  are  equal to the slope of the line joining the points $(t_{min},0)$ and $(t_{ml}, p(t_{ml}))$.

Figure 7

Substituting (2) in (3) and simplifying a bit, we obtain:

$p(t)=\frac{2(t-t_{min})}{(t_{ml}-t_{min})(t_{max}-t_{min})}\dots\ldots(4)$ for $t_{min}\leq t \leq t_{ml}$

In a similar fashion one can show that the probability for times lying between $t_{ml}$ and $t_{max}$ is given by:

$p(t)=\frac{2(t_{max}-t)}{(t_{max}-t_{ml})(t_{max}-t_{min})}\dots\ldots(5)$ for $t_{ml}\leq t \leq t_{max}$

Equations 4 and 5 together describe the probability distribution function (or PDF)  for all times between $t_{min}$ and $t_{max}$.

As it turns out, in Monte Carlo simulations, we don’t directly work with the probability distribution function. Instead we work with the cumulative distribution function (or CDF) which is the probability, $P$,  that the task is completed by time $t$. To reiterate, the PDF, $p(t)$, is the probability of the task finishing at time $t$ whereas the CDF, $P(t)$, is the probability of the task completing by time $t$. The CDF, $P(t)$,  is essentially a sum of all probabilities between $t_{min}$ and $t$. For $t < t_{min}$ this is the area under the triangle with apexes at   ($t_{min}$, 0), (t, 0) and (t, p(t)).  Using the formula for the area of a triangle (1/2 base times height) and equation (4) we get:

$P(t)=\frac{(t-t_{min})^2}{(t_{ml}-t_{min})(t_{max}-t_{min})}\ldots\ldots(6)$ for $t_{min}\leq t \leq t_{ml}$

Noting that for $t \geq t_{ml}$, the area under the curve equals the total area minus the area enclosed by the triangle with base between t and $t_{max}$, we have:

$P(t)=1- \frac{(t_{max}-t)^2}{(t_{max}-t_{ml})(t_{max}-t_{min})}\ldots\ldots(7)$ for $t_{ml}\leq t \leq t_{max}$

As expected,  $P(t)$  starts out with a value 0 at $t_{min}$ and then increases monotonically, attaining a value of 1 at $t_{max}$.

To end this section let’s plug in the numbers quoted by our estimator at the start of this section: $t_{min}=2$, $t_{ml}=4$ and $t_{max}=8$.  The resulting PDF and CDF are shown in figures 8 and 9.

Figure 8: PDF for triangular distribution (tmin=2, tml=4, tmax=8)

Figure 9 – CDF for triangular distribution (tmin=2, tml=4, tmax=8)

### Monte Carlo in a minute

Now with all that conceptual work done, we can get to the main topic of this post:  Monte Carlo estimation. The basic idea behind Monte Carlo is to simulate the entire project (all 4 tasks in this case) a large number N (say 10,000) times and thus obtain N overall completion times.  In each of the N trials, we simulate each of the tasks in the project and add them up appropriately to give us an overall project completion time for the trial.  The resulting N overall completion times will all be different, ranging from the sum of the minimum completion times to the sum of the maximum completion times.  In other words, we will obtain the PDF and CDF for the overall completion time, which will enable us to answer questions such as:

• How likely is it that the project will be completed within 17 days?
• What’s the estimated time for which I can be 90% certain that the project will be completed? For brevity, I’ll call this the 90% completion time in the rest of this piece.

“OK, that sounds great”, you say, “but how exactly do we simulate a single task”?

Good question, and I was just about to get to that…

### Simulating a single task using the CDF

As we saw earlier, the CDF for the triangular has a S shape and ranges from 0 to 1 in value. It turns out that the S shape is characteristic of all CDFs, regardless of the details underlying PDF. Why? Because, the cumulative probability must lie between 0 and 1 (remember, probabilities can never exceed 1, nor can they be negative).

OK, so to simulate a task, we:

• generate a random number between 0 and 1, this corresponds to the probability that the task will finish at time t.
• find the time, t, that this corresponds to this value of probability. This is the completion time for the task for this trial.

Incidentally, this method is called inverse transform sampling.

An example might help clarify how inverse transform sampling works.  Assume that the random number generated is 0.4905. From the CDF for the first task, we see that this value of probability corresponds to a completion time of 4.503 days, which is the completion for this trial (see Figure 10). Simple!

Figure 10: Illustrating inverse transform sampling

In this case we found the time directly from the computed CDF. That’s not too convenient when you’re simulating the project 10,000 times. Instead, we need a programmable math expression that gives us the time corresponding to the probability directly. This can be obtained by solving equations (6) and (7) for $t$. Some straightforward algebra, yields the following two expressions for $t$:

$t = t_{min} + \sqrt{P(t)(t_{ml} - t_{min})(t_{max} - t_{min})} \ldots\ldots(8)$ for $t_{min}\leq t \leq t_{ml}$

And

$t = t_{max} - \sqrt{[1-P(t)](t_{max} - t_{ml})(t_{max} - t_{min})} \ldots\ldots(9)$ for $t_{ml}\leq t \leq t_{max}$

These can be easily combined in a single Excel formula using an IF function, and I’ll show you exactly how in a minute. Yes, we can now finally get down to the Excel simulation proper and you may want to download the workbook if you haven’t done so already.

### The simulation

Open up the workbook and focus on the first three columns of the first sheet to begin with. These simulate the first task in Figure 1, which also happens to be the task we have used to illustrate the construction of the triangular distribution as well as the mechanics of Monte Carlo.

Rows 2 to 4 in columns A and B list the min, most likely and max completion times while the same rows in column C list the probabilities associated with each of the times. For $t_{min}$ the probability is 0 and for $t_{max}$ it is 1.  The probability at $t_{ml}$ can be calculated using equation (6) which, for $t=t_{max}$, reduces to

$P(t_{ml}) =\frac{(t_{ml}-t_{min})}{t_{max}-t_{min}}\ldots\ldots(10)$

Rows 6 through 10005 in column A are simulated probabilities of completion for Task 1. These are obtained via the Excel RAND() function, which generates uniformly distributed random numbers lying between 0 and 1.  This gives us a list of probabilities corresponding to 10,000 independent simulations of Task 1.

The 10,000 probabilities need to be translated into completion times for the task. This is done using equations (8) or (9) depending on whether the simulated probability is less or greater than $P(t_{ml})$, which is in cell C3 (and given by Equation (10) above). The conditional statement can be coded in an Excel formula using the IF() function.

Tasks 2-4 are coded in exactly the same way, with distribution parameters in rows 2 through 4 and simulation details in rows 6 through 10005 in the columns listed below:

• Task 2 – probabilities in column D; times in column F
• Task 3 – probabilities in column H; times in column I
• Task 4 – probabilities in column K; times in column L

That’s basically it for the simulation of individual tasks. Now let’s see how to combine them.

For tasks in series (Tasks 1 and 2), we simply sum the completion times for each task to get the overall completion times for the two tasks.  This is what’s shown in rows 6 through 10005 of column G.

For tasks in parallel (Tasks 3 and 4), the overall completion time is the maximum of the completion times for the two tasks. This is computed and stored in rows 6 through 10005 of column N.

Finally, the overall project completion time for each simulation is then simply the sum of columns G and N (shown in column O)

Sheets 2 and 3 are plots of the probability and cumulative probability distributions for overall project completion times. I’ll cover these in the next section.

### Discussion – probabilities and estimates

The figure on Sheet 2 of the Excel workbook (reproduced in Figure 11 below) is the probability distribution function (PDF) of completion times. The x-axis shows the elapsed time in days and the y-axis the number of Monte Carlo trials that have a completion time that lie in the relevant time bin (of width 0.5 days). As an example, for the simulation shown in the Figure 11, there were 882 trials (out of 10,000) that had a completion time that lie between 16.25 and 16.75 days. Your numbers will vary, of course, but you should have a maximum in the 16 to 17 day range and a trial number that is reasonably close to the one I got.

Figure 11: Probability distribution of completion times (N=10,000)

I’ll say a bit more about Figure 11 in the next section. For now, let’s move on to Sheet 3 of workbook which shows the cumulative probability of completion by a particular day (Figure 12 below).  The figure shows the cumulative probability function (CDF), which is the sum of all completion times from the earliest possible completion day to the particular day.

Figure 12: Probability of completion by a particular day (N=10,000)

To reiterate a point made earlier,  the reason we work with the CDF  rather than the PDF is that we are interested in knowing the probability of completion by a particular date (e.g. it is 90% likely that we will finish by April 20th) rather than the probability of completion on a particular date (e.g. there’s a 10% chance we’ll finish on April 17th). We can now answer the two questions we posed earlier. As a reminder, they are:

• How likely is it that the project will be completed within 17 days?
• What’s the 90% likely completion time?

Both questions are easily answered by using the cumulative distribution chart on Sheet 3 (or Fig 12).  Reading the relevant numbers from the chart, I see that:

• There’s a 60% chance that the project will be completed in 17 days.
• The 90% likely completion time is 19.5 days.

How does the latter compare to the sum of the 90% likely completion times for the individual tasks? The 90% likely completion time for a given task can be calculated by solving Equation 9 for $t$, with appropriate values for the parameters $t_{min}$, $t_{max}$ and $t_{ml}$ plugged in, and $P(t)$ set to 0.9. This gives the following values for the 90% likely completion times:

• Task 1 – 6.5 days
• Task 2 – 8.1 days
• Task 3 – 7.7 days
• Task 4 – 5.8 days

Summing up the first three tasks (remember, Tasks 3 and 4 are in parallel) we get a total of 22.3 days, which is clearly an overestimation. Now, with the benefit of having gone through the simulation, it is easy to see that the sum of 90% likely completion times for individual tasks does not equal the 90% likely completion time for the sum of the relevant individual tasks – the first three tasks in this particular case. Why? Essentially because a Monte Carlo run in which the first three tasks tasks take as long as their (individual) 90% likely completion times is highly unlikely. Exercise:  use the worksheet to estimate how likely this is.

There’s much more that can be learnt from the CDF. For example, it also tells us that the greatest uncertainty in the estimate is in the 5 day period from ~14 to 19 days because that’s the region in which the probability changes most rapidly as a function of elapsed time. Of course, the exact numbers are dependent on the assumed form of the distribution. I’ll say more about this in the final section.

To close this section, I’d like to reprise a point I mentioned earlier: that uncertainty is a shape, not a number. Monte Carlo simulations make the uncertainty in estimates explicit and can help you frame your estimates in the language of probability…and using a tool like Excel can help you explain these to non-technical people like your manager.

### Closing remarks

We’ve covered a fair bit of ground: starting from general observations about how long a task takes, saw how to construct simple probability distributions and then combine these using Monte Carlo simulation.  Before I close, there are a few general points I should mention for completeness…and as warning.

First up, it should be clear that the estimates one obtains from a simulation depend critically on the form and parameters of the distribution used. The parameters are essentially an empirical matter; they should be determined using historical data. The form of the function, is another matter altogether: as pointed out in an earlier section, one cannot determine the shape of a function from a finite number of data points. Instead, one has to focus on the properties that are important. For example, is there a small but finite chance that a task can take an unreasonably long time? If so, you may want to use a lognormal distribution…but remember, you will need to find a sensible way to estimate the distribution parameters from your historical data.

Second, you may have noted from the probability distribution curve (Figure 11)  that despite the skewed distributions of the individual tasks, the distribution of the overall completion time is somewhat symmetric with a minimum of ~9 days, most likely time of ~16 days and maximum of 24 days.  It turns out that this is a general property of distributions that are generated by adding a large number of independent probabilistic variables. As the number of variables increases, the overall distribution will tend to the ubiquitous Normal distribution.

The assumption of independence merits a closer look.  In the case it hand,  it implies that the completion times for each task are independent of each other. As most project managers will know from experience, this is rarely the case: in real life,  a task that is delayed will usually have knock-on effects on subsequent tasks. One can easily incorporate such dependencies in a Monte Carlo simulation. A formal way to do this is to introduce a non-zero correlation coefficient between tasks as I have done here. A simpler and more realistic approach is to introduce conditional inter-task dependencies As an example, one could have an inter-task delay that kicks in only if the predecessor task takes more than 80%  of its maximum time.

Thirdly, you may have wondered why I used 10,000 trials: why not 100, or 1000 or 20,000. This has to do with the tricky issue of convergence. In a nutshell, the estimates we obtain should not depend on the number of trials used.  Why? Because if they did, they’d be meaningless.

Operationally, convergence means that any predicted quantity based on aggregates should not vary with number of trials.  So, if our Monte Carlo simulation has converged, our prediction of 19.5 days for the 90% likely completion time should not change substantially if I increase the number of trials from ten to twenty thousand. I did this and obtained almost the same value of 19.5 days. The average and median completion times (shown in cell Q3 and Q4 of Sheet 1) also remained much the same (16.8 days). If you wish to repeat the calculation, be sure to change the formulas on all three sheets appropriately. I was lazy and hardcoded the number of trials. Sorry!

Finally, I should mention that simulations can be usefully performed at a higher level than individual tasks. In their highly-readable book,  Waltzing With Bears: Managing Risk on Software Projects, Tom De Marco and Timothy Lister show how Monte Carlo methods can be used for variables such as  velocity, time, cost etc.  at the project level as opposed to the task level. I believe it is better to perform simulations at the lowest possible level, the main reason being that it is easier, and less error-prone, to estimate individual tasks than entire projects. Nevertheless, high level simulations can be very useful if one has reliable data to base these on.

There are a few more things I could say about the usefulness of the generated distribution functions and Monte Carlo in general, but they are best relegated to a future article. This one is much too long already and I think I’ve tested your patience enough. Thanks so much for reading, I really do appreciate it and hope that you found it useful.

Acknowledgement: My thanks to Peter Holberton for pointing out a few typographical and coding errors in an earlier version of this article. These have now been fixed. I’d be grateful if readers could bring any errors they find to my attention.

Written by K

March 27, 2018 at 4:11 pm

Tagged with

## The drunkard’s dartboard revisited: yet another Excel-based example of Monte Carlo simulation

(Note: An Excel sheet showing sample calculations and plots discussed in this post can be downloaded here.)

### Introduction

Some months ago, I wrote a post explaining the basics of Monte Carlo simulation using the example of a drunkard throwing darts at a board. In that post I assumed that the darts could land anywhere on the dartboard with equal probability. In other words, the hit locations were assumed to be uniformly distributed. In a comment on the piece, George Gkotsis challenged this assumption, arguing that that regardless of the level of inebriation of the thrower, a dart would be more likely to land near the centre of the board than away from it (providing the player is at least moderately skilled). He also suggested using the Normal Distribution to model the spread of hits, with the variance of the distribution serving as a rough measure of the inaccuracy (or drunkenness!) of the drunkard. In George’s words:

I would propose to introduce a ‘skill’ factor, which represents the circle/square ratio (maybe a normal-Gaussian distribution). Of course, this skill factor would be very low (high variance) for a drunken player, but would still take into account the fact that throwing darts into a square is not purely random.

In this post I revisit the drunkard’s dartboard, taking into account George’s suggestions.

### Setting the stage

To keep things simple, I’ll make the following assumptions:

Figure 1: The dartboard

1. The dartboard is a circle of radius 0.5 units centred at the origin (see Figure 1)
2. The chance of a hit is greatest at the centre of the dartboard and falls off as one moves away from it.
3. The distribution of hits is a function of distance from the centre but does not depend on direction. In mathematical terms, for a given distance $r$ from the centre of the dartboard, the dart can land at any angle $\theta$ with equal probability, $\theta$ being the angle between the line joining the centre of the board to the dart and the x axis. See Figure 2 for graphical representations of a hit location in terms of $r$ and $\theta$. Note that that the $x$ and $y$ coordinates can be obtained using the formulas $x = r\cos\theta$ and $y= r\sin\theta$ as s shown in Figure 2.
4. Hits are distributed according to the Normal distribution with maximum at the centre of the dartboard.
5. The variance of the Normal distribution is a measure of inaccuracy/drunkenness of the drunkard: the more drunk the drunk, the greater the variation in his aim.

Figure 2: The coordinates of a hit location

These assumptions are consistent with George’s suggestions.

### The simulation

[Note to the reader: you may want to download the demo before continuing.]

The steps of a simulation run are as follows:

1. Generate a number that is normally distributed with a zero mean and a specified standard deviation. This gives the distance, $r$, of a randomly thrown dart from the centre of the board for a player with a “inaccuracy factor” represented by the standard deviation. Column A in the demo contains normally distributed random numbers with zero mean and a standard deviation of 0.2 . Note that I selected the latter number for no other reason than the results show up clearly on a fixed-axis plot shown in Figure 2.
2. Generate a uniformly distributed random number lying between 0 and $2\pi$. This represents the angle $\theta$. This is the content of column B of the demo.
3. The numbers obtained from steps 1 and 2 for completely specify the location of a hit. The location’s $x$ and $y$ coordinates can be worked out using the formulas $x = r\cos\theta$ and $y= r\sin\theta$. These are listed in columns C and D in the Excel demo.
4. Re-run steps 1 through 4 as many times as needed. Note that the demo is set up for 5000 runs. You can change this manually or, better yet, automate it. The latter is left as an exercise for you.

It is instructive to visualize the resulting hits using a scatter plot. Among other things this can tell you, at a glance, if the results make sense. For example, we would expect hits to be symmetrically distributed about the origin because the drunkard’s throws are not biased in any particular direction around the centre). A non-symmetrical distribution is thus an indication that there is an error in the calculations.

Now, any finite collection of hits is unlikely to be perfectly symmetrical because of outliers. Nevertheless, the distributions should be symmetrical on average. To test this, run the demo a few times (hit F9 with the demo open). Notice how the position of outliers and the overall shape of the distribution of points changes randomly from simulation to simulation. In all cases, however, there is a clear maximum at the centre of the dartboard with the probability of a hit falling with distance from the centre.

Figure 3: Scatter plot for standard deviation=0.2

Figure 3 shows the results of simulations for a standard deviation of 0.2. Figures 4 and 5 show the results of simulations for standard deviations of 0.1 and 0.4.

Figure 4: Scatter plot for standard deviation=0.1

Note that the plot has fixed axes- i.e. the area depicted is the 1×1 square that encloses the dartboard, regardless of the standard deviation. Consequently, for larger standard deviations (such as 0.4) many hits will be out of range and will not show up on the plot.

Figure 5: Scatter plot for standard deviation=0.4

### Closing remarks

As I have stressed in my previous posts on Monte Carlo simulation, the usefulness of a simulation depends on the choice of an appropriate distribution. If the selected distribution does not reflect reality, neither will the simulation. This is true regardless of whether one is simulating a drunkard’s wayward aim or the duration of project task. You may have noted that the assumption of normally-distributed hits has no justification whatsoever; it is just as arbitrary as my original assumption of uniformity. In fact, the hit locations of drunken dart throws is highly unlikely to be either uniform or Normal. Nevertheless, I hope that some of my readers will find the above example to be of pedagogical value.

### Acknowledgement

Thanks to George Gkotsis for his comment which got me thinking about this post.

Written by K

November 3, 2011 at 4:59 am