# Eight to Late

Sensemaking and Analytics for Organizations

## A gentle introduction to text mining using R

### Preamble

This article is based on my exploration of the basic text mining capabilities of  R, the open source statistical software. It is intended primarily as a tutorial  for novices in text mining as well as R. However, unlike conventional tutorials,  I spend a good bit of time setting the context by describing the problem that led me to text mining and thence  to R. I also talk about the limitations of  the techniques I describe,  and point out directions for further exploration for those who are interested. Indeed, I’ll likely explore some of these myself in future articles.

If you have no time and /or wish to cut to the chase,  please go straight to the section entitled, Preliminaries – installing R and RStudio. If you have already installed R and have worked with it, you may want to stop reading as I  doubt there’s anything I can tell you that you don’t already know 🙂

A couple of warnings are in order before we proceed. R and the text mining options we explore below are open source software. Version differences in open source can be significant and are not always documented in a way that corporate IT types are used to. Indeed, I was tripped up by such differences in an earlier version of this article (now revised). So, just for the record, the examples below were run on version 3.2.0 of R and version 0.6-1 of the tm (text mining) package for R. A second point follows from this: as is evident from its version number, the tm package is still in the early stages of its evolution. As a result – and we will see this below – things do not always work as advertised. So assume nothing, and inspect the results in detail at every step. Be warned that I do not always do this below, as my aim is introduction rather than accuracy.

### Background and motivation

Traditional data analysis is based on the relational model in which data is stored in tables. Within tables, data is stored in rows – each row representing a  single record of an entity of interest (such as a customer or an account). The columns represent attributes of the entity. For example, the customer table might consist of columns such as name, street address, city, postcode, telephone number .  Typically these  are defined upfront, when the data model is created. It is possible to add columns after the fact, but this tends to be messy because one also has to update existing rows with information pertaining to the added attribute.

As long as one asks for information that is based  only on existing attributes – an example being,  “give me a list of customers based  in Sydney” –  a database analyst can use Structured Query Language (the defacto language of relational databases ) to  get an answer.  A problem arises, however, if one asks for information that is based on attributes that are not included in the database. An example  in the above case would be: “give me a list of  customers who have made a complaint in the last twelve months.”

As a result of the above, many data modelers will include  a “catch-all”  free text column that can be used to capture additional information in an ad-hoc way. As one might imagine, this column will often end up containing several lines, or even paragraphs of text that are near impossible to analyse with the tools available in relational databases.

(Note: for completeness I should add that most database vendors have incorporated text mining capabilities into their products. Indeed, many  of them now include R…which is another good reason to learn it.)

### My story

Over the last few months, when time permits, I’ve been doing an in-depth exploration of  the data captured by my organisation’s IT service management tool.  Such tools capture all support tickets that are logged, and track their progress until they are closed.  As it turns out,  there are a number of cases where calls are logged against categories that are too broad to be useful – the infamous catch-all category called “Unknown.”  In such cases, much of the important information is captured in a free text column, which is difficult to analyse unless one knows what one is looking for. The problem I was grappling with was to identify patterns and hence define sub-categories that would enable support staff to categorise these calls meaningfully.

One way to do this  is  to guess what the sub-categories might be…and one can sometimes make pretty good guesses if one knows the data well enough.  In general, however, guessing is a terrible strategy because one does not know what one does not know. The only sensible way to extract subcategories is to  analyse  the content of the free text column systematically. This is a classic text mining problem.

Now, I knew a bit about the theory of text mining, but had little practical experience with it. So the logical place for me to  start was to look for a suitable text mining tool. Our vendor (who shall remain unnamed) has a “Rolls-Royce” statistical tool that has a good text mining add-on. We don’t have licenses for the tool, but the vendor was willing to give us a trial license for a few months…with the understanding that this was on an intent-to-purchase basis.

I therefore started looking at open source options. While doing so, I stumbled on an interesting paper by Ingo Feinerer that describes a text mining framework for the R environment. Now, I knew about R, and was vaguely aware that it offered text mining capabilities, but I’d not looked into the details.  Anyway, I started reading the paper…and kept going until I finished.

As I read, I realised that this could be the answer to my problems. Even better, it would not require me trade in assorted limbs for a license.

I decided to give it a go.

### Preliminaries – installing R and RStudio

R can be downloaded from the R Project website. There is a Windows version available, which installed painlessly on my laptop. Commonly encountered installation issues are answered in the (very helpful)  R for Windows FAQ.

RStudio is an integrated development environment (IDE) for R. There is a commercial version of the product, but there is also a free open source version. In what follows, I’ve used the free version. Like R, RStudio installs painlessly and also detects your R installation.

RStudio has  the following panels:

• A script editor in which you can create R scripts (top left). You can also open a new script editor window by going to File > New File > RScript.
• The console where you can execute R commands/scripts (bottom left)
• Environment and history (top right)
• Files in the current working directory, installed R packages, plots and a help display screen (bottom right).

Check out this short video for a quick introduction to RStudio.

You can access help anytime (within both R and RStudio) by typing a question mark before a command. Exercise: try this by typing ?getwd() and ?setwd() in the console.

I should reiterate that the installation process for both products was seriously simple…and seriously impressive.  “Rolls-Royce” business intelligence vendors could take a lesson from  that…in addition to taking a long hard look at the ridiculous prices they charge.

There is another small step before we move on to the fun stuff.   Text mining  and certain plotting packages are not installed by default so one has to install them manually The relevant packages are:

1. tm – the text mining package (see documentation). Also check out this excellent introductory article on tm.
2. SnowballC – required for stemming (explained below).
3. ggplot2 – plotting capabilities (see documentation)
4. wordcloud – which is self-explanatory (see documentation) .

(Warning for Windows users: R is case-sensitive so Wordcloud != wordcloud)

The simplest way to install packages is to use RStudio’s built in capabilities (go to Tools > Install Packages in the menu). If you’re working on Windows 7 or 8, you might run into a permissions issue when installing packages. If you do, you might find this advice from the R for Windows FAQ helpful.

### Preliminaries – The example dataset

The data I had from our service management tool isn’t  the best dataset to learn with as it is quite messy. But then, I  have a reasonable data source in my virtual backyard:  this blog. To this end, I converted all posts I’ve written since Dec 2013 into plain text form (30 posts in all). You can download the zip file of these here .

I suggest you create a new folder called – called, say, TextMining – and unzip the files in that folder.

That done, we’re good to start…

A few things to note before we proceed:

• In what follows, I enter the commands directly in the console. However,  here’s a little RStudio tip that you may want to consider: you can enter an R command or code fragment in the script editor and then hit Ctrl-Enter  (i.e. hit the Enter key while holding down the Control key) to copy the line to the console.  This will enable you to save the script as you go along.
• In the code snippets below, the functions / commands to be typed in the R console are in blue font.   The output is in black. I will also denote references to  functions / commands in the body of the article by italicising them as in “setwd()”.  Be aware that I’ve omitted the command prompt “>” in the code snippets below!
• It is best not to cut-n-paste commands directly from the article as quotes are sometimes not rendered correctly. A text file of all the code in this article is available here.

The > prompt in the RStudio console indicates that R is ready to process commands.

To see the current working directory type in getwd() and hit return. You’ll see something like:

