Thursday, September 13, 2012

Estimating Pi with R via MCS-Dart: A very simple example of numerical integration, illustrated and computed in R.

Have you ever played Monte Carlo Dart? If not, read this post and learn how to do it with R and what it can be used for. In fact it is a very easy and prevalent example (which I have come across in a  computational economics course last spring semester) that demonstrates the idea behind numerical integration. The aim of this post is not to present a programmatically swell implementation but rather to give an illustrative explanation. Having that said, here's the problem:
You know of the existence of the number Pi and its role in determining the area of a circle. But, you don't know its value. How can you estimate the value of Pi using Monte Carlo Simulation, or in other words using random numbers?

Square and circle

First, draw a circle with radius r=1 placed in a square with side length a=2r=2 to illustrate the initial setting. Note that I use the shape package to draw the circle.




library(shape)

r <- 1
par(pty="s")
emptyplot(c(0, 0),frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol="blue")
lines(x=c(0,1),y=c(0,0), lty=1, col="red")
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)
text(x=0.5,y=0.2, labels=expression("r=1"), col="red")


Generally, if a circle with radius r is surrounded by a square with side length a=2r, Pi is equal to the area of the circle divided by the area of the square times 4. Hence, by estimating the fraction of the two areas multiplied by 4 we estimate the value of Pi.


# in case you want to plot the math in R...

emptyplot( xlim=c(0,10), ylim=c(0,11))
text(5,10, expression(paste(A[circle] == r^2 %*% pi, "  and  ", A[square] == (2 * r)^{2},"; ", sep="") ))

text(5,7, expression( paste( frac(A[circle], A[square]) == frac(r^2 %*% pi, (2 * r)^{2}), " = ", frac(pi, 4), ";  thus  ", 4 %*% frac(A[circle], A[square]) == pi, sep="") ))


Playing Monte Carlo Dart

To estimate the fraction of the two areas, we use Monte Carlo Simulation. Imagine it (in this context) as a simulation of randomly throwing darts at the graph above and check for each dart whether it is inside or outside of the circle. The fraction of darts inside the circle is approximately the area of the circle relative to the area of the square (and thus, multiplied by 4, approximately the value of Pi). Intuitively this approximation gets better the greater the number of darts (random numbers) is. But even with a high number of darts the estimated value might in one execution of the experiment be above and in another under the true value of Pi. On average the estimated value is, though, unbiased and pretty close to the real value. Hence, to get a fair estimation we repeat the experiment many times and take the mean of the results.
In fact, the approach described above can be even more simplified. Instead of looking at the whole circle/square it is sufficient to only look in the upper right fourth of the circle/square (as illustrated below). With that target we do exactly the same "dart-throwing" and multiply the resulting fraction by 4.





r <- 1
par(pty="s")
emptyplot(frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol="blue")
lines(x=c(0,1),y=c(0,0), lty=1, col="red")
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)

How to "throw the darts"?

We draw two vectors of random numbers that are ~U(0,1). One for the x coordinates and one for the y coordinates. Each combination of the elements of these two vectors gives us the location where a "dart hit the target". The remaining question is how to determine whether the dart hit the square inside or outside the circle. The answer comes from the Pythagorean theorem. Imagine the distances from the axes to the point (0.8,0.58) as side lengths a and b of a triangle. If the value of the unknown side length c is greater than the radius r (1), the point is outside the circle; otherwise it is inside (or on) the circle (see illustration below).



 

# Illustration Pythagorean
par(pty="s")
emptyplot(frame.plot=r)
plotcircle(r=r, mid=c(0,0), lwd=1,lcol="blue")
lines(x=c(0,1),y=c(0,0), lty=1, col="red")
lines(x=c(-1,1),y=c(0,0), lty=2)
lines(y=c(-1,1),x=c(0,0), lty=2)
axis(1)
axis(2)

points(x=0.8,y=0.58,col="red", lwd=2)
lines(x=c(0,0.8), y=c(0,0), col="red", lwd=2)
lines(x=c(0.8,0.8), y=c(0,0.58), col="red", lwd=2)
lines(x=c(0,0.8), y=c(0,0.58), col="red", lty=2,lwd=2)
text(x=0.4,y=0.4,"c",col="red")
title(expression(c == sqrt(a^2 + b^2)))

