Category Archives: R

Monty Hall problem – a small simulation in R

The Monty Hall problem is an interesting example for how much intuition can mislead us in some statistical contexts. Even more disturbing though is, for how long we are prepared to debate and defend an expected result before actually checking our initial guesses using a simple Monte Carlo simulation.

Here is simple simulation implementation the Monty Hall game show problem:

In the TV show “Let’s Make a Deal” the host Monty Hall would offer a game participant the choice of three doors. One of them was hiding a valuable price (eg a car) – behind the other two doors were only two less desirable goats.

After the participant has choose a door, which is still kept closed, Monty would open one of the other doors showing that it did not hide the price. After that door was open the participant was asked if he/she wanted to maintain their initial choice or rather change their mind and now swap to the other of the two still closed doors.

It turns out, that swapping the initial choice of door is in fact a better strategy in this game, which is considered by many counter-intuitive. The public discussion about this problem in a print magazine revealed a lot about the excitement and sometime rather unscientific culture when debating a “scientific” question via public media. And twitter was not even available.


Here a small Monte Carlo toy in R to compare between the two game strategies. It uses mainly two features:

  • the sample() function to pick randomly – here from a vector of doors.
  • vectors indexed with negative indices: evaluate to all but the “subtracted” elements
doors     <- 1:3  # 3 doors named 1 to 3
N         <- 1000 # we play N times

swap_wins <- 0    # set win counter to 0

for (i in 1:N) {
  price  <- sample(doors,1)                   # randomly pick a door for the price
  player <- sample(doors,1)                   # player randomly picks a door
  monty  <- sample(doors[-c(player,price)],1) # monty picks from remaining doors
                                              # eg price and player doors removed

  swap_wins <- (price != player) + swap_wins  # did swap strategy win this time?
}                                             # (price was not behind initial choice)

message("fraction of swap wins: f",round(swap_wins/N,3))

The result of around 0.66 confirms that indeed swapping the door after Monty has opened another one is a better strategy than sticking to the initial choice, which amounts only to 1/3 probability of success.

Such a simulation is of course neither the only nor the most elegant way to arrive at this result. But in case of doubt it looks definitely easier to run a quick test than to embark in a public discussions. At the time this involved the exchange of 10,000 reader letters in the US and proportinal excitement in similar thread in Germany. For more details check out the wikipedia article here or a small book about this topic (in German) here.

PS: Did you notice that the simulation program does not actually use the line that determines Monty’s pick? A hint to an analytical solution.

Slides and blog posts with R and emacs org-mode

Preparing a larger number of slides with R code and plots can be a bit tedious with standard desktop presentation software like powerpoint or keynote. The manual effort to change the example code, run the analysis and then cut and paste updated graphs, tables and code is high. Sooner or later one is bound to create inconsistencies between code and expected results or even syntax errors

Using emacs and its swiss army knife org-mode there is another elegant and reproducible solution: just export the consistent code and output from an org-mode code block to generate either slides or bog pages (or both) from a single source. As a free benefit freen this tool chain can also produce high quality PDF handouts with a managed table of content and an index.

Org Mode source blocks

Here is a standard statistics “hello world” example using ggplot2:


df <- data.frame( norm = rnorm(10000))
ggplot(df, aes(x = norm)) + geom_histogram()


Continue reading

Using Data Frames in Feather format (Apache Arrow)

Triggered by the RStudio blog article about feather I did the one line install and compared the results on a data frame of 19 million rows. First results look indeed promising:

# build the package
> devtools::install_github("wesm/feather/R")

# load an existing data frame (19 million rows with batch job execution results)
> load("batch-12-2015.rda")

# write it in feather format...
> write_feather(dt,"batch-12-2015.feather")

# ... which is not compressed, hence larger on disk
> system("ls -lh batch-12-2015.*")
-rw-r--r-- 1 dirkd staff 813M 7 Apr 11:35 batch-12-2015.feather
-rw-r--r-- 1 dirkd staff 248M 27 Jan 22:42 batch-12-2015.rda