getwd()
[1] “C:/Users/Documents”

The exact output will of course depend on your working directory.  Note the forward slashes in the path. This is because of R’s Unix heritage (backslash is an escape character in R.). So, here’s how would change the working directory to C:\Users:

setwd(“C:/Users”)

You can now use getwd()to check that setwd() has done what it should.

getwd()
[1]”C:/Users”

I won’t say much more here about R  as I want to get on with the main business of the article.  Check out this very short introduction to R for a quick crash course.

Start RStudio and open the TextMining project you created earlier.

The next step is to load the tm package as this is not loaded by default.  This is done using the library() function like so:

library(tm)

Dependent packages are loaded automatically – in this case the dependency is on the NLP (natural language processing) package.

Next, we need to create a collection of documents (technically referred to as a Corpus) in the R environment. This basically involves loading the files created in the TextMining folder into a Corpus object. The tm package provides the Corpus() function to do this. There are several ways to  create a Corpus (check out the online help using ? as explained earlier). In a nutshell, the  Corpus() function can read from various sources including a directory. That’s the option we’ll use:

#Create Corpus
docs <- Corpus(DirSource(“C:/Users/Kailash/Documents/TextMining”))

At the risk of stating the obvious, you will need to tailor this path as appropriate.

A couple of things to note in the above. Any line that starts with a # is a comment, and the “<-“ tells R to assign the result of the command on the right hand side to the variable on the left hand side. In this case the Corpus object created is stored in a variable called docs.  One can also use the equals sign (=)  for assignment if one wants to.

Type in docs to see some information about the newly created corpus:

docs
<<VCorpus>>
Metadata: corpus specific: 0, document level (indexed): 0
Content: documents: 30

The summary() function gives more details, including a complete listing of files…but it isn’t particularly enlightening.  Instead, we’ll examine a particular document in the corpus.

#inspect a particular document
writeLines(as.character(docs[[30]]))
…output not shown…

Which prints the entire content of 30th document in the corpus to the console.

### Pre-processing

Data cleansing, though tedious, is perhaps the most important step in text analysis.   As we will see, dirty data can play havoc with the results.  Furthermore, as we will also see, data cleaning is invariably an iterative process as there are always problems that are overlooked the first time around.

The tm package offers a number of transformations that ease the tedium of cleaning data. To see the available transformations  type getTransformations() at the R prompt:

> getTransformations()
[1] “removeNumbers” “removePunctuation” “removeWords” “stemDocument” “stripWhitespace”

Most of these are self-explanatory. I’ll explain those that aren’t as we go along.

There are a few preliminary clean-up steps we need to do before we use these powerful transformations. If you inspect some documents in the corpus (and you know how to do that now), you will notice that I have some quirks in my writing. For example, I often use colons and hyphens without spaces between the words separated by them. Using the removePunctuation transform  without fixing this will cause the two words on either side of the symbols  to be combined. Clearly, we need to fix this prior to using the transformations.

To fix the above, one has to create a custom transformation. The tm package provides the ability to do this via the content_transformer function. This function takes a function as input, the input function should specify what transformation needs to be done. In this case, the input function would be one that replaces all instances of a character by spaces. As it turns out the gsub() function does just that.

Here is the R code to build the content transformer, which  we will call toSpace:

#create the toSpace content transformer
toSpace <- content_transformer(function(x, pattern) {return (gsub(pattern, ” “, x))})

Now we can use  this content transformer to eliminate colons and hypens like so:

docs <- tm_map(docs, toSpace, “-“)
docs <- tm_map(docs, toSpace, “:”)

Inspect random sections f corpus to check that the result is what you intend (use writeLines as shown earlier). To reiterate something I mentioned in the preamble, it is good practice to inspect the a subset of the corpus after each transformation.

If it all looks good, we can now apply the removePunctuation transformation. This is done as follows:

#Remove punctuation – replace punctuation marks with ” “
docs <- tm_map(docs, removePunctuation)

Inspecting the corpus reveals that several  “non-standard” punctuation marks have not been removed. These include the single curly quote marks and a space-hyphen combination. These can be removed using our custom content transformer, toSpace. Note that you might want to copy-n-paste these symbols directly from the relevant text file to ensure that they are accurately represented in toSpace.

docs <- tm_map(docs, toSpace, “’”)
docs <- tm_map(docs, toSpace, “‘”)
docs <- tm_map(docs, toSpace, ” -“)

Inspect the corpus again to ensure that the offenders have been eliminated. This is also a good time to check for any other special symbols that may need to be removed manually.

If all is well, you can move  to the next step which is  to:

• Convert the corpus to lower case
• Remove all numbers.

Since R is case sensitive, “Text” is not equal to “text” – and hence the rationale for converting to a standard case.  However, although there is a tolower transformation, it is not a part of the standard tm transformations (see the output of getTransformations() in the previous section). For this reason, we have to convert tolower into a transformation that can handle a corpus object properly. This is done with the help of our new friend, content_transformer.

Here’s the relevant code:

#Transform to lower case (need to wrap in content_transformer)
docs <- tm_map(docs,content_transformer(tolower))

Text analysts are typically not interested in numbers since these do not usually contribute to the meaning of the text. However, this may not always be so. For example, it is definitely not the case if one is interested in getting a count of the number of times a particular year appears in a corpus. This does not need to be wrapped in content_transformer as it is a standard transformation in tm.

#Strip digits (std transformation, so no need for content_transformer)
docs <- tm_map(docs, removeNumbers)

Once again, be sure to inspect the corpus before proceeding.

The next step is to remove common words  from the text. These  include words such as articles (a, an, the), conjunctions (and, or but etc.), common verbs (is), qualifiers (yet, however etc) . The tm package includes  a standard list of such stop words as they are referred to. We remove stop words using the standard removeWords transformation like so:

#remove stopwords using the standard list in tm
docs <- tm_map(docs, removeWords, stopwords(“english”))

Finally, we remove all extraneous whitespaces using the stripWhitespace transformation:

#Strip whitespace (cosmetic?)
docs <- tm_map(docs, stripWhitespace)

### Stemming

Typically a large corpus will contain  many words that have a common root – for example: offer, offered and offering.  Stemming is the process of reducing such related words to their common root, which in this case would be the word offer.

Simple stemming algorithms (such as the one in tm) are relatively crude: they work by chopping off the ends of words. This can cause problems: for example, the words mate and mating might be reduced to mat instead of mate.  That said, the overall benefit gained from stemming more than makes up for the downside of such special cases.

To see what stemming does, let’s take a look at the  last few lines  of the corpus before and after stemming.  Here’s what the last bit looks  like prior to stemming (note that this may differ for you, depending on the ordering of the corpus source files in your directory):

writeLines(as.character(docs[[30]]))
flexibility eye beholder action increase organisational flexibility say redeploying employees likely seen affected move constrains individual flexibility dual meaning characteristic many organizational platitudes excellence synergy andgovernance interesting exercise analyse platitudes expose difference espoused actual meanings sign wishing many hours platitude deconstructing fun

Now let’s stem the corpus and reinspect it.