A function for one round of Monte Carlo Dart

The following function estimates Pi with the approach described above. The only input to the function is the number of "darts" (random coordinates) n. Note that this implementation is programmatically not smart because it uses a loop instead of vectorization. But, I think it is more illustrative for (R-) programming beginners this way.
  

est.pi <- function(n){
inside <- 0
outside <-0
pi.est <- 0

# draw random numbers for the coordinates of the "dart-hits"
a <- runif(n,0,1)
b <- runif(n,0,1)
# use the pythagorean theorem
c <- sqrt((a^2) + (b^2) )

# loop: check for each point whether it is inside or outside the circle:

for (i in 1:n) {
   
   if (c[i] > 1) {
     
      outside <- outside + 1
     
    } else {
     
      inside <- inside + 1
     
    }
       
    if (i==n) {pi.est <- (inside/(inside+outside))*4} # compute an estimation of the value of pi
    
  }
 return(pi.est)
}


Try that function several times. For the same n you will every time get different estimations of the value of Pi. That is due to the fact, that the whole approach is based on random numbers.
   
est.pi(n=1000)
est.pi(n=1000)
est.pi(n=1000)


Now, repeat this experiment 1000 times and take the mean of the results as an unbiased estimation of the value of pi. Compare the result with the value of Pi defined in R.
   

rep <- 1000

pis <- sapply(1:rep, function(x){
  est.pi(n=1000)
})

mean(pis)
pi # value of pi defined in R


Finally, get a graphical overview of the estimation.
  


plot(density(pis), main=expression(paste("Distribution of estimated ",pi,"-values")),cex.main=0.7, cex.axis=0.7, cex.lab=0.7)

# the value of pi defined in R
lines(x=c(pi,pi),y=c(0,10), lty=1,col="blue")
# the estimation
lines(x=c(mean(pis),mean(pis)),y=c(0,10), lty=1,col="red")


What has this to do with numerical integration?

The link between the application of "MCS-Dart" to estimate the value of Pi and the general idea behind numerical integration is straightforward. The core of the approach described above is to draw random numbers (of a certain interval) for coordinates and check, whether the resulting points are under or above a curve. If one is interested in the (approximate) area under a curve of a given function y=f(x) =... in a certain interval and it is very hard to find the integral analytically one can apply exactly the same approach. Draw random numbers for x and y, check if y is greater than f(x), if so, the point (x,y) must be above the curve and otherwise under the curve. Compute the fraction of points under the curve and multiply it by the total area of the rectangle with the x- and y-intervals as side lengths...

Thursday, August 23, 2012

R and the web (for beginners), Part III: Scraping MPs' expenses in detail from the web

In this last post of my little series (see my latest post) on R and the web I explain how to extract data of a website (web scraping/screen scraping) with R. If the data you want to analyze are a part of a web page, for example a HTML-table (or hundreds of them) it might be very time-consuming (and boring!) to manually copy/paste all of its content or even typewrite it to a spreadsheet table or data frame. Instead, you can let R do the job for you!

This post is really aimed at beginners. Thus, to keep things simple it only deals with scraping one data table from one web page: a table published by BBC NEWS containing the full range of British Members of Parliament' expenses in 2007-2008. Quite an interesting data set if you are into political scandals...


Web scraping with R

There are several R packages that might be helpful for web scraping, such as XML, RCurl, and scrapeR. In this example only the XML package is used. As a fist step, you parse the whole HTML-file and extract all HTML-tables in it:

library(XML)

# URL of interest:
mps <- "http://news.bbc.co.uk/2/hi/uk_politics/8044207.stm" 

# parse the document for R representation:
mps.doc <- htmlParse(mps)

# get all the tables in mps.doc as data frames
mps.tabs <- readHTMLTable(mps.doc) 

mps.tabs is a list containing in each element a HTML-table from the parsed website (mps.doc) as data.frame. The website contains several HTML-tables (some are rather used to structure the website and not to present data). The list mps.tabs actually has seven entries, hence there were seven HTML-tables in the parsed document:

length(mps.tabs)