# a few repeat reads on an older macbook with sdd
> system.time(load("batch-12-2015.rda"))
user system elapsed
8.984 0.332 9.331
> system.time(dt1 <- read_feather("batch-12-2015.feather"))
user system elapsed
1.103 1.094 7.978
> system.time(load("batch-12-2015.rda"))
user system elapsed
9.045 0.352 9.418
> system.time(dt1 <- read_feather("batch-12-2015.feather"))
user system elapsed
1.110 0.658 3.997
> system.time(load("batch-12-2015.rda"))
user system elapsed
9.009 0.356 9.393
> system.time(dt1 <- read_feather("batch-12-2015.feather"))
user system elapsed
1.099 0.711 4.548

So, around half the elapsed time and about 1/10th of the user cpu time (uncompressed) ! Of course these measurements are from file system cache rather than the laptop SSD, but the reduction in wall time is nice for larger volume loads.

More important though is the cross-language support for R, Python, Scala/Spark and others, which could make feather the obvious exchange format within a team or between workflow steps that use different implementation languages.

Setting up an RStudio server for iPad access

Sometimes it can be convenient to run RStudio remotely from an iPad or another machine with little RAM or disk space. This can be done quite easily using the free RStudio Server on OSX via docker. To do this:

  • Find the rocker/rstudio image on docker hub and follow the setup steps here [github].
  • Once the image is running, you should be able to connect with Safari on the host Mac to the login page eg at
    $ open
  • Now there is is only a small last step needed. You need to expose the server port from the host on the local network using the OSX firewall. In the somewhat explicit language of the “new” OSX firewall this can be done using:

    $ echo "rdr pass inet proto tcp from any to any port 8787 -> port 8787" | sudo pfctl -ef -

    At this point you should be able to connect remotely from your iPad to


    and continue your R session where you left it before eg on your main machine.

    BTW: If your network can not be trusted then you should probably change the default login credentials as described in the image docs.

    Cached, asychronous IP resolving

    Resolving IP addresses to host names is quite helpful for getting a quick overview of who is connecting from where. This may need some care to not put too much strain on your DNS server with a large number of repeated lookups. Also you may not want to wait for timeouts on IPs that do not resove. R itself is not supporting this specifically but can easily exploit asyncronous DNS lookup tools like adns (on OSX from homebrew) and provide a cache to speed things up. Here is an simple example for a vertorised lookup using a data.table as persistent cache.

    ## this basic aysnc lookup is a modified version of an idea described in
    ## <- function(ips) {
      ## store ip list in a temp file
      tf <- tempfile()
      cat(ips, sep='\n', file=tf)
      ## use the adns filter to resolve them asynchronously (see man page for timeouts and other options)
      host.names <- system(paste("adnsresfilter <", tf) ,intern=TRUE, ignore.stderr=TRUE)
      ## cleanup the temp file
    ## now extend the above to implement a  ip to name cache
    ip.cached.lookup <- function(ips, reset.cache=FALSE) {
      cache.file <- "~/.ip.cache.rda"
      ## if the cache file exists: load it
      if (!reset.cache & !file.access(cache.file,4)){
          message("ip cache entries loaded :", nrow(host))
      } else {
          ## create an empty table (with just locahost)
          host <- data.table(hip="", hname="localhost")
      ## prepare a table of query ip and name
      qh <- data.table(hip=as.character(ips),hname=NA)
      ## keep them sorted by ip to speedup data.table lookups
      ## resolve all known host name from the cache
      qh$hname <- host[qh]$hname
      ## collect the list of unique ips which did not get resolved yet
      new.ips <- unique(qh[$hname)]$hip)
      ## if not empty, resolve the rest
      if (length(new.ips) > 0) {
        ## add the new ips to the cache table
        host <- rbind(host, list(hip=new.ips,hname=NA))
        ## find locations which need resolution (either new or expired)
        need.resolving <-$hname)
        message("new ips to resolve: ", sum(need.resolving))
        ## and resolve them
        host$hname[need.resolving] <-[need.resolving]$hip)
        ## need to set key again after rbind above..
        ## .. to do the remaining lookups
        qh$hname <- host[qh]$hname
        ## save the new cache status
        save(host, file = cache.file)
    ## with this function you can easily add a column to your
    ## weblog data.table from the previous posts to get started with
    ## the real log analysis
    w$ <- ip.cached.lookup(w$host)

    Getting hold of remote weblogs

    The last post was assuming that the weblogs to analyse are directly accessible by the R session which may not be the case if your analysis is running on a remote machine. Also in some cases you may want to filter out some uninteresting log records (eg local clients on the web server or local area accesses from known clients). The next examples show how to modify the previous R script using the R pipe function to take this into account:

    ## read the last 100K log entries from svr via a ssh connection
    ## (this assumes you have setup the ssh keys correctly beforehand)
    w <- data.table(read.table(pipe("ssh svr 'tail -n 100000 /var/log/apache2/access_log'")))
    ## in addition filter out all accesses from local clients on the web
    ## server or the local subnet (in this case
    w <- data.table(read.table(pipe("ssh svr 'tail -n 100000 /var/log/apache2/access_log | awk \"\\$2 !~ /127\\.0\\.0\\.1|192\\.168\\.10\\./\"'")))
    ## note: the proper quoting/escaping of R and shell strings on this one takes
    ## more effort than the processing. There must be an R function which does this...

    In a similar way you could concatenate multiple (eg already logrotated) logs and/or unzip logfiles. As this pre-filtering takes place locally on the server machine holding the log files this helps to bring down the data amount to be transfered and analysed: always a good start to avoid the popular ‘unecessarily big data’ syndrome…

    Using R for weblog analysis

    Apache Weblog Analysis

    Whether you run your own blog or web server or use some hosted service – at some point you may be interested in some information on how well your server or your users are doing. Many infos like hit frequency, geolocation of users and distribution of spent bandwidth are very useful for this and can be obtained in different ways:

    • by instrumenting the page running inside the client browser (eg piwik)
    • by analysis of the web server logs (eg webalizer)

    For the latter I have been using for several years webalizer, which does nice web based analysis plots. More recently I moved to a more complicated server environment with several virtual web services and I found the configuration and data selection options a bit limting. Hence I started as a toy project to implement the same functionality with a set of simple R scripts, which I will progressively share here.

    As a first step some simple examples for the data import, cleaning and overview plots. We’ll then add anychronous IP resolution, add and analyse goelocation information and as a last step wrap the analysis output tables and plots into a web application, which can be consulted from a remote browser.

    data.table vs. dplyr

    One of my favourite R packages for data handling, which I will use also here is the ‘data.table’ package. Note: Most of the results can be obtained in a similar way also using the excellent ‘dplyr’ package, but for some of my other (larger volume) projects data.table has some performance and memory efficiency advantages, so I’ll stick to data.table. If you are using R for data handling/aggregating and are not familar with either packages – take a look at both and make your own choice.

    Importing the logs into R

    Well, this part is rather simple since apache logs can be read via the standard read.table function:

    ## read the complete log - your file name is likely different
    w <- data.table(read.table("/var/log/apache2/access_log"))
    ## there are a few different log types which vary in the number and sequence
    ## if log items. Have a look at the apache configuration or just the file.
    ## In my case I get a so called 'combinedvhost' file which lists in the first
    ## two columns the website (out of several virtual sites on the some server)
    ## and as second field the client host which accessed the server.
    ## There is a good chance that your server config does omit the first field
    ## so you may try to drop the 'vhost' string below.
    ## try the following command to see if data and field names match:
    ## btw: already this summary shows a lot of interesting info