library(SnowballC)
#Stem document
docs <- tm_map(docs,stemDocument)
writeLines(as.character(docs[[30]]))
flexibl eye behold action increas organis flexibl say redeploy employe like seen affect move constrain individu flexibl dual mean characterist mani organiz platitud excel synergi andgovern interest exercis analys platitud expos differ espous actual mean sign wish mani hour platitud deconstruct fun

A careful comparison of the two paragraphs reveals the benefits and tradeoff of this relatively crude process.

There is a more sophisticated procedure called lemmatization that takes grammatical context into account. Among other things, determining the lemma of a word requires a knowledge of its part of speech (POS) – i.e. whether it is a noun, adjective etc. There are POS taggers that automate the process of tagging terms with their parts of speech. Although POS taggers are available for R (see this one, for example), I will not go into this topic here as it would make a long post even longer.

On another important note, the output of the corpus also shows up a problem or two. First, organiz and organis are actually variants of the same stem organ. Clearly, they should be merged. Second, the word andgovern should be separated out into and and govern (this is an error in the original text).  These (and other errors of their ilk) can and should be fixed up before proceeding.  This is easily done using gsub() wrapped in content_transformer. Here is the code to  clean up these and a few other issues  that I found:

docs <- tm_map(docs, content_transformer(gsub), pattern = “organiz”, replacement = “organ”)
docs <- tm_map(docs, content_transformer(gsub), pattern = “organis”, replacement = “organ”)
docs <- tm_map(docs, content_transformer(gsub), pattern = “andgovern”, replacement = “govern”)
docs <- tm_map(docs, content_transformer(gsub), pattern = “inenterpris”, replacement = “enterpris”)
docs <- tm_map(docs, content_transformer(gsub), pattern = “team-“, replacement = “team”)

Note that I have removed the stop words and and in in the 3rd and 4th transforms above.

There are definitely other errors that need to be cleaned up, but I’ll leave these for you to detect and remove.

### The document term matrix

The next step in the process is the creation of the document term matrix  (DTM)– a matrix that lists all occurrences of words in the corpus, by document. In the DTM, the documents are represented by rows and the terms (or words) by columns.  If a word occurs in a particular document, then the matrix entry for corresponding to that row and column is 1, else it is 0 (multiple occurrences within a document are recorded – that is, if a word occurs twice in a document, it is recorded as “2” in the relevant matrix entry).

A simple example might serve to explain the structure of the TDM more clearly. Assume we have a simple corpus consisting of two documents, Doc1 and Doc2, with the following content:

Doc1: bananas are good

Doc2: bananas are yellow

The DTM for this corpus would look like:

 bananas are yellow good Doc1 1 1 1 0 Doc2 1 1 0 1

Clearly there is nothing special about rows and columns – we could just as easily transpose them. If we did so, we’d get a term document matrix (TDM) in which the terms are rows and documents columns. One can work with either a DTM or TDM. I’ll use the DTM in what follows.

There are a couple of general points worth making before we proceed. Firstly, DTMs (or TDMs) can be huge – the dimension of the matrix would be number of document  x the number of words in the corpus.  Secondly, it is clear that the large majority of words will appear only in a few documents. As a result a DTM is invariably sparse – that is, a large number of its entries are 0.

The business of creating a DTM (or TDM) in R is as simple as:

dtm <- DocumentTermMatrix(docs)

This creates a term document matrix from the corpus and stores the result in the variable dtm. One can get summary information on the matrix by typing the variable name in the console and hitting return:

dtm
<<DocumentTermMatrix (documents: 30, terms: 4209)>>
Non-/sparse entries: 14252/112018
Sparsity : 89%
Maximal term length: 48
Weighting : term frequency (tf)

This is a 30 x 4209 dimension matrix in which 89% of the rows are zero.

One can inspect the DTM, and you might want to do so for fun. However, it isn’t particularly illuminating because of the sheer volume of information that will flash up on the console. To limit the information displayed, one can inspect a small section of it like so:

inspect(dtm[1:2,1000:1005])
<<DocumentTermMatrix (documents: 2, terms: 6)>>
Non-/sparse entries: 0/12
Sparsity : 100%
Maximal term length: 8
Weighting : term frequency (tf)
Docs                               creation creativ credibl credit crimin crinkl
BeyondEntitiesAndRelationships.txt        0        0      0      0      0      0
bigdata.txt                               0        0      0      0      0      0

This command displays terms 1000 through 1005 in the first two rows of the DTM. Note that your results may differ.

### Mining the corpus

Notice that in constructing the TDM, we have converted a corpus of text into a mathematical object that can be analysed using quantitative techniques of matrix algebra.  It should be no surprise, therefore, that the TDM (or DTM) is the starting point for quantitative text analysis.

For example, to get the frequency of occurrence of each word in the corpus, we simply sum over all rows to give column sums:

freq <- colSums(as.matrix(dtm))

Here we have  first converted the TDM into a mathematical matrix using the as.matrix() function. We have then summed over all rows to give us the totals for each column (term). The result is stored in the (column matrix) variable freq.

Check that the dimension of freq equals the number of terms:

#length should be total number of terms
length(freq)
[1] 4209

Next, we sort freq in descending order of term count:

#create sort order (descending)
ord <- order(freq,decreasing=TRUE)

Then list the most and least frequently occurring terms:

#inspect most frequently occurring terms
one organ can manag work system
314 268    244  222   202   193
#inspect least frequently occurring terms
freq[tail(ord)]
yield yorkshir   youtub     zeno     zero    zulli
1        1        1        1        1        1

The  least frequent terms can be more interesting than one might think. This is  because terms that occur rarely are likely to be more descriptive of specific documents. Indeed, I can recall the posts in which I have referred to Yorkshire, Zeno’s Paradox and  Mr. Lou Zulli without having to go back to the corpus, but I’d have a hard time enumerating the posts in which I’ve used the word system.

There are at least a couple of ways to simple ways to strike a balance between frequency and specificity. One way is to use so-called  inverse document frequencies. A simpler approach is to  eliminate words that occur in a large fraction of corpus documents.   The latter addresses another issue that is evident in the above. We deal with this now.

Words like “can” and “one”  give us no information about the subject matter of the documents in which they occur. They can therefore be eliminated without loss. Indeed, they ought to have been eliminated by the stopword removal we did earlier. However, since such words occur very frequently – virtually in all documents – we can remove them by enforcing bounds when creating the DTM, like so:

dtmr <-DocumentTermMatrix(docs, control=list(wordLengths=c(4, 20),
bounds = list(global = c(3,27))))

Here we have told R to include only those words that occur in  3 to 27 documents. We have also enforced  lower and upper limit to length of the words included (between 4 and 20 characters).

Inspecting the new DTM:

dtmr
<<DocumentTermMatrix (documents: 30, terms: 1290)>>
Non-/sparse entries: 10002/28698
Sparsity : 74%
Maximal term length: 15
Weighting : term frequency (tf)

The dimension is reduced to 30 x 1290.

Let’s calculate the cumulative frequencies of words across documents and sort as before:

freqr <- colSums(as.matrix(dtmr))
#length should be total number of terms
length(freqr)
[1] 1290
#create sort order (asc)
ordr <- order(freqr,decreasing=TRUE)
#inspect most frequently occurring terms
organ manag work system project problem
268     222  202    193     184     171
#inspect least frequently occurring terms
freqr[tail(ordr)]
wait warehous welcom whiteboard wider widespread
3           3      3          3     3          3