To proceed you need to check which of these data frames (list entries) contains the table you want (the MPs' expenses). You can do that "manually" by checking how the data frame starts and ends and compare it with the original table of the website:

head(mps.tabs[[1]])  #and
tail(mps.tabs[[1]])  #for 1 to 7

With only seven entries this is quite fast. But alternatively you could also write a little loop to do the job for you. The loop checks each data frame for certain conditions. In this case: the string of the first row and first column and the string in the last row and column. According to the original table from the website that should be:
first <- "Abbott, Ms Diane"
last <- "157,841"

# ... and the loop:

for (i in 1:length(mps.tabs)) {
 
  lastrow <- nrow(mps.tabs[[i]]) # get number of rows
  lastcol <- ncol(mps.tabs[[i]])
 
  if (as.character(mps.tabs[[i]][1,1])==first & as.character(mps.tabs[[i]][lastrow,lastcol])==last) {
   
    tabi <- i
     
    }
  }

Check if that is realy what you want and extract the relevant table as data frame.

head(mps.tabs[[tabi]])
tail(mps.tabs[[tabi]])
mps <- mps.tabs[[tabi]] 

Before you can properly analyze this data set we have to remove the commas in the columns with expenses and format them as numeric:

money <- sapply(mps[,-1:-3], FUN= function(x) as.numeric(gsub(",", "", as.character(x), fixed = TRUE) ))

mps2 <- cbind(mps[,1:3],money)

Now you are ready to go... For example, you could compare how the total expenses are distributed for each of the five biggest parties:

# which are the five biggest parties by # of mps?
nbig5 <- names(summary(mps2$Party)[order(summary(mps2$Party)*-1)][1:5])

#subset of mps only with the five biggest parties:
big5 <- subset(mps2, mps$Party%in%nbig5)

# load the lattice package for a nice plot

library(lattice)

bwplot(Total ~  Party, data=big5, ylab="Total expenses per MP (in £)")


Here is the resulting plot: 



 And the relevant R code in one piece:

library(XML)

# URL of interest:
mps <- "http://news.bbc.co.uk/2/hi/uk_politics/8044207.stm" 

# parse the document for R representation:
mps.doc <- htmlParse(mps)


# get all the tables in mps.doc as data frames
mps.tabs <- readHTMLTable(mps.doc)
# loop to find relevant table:

first <- "Abbott, Ms Diane"
last <- "157,841"

for (i in 1:length(mps.tabs)) {
 
  lastrow <- nrow(mps.tabs[[i]]) # get number of rows
  lastcol <- ncol(mps.tabs[[i]])
 
  if (as.character(mps.tabs[[i]][1,1])==first & as.character(mps.tabs[[i]][lastrow,lastcol])==last) {
   
    tabi <- i
     
    }
  }


# extract the relevant table and format it:

mps <- mps.tabs[[tabi]]  

money <- sapply(mps[,-1:-3], FUN= function(x) as.numeric(gsub(",", "", as.character(x), fixed = TRUE) ))

mps2 <- cbind(mps[,1:3],money)


# which are the five biggest parties by # of mps?
nbig5 <- names(summary(mps2$Party)[order(summary(mps2$Party)*-1)][1:5])

#subset of mps only with the five biggest parties:
big5 <- subset(mps2, mps$Party%in%nbig5)

# load the lattice package for a nice plot

library(lattice)

bwplot(Total ~  Party, data=big5, ylab="Total expenses per MP (in £)")

Friday, June 22, 2012

R and the web (for beginners), Part II: XML in R


This second post of my little series on R and the web deals with how to access and process XML-data with R. XML is a markup language that is commonly used to interchange data over the Internet. If you want to access some online data over a webpage's API you are likely to get it in XML format. So here is a very simple example of how to deal with XML in R.
Duncan Temple Lang wrote a very helpful R-package which makes it quite easy to parse, process and generate XML-data with R. I use that package in this example. The XML document (taken from w3schools.com) used in this example describes a fictive plant catalog. Not that thrilling, I know, but the goal of this post is not to analyze the given data but to show how to parse it and transform it to a data frame. The analysis is up to you...

How to parse/read this XML-document into R?
 
# install and load the necessary package

install.packages("XML")
library(XML)


# Save the URL of the xml file in a variable

xml.url <- "http://www.w3schools.com/xml/plant_catalog.xml"

# Use the xmlTreePares-function to parse xml file directly from the web
 
xmlfile <- xmlTreeParse(xml.url)


# the xml file is now saved as an object you can easily work with in R:

class(xmlfile)



# Use the xmlRoot-function to access the top node

xmltop = xmlRoot(xmlfile)

# have a look at the XML-code of the first subnodes:

print(xmltop)[1:2]

This should look more or less like:


$PLANT
<PLANT>
 <COMMON>Bloodroot</COMMON>
 <BOTANICAL>Sanguinaria canadensis</BOTANICAL>
 <ZONE>4</ZONE>
 <LIGHT>Mostly Shady</LIGHT>
 <PRICE>$2.44</PRICE>
 <AVAILABILITY>031599</AVAILABILITY>
</PLANT>

$PLANT
<PLANT>
 <COMMON>Columbine</COMMON>
 <BOTANICAL>Aquilegia canadensis</BOTANICAL>
 <ZONE>3</ZONE>
 <LIGHT>Mostly Shady</LIGHT>
 <PRICE>$9.37</PRICE>
 <AVAILABILITY>030699</AVAILABILITY>
</PLANT>

attr(,"class")
[1] "XMLNodeList"

One can already assume how this data should look like in a matrix or data frame. The goal is to extract the XML-values from each XML-tag <> for all $PLANT nodes and save them in a data frame with a row for each plant ($PLANT-node) and a column for each tag (variable) describing it. How can you do that?


# To extract the XML-values from the document, use xmlSApply:

plantcat <- xmlSApply(xmltop, function(x) xmlSApply(x, xmlValue))


# Finally, get the data in a data-frame and have a look at the first rows and columns

plantcat_df <- data.frame(t(plantcat),row.names=NULL)
plantcat_df[1:5,1:4]

The first rows and columns of that data frame should look like this:
 
               COMMON              BOTANICAL ZONE        LIGHT
1           Bloodroot Sanguinaria canadensis    4 Mostly Shady
2           Columbine   Aquilegia canadensis    3 Mostly Shady
3      Marsh Marigold       Caltha palustris    4 Mostly Sunny
4             Cowslip       Caltha palustris    4 Mostly Shady
5 Dutchman's-Breeches    Dicentra cucullaria    3 Mostly Shady
Which is exactly what we need to analyze this data in R.



Thursday, June 21, 2012

R and the web (for beginners), Part I: How is the local nuclear plant doing?



One of the things I especially like about R is its ability to easily access and process data from the web. If you are new to R, or never have used it to access data from the Internet, here is the first part of a little series of posts with examples to get you started. This first post gives a very simple example of how to access a data set that is saved online.

This might be particularly useful if the data set at hand is frequently updated and you want to repeatedly generate a statistical report of some kind based on this data set. Hence, having the analysis and the link to the data in one R-script means you only have to rerun the script whenever you want to update your report or anybody else wants to reproduce it. The data I'm using for this example is exactly of that type. It's a file published by the United States Nuclear Regulatory Commission (U.S. NRC) reporting the power reactor status of U.S. nuclear power plants for the last 365 days, thus it is updated every day. 

How to access that data directly through R? 

# First: save the url of the data file as character string

url.npower <- "http://www.nrc.gov/reading-rm/doc-collections/event-status/reactor-status/PowerReactorStatusForLast365Days.txt"


# then read the data file into R (in this case the data file is a text file with "|" separating the columns) 

npower <- read.table(url.npower, sep="|", header=TRUE)

# and format the date column

npower$ReportDt <- as.Date(npower$ReportDt, format="%m/%d/%Y")


The data set is now ready for analysis. For example: a graphical analysis of the recent power reactor status of some of the nuclear power plants:

# load the necessary lattice package
# (if it isn't installed yet, run: install.packages("lattice")
library(lattice)

# take a subset of the data
sample <- npower[npower$Unit==as.character(unique(npower$Unit)[1:24]),]

# get a graphical overview
xyplot(Power~ReportDt | Unit, data=sample, type="l",col.line="black", xlab="Time",ylab="Power" )



Save the code above in an R-script, rerun it some days later, and your graphical analysis will be up to date.