The results make sense: the top 6 keywords are pretty good descriptors of what my blogs is about – projects, management and systems. However, not all high frequency words need be significant. What they do, is give you an idea of potential classification terms.

That done, let’s take get a list of terms that occur at least a  100 times in the entire corpus. This is easily done using the findFreqTerms() function as follows:

findFreqTerms(dtmr,lowfreq=80)
[1] “action” “approach” “base” “busi” “chang” “consult” “data” “decis” “design”
[10] “develop” “differ” “discuss” “enterpris” “exampl” “group” “howev” “import” “issu”
[19] “like” “make” “manag” “mani” “model” “often” “organ” “peopl” “point”
[28] “practic” “problem” “process” “project” “question” “said” “system” “thing” “think”
[37] “time” “understand” “view” “well” “will” “work”

Here I have asked findFreqTerms() to return all terms that occur more than 80 times in the entire corpus. Note, however, that the result is ordered alphabetically, not by frequency.

Now that we have the most frequently occurring terms in hand, we can check for correlations between some of these and other terms that occur in the corpus.  In this context, correlation is a quantitative measure of the co-occurrence of words in multiple documents.

The tm package provides the findAssocs() function to do this.  One needs to specify the DTM, the term of interest and the correlation limit. The latter is a number between 0 and 1 that serves as a lower bound for  the strength of correlation between the  search and result terms. For example, if the correlation limit is 1, findAssocs() will return only  those words that always co-occur with the search term. A correlation limit of 0.5 will return terms that have a search term co-occurrence of at least  50% and so on.

Here are the results of  running findAssocs() on some of the frequently occurring terms (system, project, organis) at a correlation of 60%.

findAssocs(dtmr,”project”,0.6)
project
inher 0.82
handl 0.68
manag 0.68
occurr 0.68
manager’ 0.66
findAssocs(dtmr,”enterpris”,0.6)
enterpris
agil       0.80
realist    0.78
increment  0.77
upfront    0.77
technolog  0.70
neither    0.69
solv       0.69
architectur 0.67
happi      0.67
movement   0.67
architect  0.66
chanc      0.65
fine       0.64
featur     0.63
findAssocs(dtmr,”system”,0.6)
system
design  0.78
subset  0.78
user    0.77
involv  0.71
specifi 0.71
function 0.70
intend  0.67
softwar 0.67
step    0.67
compos  0.66
intent  0.66
specif  0.66
depart  0.65
phone  0.63
frequent 0.62
today  0.62
pattern 0.61
cognit 0.60
wherea 0.60

An important point to note that the presence of a term in these list is not indicative of its frequency.  Rather it is a measure of the frequency with which the two (search and result term)  co-occur (or show up together) in documents across . Note also, that it is not an indicator of nearness or contiguity. Indeed, it cannot be because the document term matrix does not store any information on proximity of terms, it is simply a “bag of words.”

That said, one can already see that the correlations throw up interesting combinations – for example, project and manag, or enterpris and agil or architect/architecture, or system and design or adopt. These give one further insights into potential classifications.

As it turned out,  the very basic techniques listed above were enough for me to get a handle on the original problem that led me to text mining – the analysis of free text problem descriptions in my organisation’s service management tool.  What I did was to work my way through the top 50 terms and find their associations. These revealed a number of sets of keywords that occurred in multiple problem descriptions,  which was good enough for me to define some useful sub-categories.  These are currently being reviewed by the service management team. While they’re busy with that that, I’m looking into refining these further using techniques such as  cluster analysis and tokenization.   A simple case of the latter would be to look at two-word combinations in the text (technically referred to as bigrams). As one might imagine, the dimensionality of the DTM will quickly get out of hand as one considers larger multi-word combinations.

Anyway,  all that and more will topics have to wait for future  articles as this piece is much too long already. That said, there is one thing I absolutely must touch upon before signing off. Do stay, I think you’ll find  it interesting.

### Basic graphics

One of the really cool things about R is its graphing capability. I’ll do just a couple of simple examples to give you a flavour of its power and cool factor. There are lots of nice examples on the Web that you can try out for yourself.

Let’s first do a simple frequency histogram. I’ll use the ggplot2 package, written by Hadley Wickham to do this. Here’s the code:

wf=data.frame(term=names(freqr),occurrences=freqr)
library(ggplot2)
p <- ggplot(subset(wf, freqr>100), aes(term, occurrences))
p <- p + geom_bar(stat=”identity”)
p <- p + theme(axis.text.x=element_text(angle=45, hjust=1))
p

Figure 1 shows the result.

Fig 1: Term-occurrence histogram (freq>100)

The first line creates a data frame – a list of columns of equal length. A data frame also contains the name of the columns – in this case these are term and occurrence respectively.  We then invoke ggplot(), telling it to consider plot only those terms that occur more than 100 times.  The aes option in ggplot describes plot aesthetics – in this case, we use it to specify the x and y axis labels. The stat=”identity” option in geom_bar () ensures  that the height of each bar is proportional to the data value that is mapped to the y-axis  (i.e occurrences). The last line specifies that the x-axis labels should be at a 45 degree angle and should be horizontally justified (see what happens if you leave this out). Check out the voluminous ggplot documentation for more or better yet, this quick introduction to ggplot2 by Edwin Chen.

Finally, let’s create a wordcloud for no other reason than everyone who can seems to be doing it.  The code for this is:

#wordcloud
library(wordcloud)
#setting the same seed each time ensures consistent look across clouds
set.seed(42)
#limit words by specifying min frequency
wordcloud(names(freqr),freqr, min.freq=70)

The result is shown Figure 2.

Fig 2: Wordcloud (freq>70)

Here we first load the wordcloud package which is not loaded by default. Setting a seed number ensures that you get the same look each time (try running it without setting a seed). The arguments of the wordcloud() function are straightforward enough. Note that one can specify the maximum number of words to be included instead of the minimum frequency (as I have done above).  See the word cloud  documentation for more.

This word cloud also makes it clear that stop word removal has not done its job well, there are a number of words it has missed (also and however, for example). These can be removed by augmenting the built-in stop word list with a custom one. This is left as an exercise for the reader :-).

Finally, one can make the wordcloud more visually appealing by adding colour as follows:

wordcloud(names(freqr),freqr,min.freq=70,colors=brewer.pal(6,”Dark2″))

The result is shown Figure 3.

Fig 3: Wordcloud (freq > 70)

You may need to load the RColorBrewer package to get this to work. Check out the brewer documentation to experiment with more colouring options.

### Wrapping up

This brings me to the end of this rather long  (but I hope, comprehensible) introduction to text mining R.  It should be clear that despite the length of the article, I’ve covered only the most rudimentary basics.  Nevertheless, I hope I’ve succeeded in conveying  a sense of the possibilities in the vast and rapidly-expanding discipline of text analytics.

Note added on July 22nd, 2015:

If you liked this piece, you may want to check out the sequel – a gentle introduction to cluster analysis using R.

Note added on September 29th 2015:

Note added on December 3rd 2015:

Written by K

May 27, 2015 at 8:08 pm

## On the statistical downsides of blogging

### Introduction

The stats on the 200+ posts I’ve written since I started blogging make it pretty clear that:

1. Much of what I write does not get much attention – i.e. it is not of interest to most readers.
2. An interesting post – a rare occurrence in itself  – is invariably followed by a series of uninteresting ones.

In this post, I ignore the very real possibility that my work is inherently uninteresting and discuss how the above observations can be explained via concepts of probability.

### Base rate of uninteresting ideas

A couple of years ago I wrote a piece entitled, Trumped by Conditionality, in which I used conditional probability to show that majority of the posts on this blog will be uninteresting despite my best efforts.  My argument was based on the following observations:

1. There are many more uninteresting ideas than interesting ones.  In statistical terminology one would say that the base rate of uninteresting ideas is high.   This implies that if I write posts without filtering out bad ideas, I will write uninteresting posts far more frequently than interesting ones.
2. The base rate as described above is inapplicable in real life because I do attempt to filter out the bad ideas. However, and this is the key:  my ability to distinguish between interesting and uninteresting topics is imperfect. In other words, although I can generally identify an interesting idea correctly , there is a small (but significant) chance that I will incorrectly identify an uninteresting topic as being interesting.

Now,  since uninteresting ideas vastly outnumber interesting ones and my ability to filter out uninteresting ideas is imperfect, it follows that  the majority of the topics I choose to write about will be uninteresting.   This is essentially the first point I made in the introduction.

### Regression to the mean

The observation that good (i.e. interesting) posts are generally followed by a series of not so good ones is a consequence of a statistical phenomenon known as regression to the mean.  In everyday language this refers to the common observation that an extreme event is generally followed by a less extreme one.   This is simply a consequence of the fact that for many commonly encountered phenomena extreme events are much less likely to occur than events that are close to the average.

In the case at hand we are concerned with the quality of writing. Although writers might improve through practice, it is pretty clear that they cannot write brilliant posts every time they put fingers to keyboard. This is particularly true of bloggers and syndicated columnists who have to produce pieces according to a timetable – regardless of practice or talent, it is impossible to produce high quality pieces on a regular basis.

It is worth noting that people often incorrectly ascribe causal explanations to phenomena that can be explained by regression to the mean.  Daniel Kahneman and Amos Tversky describe the following example in their classic paper on decision-related cognitive biases:

…In a discussion of flight training, experienced instructors noted that praise for an exceptionally smooth landing is typically followed by a poorer landing on the next try, while harsh criticism after a rough landing is usually followed by an improvement on the next try. The instructors concluded that verbal rewards are detrimental to learning, while verbal punishments are beneficial, contrary to accepted psychological doctrine. This conclusion is unwarranted because of the presence of regression toward the mean. As in other cases of repeated examination, an improvement will usually follow a poor performance and a deterioration will usually follow an outstanding performance, even if the instructor does not respond to the trainee’s achievement on the first attempt

So, although I cannot avoid the disappointment that follows the high of writing a well-received post, I can take (perhaps, false) comfort in the possibility that I’m a victim of statistics.

### In closing

Finally, l would be remiss if I did not consider an explanation which, though unpleasant, may well be true: there is the distinct possibility that everything I write about is uninteresting. Needless to say, I reckon the explanations (rationalisations?) offered above are far more likely to be correct 🙂

Written by K

June 1, 2012 at 6:12 am

Posted in Probability, Statistics, Writing

Tagged with

## On the accuracy of group estimates

### Introduction

The essential idea behind group estimation is that an estimate made by a group is likely to be more accurate than one made by an individual in the group. This notion is the basis for the Delphi method and its variants. In this post, I use arguments involving probabilities to gain some insight into the conditions under which group estimates are more accurate than individual ones.

### An insight from conditional probability

Let’s begin with a simple group estimation scenario.

Assume we have two individuals of similar skill who have been asked to provide independent estimates of some quantity, say  a project task duration. Further, let us assume that each individual has a probability $p$ of making a correct estimate.

Based on the above, the probability that they both make a correct estimate, $P(\textnormal{both correct})$,  is:

$P(\textnormal{both correct}) = p*p = p^2$,

This is a consequence of our assumption that the individual estimates are independent of each other.

Similarly,  the probability that they both get it wrong, $P(\textnormal{both wrong})$, is:

$P(\textnormal{both wrong}) = (1-p)*(1-p) = (1-p)^2$,

Now we can ask the following question:

What is the probability that both individuals make the correct estimate if we know that they have both made the same estimate?

This can be figured out using Bayes’ Theorem, which in the context of the question can be stated as follows:

$P(\textnormal{both correct\textbar same estimate})= \displaystyle{\frac{ P(\textnormal{same estimate\textbar both correct})*P(\textnormal{both correct})}{ P(\textnormal{same estimate})}}$

In the above equation, $P(\textnormal{both correct\textbar same estimate})$ is the probability that both individuals get it right given that they have made the same estimate (which is  what we want to figure out). This is an example of a conditional probability – i.e.  the probability that an event occurs given that another, possibly related event has already occurred.  See this post for a detailed discussion of conditional probabilities.

Similarly, $P(\textnormal{same estimate\textbar both correct})$ is the conditional probability that both estimators make the same estimate given that they are both correct. This probability is 1.

Question: Why?

Answer: If both estimators are correct then they must have made the same estimate (i.e. they must both within be an acceptable range of the right answer).

Finally, $P(\textnormal{same estimate})$ is the probability that both make the same estimate. This is simply the sum of the probabilities that both get it right and both get it wrong. Expressed in terms of $p$ this is, $p^2+(1-p)^2$.

Now lets apply Bayes’ theorem to the following two cases:

1. Both individuals are good estimators – i.e. they have a high probability of making a correct estimate. We’ll assume they both have a 90% chance of getting it right ($p=0.9$).
2. Both individuals are poor estimators – i.e. they have a low probability of making a correct estimate. We’ll assume they both have a 30% chance of getting it right ($p=0.3$)

Consider the first case. The probability that both estimators get it right given that they make the same estimate is:

$P(\textnormal{both correct\textbar same estimate})= \displaystyle\frac{1*0.9*0.9}{0.9*0.9+0.1*0.1}= \displaystyle \frac{0.81}{0.82}= 0.9878$

Thus we see that the group estimate has a significantly better chance of being right than the individual ones:  a probability of 0.9878 as opposed to 0.9.

In the second case, the probability that both get it right is:

$P(\textnormal{both correct\textbar same estimate})= \displaystyle \frac{1*0.3*0.3}{0.3*0.3+0.7*0.7}= \displaystyle \frac{0.09}{0.58}= 0.155$

The situation is completely reversed: the group estimate has a much smaller chance of being right than an  individual estimate!

In summary:  estimates provided by a group consisting of individuals of similar ability working independently are more likely to be right (compared to individual estimates) if the group consists of  competent estimators and more likely to be wrong (compared to individual estimates) if the group consists of  poor estimators.

### Assumptions and complications

I have made a number of simplifying assumptions in the above argument. I discuss these below with some commentary.

1. The main assumption is that individuals work independently. This assumption is not valid for many situations. For example, project estimates are often made  by a group of people working together.  Although one can’t work out what will happen in such situations using the arguments of the previous section, it is reasonable to assume that given the right conditions, estimators will use their collective knowledge to work collaboratively.   Other things being equal,  such collaboration would lead a group of skilled estimators to reinforce each others’ estimates (which are likely to be quite similar) whereas less skilled ones may spend time arguing over their (possibly different and incorrect) guesses.  Based on this, it seems reasonable to conjecture that groups consisting of good estimators will tend to make even better estimates than they would individually whereas those consisting of poor estimators have a significant chance of making worse ones.
2. Another assumption is that an estimate is either good or bad. In reality there is a range that is neither good nor bad, but may be acceptable.
3. Yet another assumption is that an estimator’s ability can be accurately quantified using a single numerical probability.  This is fine providing the number actually represents the person’s estimation ability for the situation at hand. However, typically such probabilities are evaluated on the basis of  past estimates. The problem is, every situation is unique and history may not be a good guide to the situation at hand. The best way to address this is to involve people with diverse experience in the estimation exercise.  This will almost often lead to  a significant spread of estimates which may then have to be refined by debate and negotiation.

Real-life estimation situations have a number of other complications.  To begin with, the influence that specific individuals have on the estimation process may vary – a manager who is  a poor estimator may, by virtue of his position, have a greater influence than others in a group. This will skew the group estimate by a factor that cannot be estimated.  Moreover, strategic behaviour may influence estimates in a myriad other ways. Then there is the groupthink factor  as well.

…and I’m sure there are many others.

Finally I should mention that group estimates can depend on the details of the estimation process. For example, research suggests that under certain conditions competition can lead to better estimates than cooperation.

### Conclusion

In this post I have attempted to make some general inferences regarding the validity of group estimates based on arguments involving conditional probabilities. The arguments suggest that, all other things being equal, a collective estimate from a bunch of skilled estimators will generally be better than their individual estimates whereas an estimate from a group of less skilled estimators will tend to be worse than their individual estimates. Of course, in real life, there are a host of other factors  that can come into play:  power, politics and biases being just a few. Though these are often hidden, they can  influence group estimates in inestimable ways.

### Acknowledgement

Thanks go out to George Gkotsis and Craig Brown for their comments which inspired this post.

Written by K

December 1, 2011 at 5:16 am

## 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

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

Written by K

November 3, 2011 at 4:59 am

### Introduction

More often than not, managerial decisions are made on the basis of uncertain information. To lend some rigour to the process of decision making, it is sometimes assumed that uncertainties of interest can be quantified accurately using probabilities. As it turns out, this assumption can be incorrect in many situations because the probabilities themselves can be uncertain.   In this post I discuss a couple of ways in which such uncertainty about uncertainty can manifest itself.

### The problem of vagueness

In a paper entitled, “Is Probability the Only Coherent Approach to Uncertainty?”,  Mark Colyvan made a distinction between two types of uncertainty:

1. Uncertainty about some underlying fact. For example, we might be uncertain about the cost of a project – that there will be a cost is a fact, but we are uncertain about what exactly it will be.
2. Uncertainty about situations where there is no underlying fact.  For example, we might be uncertain about whether customers will be satisfied with the outcome of a project. The problem here is the definition of customer satisfaction. How do we measure it? What about customers who are neither satisfied nor dissatisfied?  There is no clear-cut definition of what customer satisfaction actually is.

The first type of uncertainty refers to the lack of knowledge about something that we know exists. This is sometimes referred to as epistemic uncertainty – i.e. uncertainty pertaining to knowledge. Such uncertainty arises from imprecise measurements, changes in the object of interest etc.  The key point is that we know for certain that the item of  interest has well-defined properties, but we don’t know what they are and hence the uncertainty. Such uncertainty can be quantified accurately using probability.

Vagueness, on the other hand, arises from an imprecise use of language.  Specifically, the term refers to the use of criteria that cannot distinguish between borderline cases.  Let’s clarify this using the example discussed earlier.  A popular way to measure customer satisfaction is through surveys. Such surveys may be able to tell us that customer A is more satisfied than customer B. However, they cannot distinguish between borderline cases because any boundary between satisfied and not satisfied customers is arbitrary.  This problem becomes apparent when considering an indifferent customer. How should such a customer be classified – satisfied or not satisfied? Further, what about customers who choose not to respond? It is therefore clear that any numerical probability computed from such data cannot be considered accurate.  In other words, the probability itself is uncertain.

### Ambiguity in classification

Although the distinction made by Colyvan is important, there is a deeper issue that can afflict uncertainties that appear to be quantifiable at first sight. To understand how this happens, we’ll first need to take a brief look at how probabilities are usually computed.

An operational definition of probability is that it is the ratio of the number of times the event of interest occurs to the total number of events observed. For example, if my manager notes my arrival times at work over 100 days and finds that I arrive before 8:00 am on 62 days then he could infer that the probability my arriving before 8:00 am is 0.62.   Since the probability is assumed to equal the frequency of occurrence of the event of interest, this is sometimes called the frequentist interpretation of probability.

The above seems straightforward enough, so you might be asking: where’s the problem?

The problem is that events can generally be classified in several different ways and the computed probability of an event occurring can depend on the way that it is classified. This is called the reference class problem.   In a paper entitled, “The Reference Class Problem is Your Problem Too”, Alan Hajek described the reference class problem as follows:

“The reference class problem arises when we want to assign a probability to a proposition (or sentence, or event) X, which may be classified in various ways, yet its probability can change depending on how it is classified.”

Consider the situation I mentioned earlier. My manager’s approach seems reasonable, but there is a problem with it: all days are not the same as far as my arrival times are concerned. For example, it is quite possible that my arrival time is affected by the weather: I may arrive later on rainy days than on sunny ones.  So, to get a better estimate my manager should also factor in the weather. He would then end up with two probabilities, one for fine weather and the other for foul. However, that is not all: there are a number of other criteria that could affect my arrival times – for example, my state of health (I may call in sick and not come in to work at all), whether I worked late the previous day etc.

What seemed like a straightforward problem is no longer so because of the uncertainty regarding which reference class is the right one to use.

Before closing this section, I should mention that the reference class problem has implications for many professional disciplines. I have discussed its relevance to project management in my post entitled, “The reference class problem and its implications for project management”.

### To conclude

In this post we have looked at a couple of forms of uncertainty about uncertainty that have practical implications for decision makers. In particular, we have seen that probabilities used in managerial decision making can be uncertain because of  vague definitions of events and/or ambiguities in their classification.  The bottom line for those who use probabilities to support decision-making is to ensure that the criteria used to determine events of interest refer to unambiguous facts that are appropriate to the situation at hand.  To sum up: decisions made on the basis of probabilities are only as good as the assumptions that go into them, and the assumptions themselves may be prone to uncertainties such as the ones described in this article.

Written by K

September 29, 2011 at 10:34 pm

## The drunkard’s dartboard: an intuitive explanation of Monte Carlo methods

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

Monte Carlo simulation techniques have been applied to areas ranging from physics to project management. In earlier posts, I discussed how these methods can be used to simulate project task durations (see this post and this one for example). In those articles, I described simulation procedures in enough detail for readers to be able to reproduce the calculations for themselves. However, as my friend Paul Culmsee mentioned, the mathematics tends to obscure the rationale behind the technique. Indeed, at first sight it seems somewhat paradoxical that one can get accurate answers via random numbers. In this post, I illustrate the basic idea behind Monte Carlo methods through an example that involves nothing more complicated than squares and circles. To begin with, however, I’ll start with something even simpler – a drunken darts player.

Consider a sozzled soul who is throwing darts at a board situated some distance from him.  To keep things simple, we’ll assume the following:

1. The board is modeled by the circle shown in Figure 1, and our souse scores a point if the dart falls within the circle.
2. The dart board is inscribed in a square with sides 1 unit long as shown in the figure, and we’ll assume for simplicity that the  dart always  falls somewhere  within the square (our protagonist is not that smashed).
3. Given his state, our hero’s aim is not quite what it should be –  his darts fall anywhere within the square with equal probability. (Note added on 01 March 2011: See the comment by George Gkotsis below for a critique of this assumption)

Figure 1: "Dartboard" and enclosing square

We can simulate the results of our protagonist’s unsteady attempts by generating two sets of  uniformly distributed random numbers lying between 0 and 1 (This is easily done in Excel using the rand() function).  The pairs of random numbers thus generated – one from each set –  can be treated as the (x,y) coordinates of  the dart for a particular throw. The result of 1000 pairs of random numbers thus generated (representing the drunkard’s dart throwing attempts) is shown in Figure 2 (For those interested in seeing the details, an Excel sheet showing the calculations for 100o trials can be downloaded here).

Figure 2: Result for 1000 throws

A trial results in a “hit” if it lies within the circle.  That is, if it satisfies the following equation:

$(x-0.5)^2 + (y-0.5)^2 < 0.25\ldots \ldots (1)$

(Note:  if we replace “<”  by “=”  in the above expression, we get the equation for a circle of radius 0.5 units, centered at x=0.5 and y=0.5.)

Now, according to the frequency interpretation of probability, the probability of the plastered player scoring a point is approximated by the ratio of the number of hits in the circle to the total number of attempts. In this case, I get an average of 790/1000 which is 0.79 (generated from 10 sets of 1000 trials each). Your result will be different from because you will generate different sets of random numbers from the ones I did. However, it should be reasonably close to my result.

Further, the frequency interpretation of probability tells us that the approximation becomes more accurate as the number of trials increases. To see why this is so, let’s increase the number of trials and plot the results. I carried out simulations for 2000, 4000, 8000 and 16000 trials. The results of these simulations are shown in Figures 3 through 6.

Figure 3: Result for 2000 throws

Figure 4: Result for 4000 throws

Figure 5: Result for 8000 throws

Figure 6: Result for 16000 throws

Since a dart is equally likely to end up anywhere within the square, the exact probability of a hit is simply the area of the dartboard (i.e. the circle)  divided by the entire area over which the dart can land. In this case, since the area of the enclosure (where the dart must fall) is 1 square unit, the area of the dartboard is actually equal to the probability.  This is easily seen by calculating the area of the circle using the standard formula $\pi r^2$ where $r$ is the radius of the circle (0.5 units in this case). This yields 0.785398 sq units, which is reasonably close to the number that we got for the 1000 trial case.  In the 16000 trial case, I  get a number that’s closer to the exact result: an average of 0.7860 from 10 sets of 16000 trials.

As we see from Figure 6, in the 16000 trial case, the entire square is peppered with closely-spaced “dart marks” – so much so, that it looks as though the square is a uniform blue.  Hence, it seems intuitively clear that as we increase, the number of throws, we should get a better approximation of the area and, hence, the probability.

There are a couple of points worth mentioning here. First, in principle this technique can be used to calculate areas of figures of any shape. However, the more irregular the figure, the worse the approximation – simply because it becomes less likely that the entire figure will be sampled correctly by “dart throws.” Second,  the reader may have noted that although the 16000 trial case gives a good enough result for the area of the circle, it isn’t particularly accurate considering the large number of trials. Indeed, it is known that the “dart approximation” is not a very good way of calculating areas – see this note for more on this point.

Finally, let’s look at connection between the general approach used in Monte Carlo techniques and the example discussed above  (I use the steps described in the Wikipedia article on Monte Carlo methods as representative of the general approach):

1. Define a domain of possible inputs – in our case the domain of inputs is defined by the enclosing square of side 1 unit.
2. Generate inputs randomly from the domain using a certain specified probability distribution – in our example the probability distribution is a pair of independent, uniformly distributed random numbers lying between 0 and 1.
3. Perform a computation using the inputs – this is the calculation that determines whether or not a particular trial is a hit or not (i.e.  if the x,y coordinates obey inequality  (1) it is a hit, else  it’s a miss)
4. Aggregate the results of the individual computations into the final result – This corresponds to the calculation of the probability (or equivalently, the area of the circle) by aggregating the number of hits for each set of trials.

To summarise: Monte Carlo algorithms generate random variables (such as probability) according to pre-specified distributions.  In most practical applications  one will use more efficient techniques to sample the distribution (rather than the naïve method I’ve used here.)  However, the basic idea is as simple as playing drunkard’s darts.

Acknowledgements

Thanks go out to Vlado Bokan for helpful conversations while this post was being written and to Paul Culmsee for getting me thinking about a simple way to explain Monte Carlo methods.

Written by K

February 25, 2011 at 4:51 am

## Monte Carlo simulation of risk and uncertainty in project tasks

### Introduction

When developing duration estimates for a project task, it is useful to make a distinction between the  uncertainty inherent in the task and uncertainty due to known risks.  The former is uncertainty due to factors that are not known whereas the latter corresponds uncertainty due to events that are known, but may or may not occur. In this post, I illustrate how the two types of uncertainty can be combined via Monte Carlo simulation.  Readers may find it helpful to keep my introduction to Monte Carlo simulations of project tasks handy, as I refer to it extensively in the present piece.

### Setting the stage

Let’s assume that there’s a task that needs doing, and the person who is going to do it reckons it will take between 2 and 8 hours to complete it, with a most likely completion time of 4 hours. How the estimator comes up with these numbers isn’t important here – maybe there’s some guesswork, maybe some padding or maybe it is really based on experience (as it should be).  For simplicity we’ll assume the probability distribution for the task duration is triangular. It is not hard to show that, given the above mentioned estimates, the probability, $p(t)$,  that the task will finish at time $t$ is given by the equations below (see my introductory post for a detailed derivation):

$p(t)=\frac{(t-2)}{6}\dots\ldots(1)$ for  2 hours $\leq t \leq$ 4 hours

And,

$p(t)=\frac{(8-t)}{12}\dots\ldots(2)$ for  4 hours $\leq t \leq$ 8 hours

These two expressions are sometimes referred to as the probability distribution function (PDF).  The PDF described by equations (1) and (2) is illustrated in Figure 1. (Note: Please click on the Figures to view full-size images)

Figure 1: Probability distribution for task

Now, a PDF tells us the probability that the task will finish at a particular time $t$. However, we are more interested in knowing whether or not the task will be completed by time $t$. – i.e. at or before time $t$. This quantity, which we’ll denote by $P(t)$ (capital P), is sometimes known as the cumulative distribution function (CDF). The CDF  is obtained by summing up the probabilities from $t=2$ hrs to time $t$.  It is not hard to show that the CDF for the task at hand is given by the following equations:

$P(t)=\frac{(t-2)^2}{12}\ldots\ldots(3)$ for  2 hours $\leq t \leq$ 4 hours

and,

$P(t)=1- \frac{(8-t)^2}{24}\ldots\ldots(4)$ for  4 hours $\leq t \leq$ 8 hours

For a detailed derivation, please see my introductory post. The CDF for the distribution is shown in Figure 2.

Now for the complicating factor: let us assume there is a risk that has a bearing on this task.  The risk could be any known factor that has a negative impact on task duration. For example, it could be that a required resource is delayed or that the deliverable will fails a quality check and needs rework. The consequence of the risk – should it eventuate – is that the task takes longer.  How much longer the task takes depends on specifics of the risk. For the purpose of this example we’ll assume that the additional time taken is also described by a triangular distribution with a minimum, most likely and maximum time of 1, 2 and 3 hrs respectively.  The PDF $p_{r}(t)$ for the additional time taken due to the risk is:

$p_{r}(t)=(t-1)\dots\ldots(5)$ for  1 hour $\leq t \leq$ 2 hours

And

$p_{r}(t)=(3-t)\dots\ldots(6)$ for  2 hrs $\leq t \leq$ 3 hours

The figure for this distribution is shown in Figure 3.

Figure 3: Probability distribution of additional time due to risk

The CDF for the additional time taken if the risk eventuates (which we’ll denote by $P_{r}(t)$) is given by:

$P_{r}(t)=\frac{(t-1)^2}{2}\ldots\ldots(7)$ for  1 hour $\leq t \leq$ 2 hours

and,

$P_{r}(t)=1- \frac{(3-t)^2}{2}\ldots\ldots(8)$ for  2 hours $\leq t \leq$ 3 hours

The CDF for the risk consequence is shown in Figure 4.

Figure 4: CDF of additional time due to risk

Before proceeding with the simulation it is worth clarifying what all this means, and what we want to do with it.

Firstly, equations 1-4 describe the inherent uncertainty associated with the task while equations 5 through 8 describe the consequences of the risk, if it eventuates.

Secondly, we have described the task and the risk separately. In reality, we need a unified description of the two – a combined distribution function for the uncertainty associated with the task and the risk taken together.  This is what the simulation will give us.

Finally, one thing I have not yet specified is the probabilty that the risk will actually occur. Clearly, the higher the probability, the greater the potential delay. Below I carry out simulations for risk probabilities of varying from 0.1 to 0.5.

That completes the specification of the  problem – let’s move on to the simulation.

### The simulation

The simulation procedure I used  for the zero-risk case  (i.e. the task described by equations 1 and 2 ) is as follows :

1. Generate a random number between 0 and 1.  Treat this number as the cumulative probability, $P(t)$ for the simulation run. [You can generate random numbers in Excel using the  rand() function]
2. Find the time, $t$,  corresponding to $P(t)$ by solving equations (3) or (4) for $t$. The resulting value of $t$ is the time taken to complete the task.
3. Repeat steps (1) and (2)  for a sufficiently large number of trials.

The frequency distribution of completion times for the task, based on 30,000 trials is shown in Figure 5.

Figure 5: Simulation histogram for zero-risk case

As we might expect,  Figure 5 can be translated to the probability distribution shown in Figure 1 by a straightforward normalization – i.e. by dividing each bar by the total number of trials.

What remains to be done is to  incorporate the risk (as modeled in equations 5-6) into the simulation. To simulate the task with the risk, we simply do the following for each trial:

1. Simulate the task without the risk as described earlier.
2. Generate another random number between 0 and 1.
3. If the random number is less than the probability of the risk, then simulate the  risk. Note that since the risk is described by a triangular function, the procedure to simulate it is the same as that for the task (albeit with different parameters).
4. If the random number is greater than the probability of the risk, do nothing.
5. Add the results of 1 and 4. This is the outcome of the trial.
6. Repeat steps 1-5 for as many trials as required.

I performed simulations for the task with risk probabilities of 10%, 30% and 50%. The frequency distributions of completion times for these are displayed in Figures 6-8 (in increasing order of probability). As one would expect, the spread of times increases with increasing probability. Further, the distribution takes on a distinct second peak as the probability increases: the first peak is at $t=4$, corresponding to the most likely completion time of the risk-free task and the second at $t=6$ corresponding to the most likely additional time of 2 hrs if the risk eventuates.

Figure 6: Simulation histogram (10% probability of risk)

Figure 7: Simulation histogram (30% probability of risk)

Figure 8: Frequency histogram (50% probability of risk)

It is also instructive to compare average completion times for the four cases (zero-risk and 10%, 30% and 50%). The average can computed from the simulation by simply adding up the simulated completion times (for all trials) and dividing by the total number of simulation trials (30,000 in our case). On doing this, I get the following:

Average completion time for zero-risk case = 4.66 hr

Average completion time with 10% probability of risk =  4.89 hrs

Average completion time with 30% probability of risk =  5.36 hrs

Average completion time with 50% probability of risk=  5.83 hrs

No surprises here.

One point to note is that the result obtained from the simulation for the zero-risk case compares well with the exact formula for a triangular distribution (see the Wikipedia article for the triangular distribution):

$t_{av} = \frac{t_{worst}+t_{best}+t_{most likely}}{3}=\frac{8+2+4}{3}=4.67$ hrs

This serves as a sanity check on the simulation procedure.

It is also interesting to compare the cumulative probabilities of completion in the zero-risk and high risk (50% probability) case. The CDFs for the two are shown in Figure 9. The co-plotted CDFs allow for a quick comparison of completion time predictions. For example, in the zero-risk case, there is about a  90% chance that the task will be completed in a little over 6 hrs whereas when the probability of the risk is 50%, the 90% completion time increases to 8 hrs (see Figure 9).

Figure 9: CDFs for zero risk and 50% probability of risk cases

### Next steps and wrap up

For those who want to learn more about simulating project uncertainty and risk, I can recommend the UK MOD paper – Three Point Estimates And Quantitative Risk Analysis A Process Guide For Risk Practitioners.  The paper provides useful advice on how three point estimates for project variables should be constructed. It also has a good discussion of risk and how it should be combined with the inherent uncertainty associated with a variable. Indeed, the example I have described above  was inspired by the paper’s discussion of uncertainty and risk.

Of course, as with any quantitative predictions of project variables, the numbers are only as reliable as the assumptions that go into them, the main assumption here being the three point estimates that were used to derive the distributions for the task uncertainty and risk (equations 1-2 and 5-6).  Typically these are obtained from historical data. However, there are well known problems associated with history-based estimates. For one, as we can never be sure that the historical tasks are similar to the one at hand in ways that matter (this is the reference class problem).  As Shim Marom warns us in this post, all our predictions depend on the credibility of our estimates.  Quoting from his post:

Can you get credible three point estimates? Do you have access to credible historical data to support that? Do you have access to Subject Matter Experts (SMEs) who can assist in getting these credible estimates?

If not, don’t bother using Monte Carlo.

In closing, I hope my readers will find this simple example useful in understanding how uncertainty and risk can be accounted for using Monte Carlo simulations. In the end, though, one should always keep in mind that the use of sophisticated techniques does not magically render one immune to the GIGO principle.

Written by K

February 17, 2011 at 9:55 pm