mages' blog

Rasmus' post of last week on binomial testing made me think about p-values and testing again. In my head I was tossing coins, thinking about gender diversity and toast. The toast and tossing a buttered toast in particular was the most helpful thought experiment, as I didn't have a fixed opinion on the probabilities for a toast to land on either side. I have yet to carry out some real experiments.

Suppose I tossed 6 buttered toasts and observed that all but one toast landed on their buttered side.

Now I have two questions:

• Would it be reasonable to assume that there is a 50/50 chance for a toast to land on either side?
• Which probability should I assume?

If I believe the toast is fair, then the probability for landing on the buttered (B) or unbuttered (U) side is 50%.

The probability of observing one ore more B (right tail event) is:

(gt <- 1-1/2^6)
# [1] 0.984375

and the probability of observing one or fewer B (left tail event) is:

(lt <- 1/2^6*choose(6,1) + 1/2^6)
# [1] 0.109375

while the probability of either extreme event, one or fewer B (or U), is:

2*min(c(gt, lt))
# [1] 0.21875

In summary, if the toast has an equally probability to fall on either side, then there is 22% chance to observe one or fewer B (or U) in 6 tosses. That's not that unlikely and hence I would not dismiss the hypothesis that the toast is fair, despite the fact that the sample frequency is only 1/6.

Actually, the probabilities I calculated above are exactly the p-values I get from the classical binomal tests:

## Right tail event
binom.test(1, 6, alternative="greater")
## Left tail event
binom.test(1, 6, alternative="less")
## Double tail event
binom.test(1, 6, alternative="two.sided")

Additionally I can read from the tests that my assumption of a 50% probability is on the higher end of the 95 percent confidence interval. Thus, wouldn't it make sense to update my belief about the toast following my observations? In particular, I am not convinced that a 50/50 probability is a good assumption to start with. Arguably the toast is biased by the butter.

Here the concept of a conjugate prior becomes handy again. The idea is to assume that the parameter \(\theta\) of the binomial distribution is a random variable itself. Suppose I have no prior knowledge about the true probability of the toast falling on either side, then a flat prior, such as the Uniform distribution would be reasonable. However, the beta distribution with parameter \(\alpha=1\) and \(\beta=1\) has the same property and is a conjugate to the binomial distribution with parameter \(\theta\). That means there is an analytical solution, in this case the posterior distribution is beta-binomial with hyperparaemters:
\[\alpha':=\alpha + \sum_{i=1}^n x_i,\; \beta':=\beta + n - \sum_{i=1}^n x_i,\]
and the posterior predictor for one trial is given as
\[\frac{\alpha'}{\alpha' + \beta'}\]
so in my case:

alpha <- 1; beta <- 1; n <- 6; success <- 1
alpha1 <- alpha + success
beta1 <- beta + n - success
(theta <- alpha1 / ( alpha1 + beta1))
# [1] 0.25

My updated believe about the toast landing on the unbuttered side is a probability of 25%. That's lower than my prior of 50% but still higher than the sample frequency of 1/6. If I would have more toasts I could run more experiments and update my posterior predictor.

I get the same answer from Rasmus' bayes.binom.test function:

> library(BayesianFirstAid)
> bayes.binom.test(1, 6)

 Bayesian first aid binomial test

data: 1 and 6
number of successes = 1, number of trials = 6
Estimated relative frequency of success:
  0.25 
95% credible interval:
  0.014 0.527 
The relative frequency of success is more than 0.5 by a probability of 0.061 
and less than 0.5 by a probability of 0.939 

Of course I could change my view on the prior and come to a different conclusion. I could follow the Wikipedia article on buttered toast and believe that the chance of the toast landing on the buttered side is 62%. I further have to express my uncertainty, say a standard deviation of 10%, that is a variance of 1%. With that information I can update my belief of the toast landing on the unbuttered side following my observations (and transforming the variables):

x <- 0.38
v <- 0.01
alpha <- x*(x*(1-x)/v-1)
beta <- (1-x)*(x*(1-x)/v-1)
alpha1 <- alpha + success
beta1 <- beta + n - success
(theta <- alpha1 / ( alpha1 + beta1))
# [1] 0.3351821

I would conclude that for my toasts / tossing technique the portability is 34% to land on the unbuttered side.

In summary, although their is no sufficient evidence to reject the hypothesis that the 'true' probability is not 50% (at the typical 5% level), I would work with 34% until I have more data. Toast and butter please!

Session Info

R version 3.0.1 (2013-05-16)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] BayesianFirstAid_0.1 rjags_3-12           coda_0.16-1         
[4] lattice_0.20-24     

loaded via a namespace (and not attached):
[1] grid_3.0.1  tools_3.0.1

Since Christmas I have been playing around with a Raspberry Pi. It is certainly not the fastest computer, but what a great little toy! Here are a few experiences and online resources that I found helpful.

Setup

Initially I connected the Raspberry Pi via HDMI to a TV; together with keyboard, mouse and an old USB Wifi adapter. Everything worked out of the box and I could install Raspbian and set up the network.

HDMI to VGA

Using an old VGA computer monitor via an adapter required changes to the file /boot/config.txt. You can find the parameters that match your monitor on Raspberry Pi StackExchange. In my case I had to set:

hdmi_group=2
hdmi_mode=35 # 1280x1024 @ 60Hz

Remote Access

But who needs a monitor when you can access the Pi remotely anyway? The command ifconfig tells me the local IP address of the Raspberry Pi.

With XQuartz on my Mac running I can connect to the Raspberry Pi via ssh with the X session forwarded:

ssh -X pi@your.ip.address.here

However, the performance is a bit sluggish and online comments suggest to use VNC instead. Nothing easier than that, install the VNC server on the Pi and use Screen Sharing on your Mac to access the Pi. Mitch Malone has a great post on this subject. Following the VNC setup on the Raspberry Pi I can type:

vnc://pi@your.ip.address.here:5901

into Safari and that will bring up the Screen Sharing App; see screen shot below.

AirPrint and AirPlay

Ok, let's give the Pi something to do: Rohan Kapoor explains how to set up the Raspberry Pi as a print server with AirPrint.

How about AirPlay as well? Follow Thorin Klosowski's steps on Lifehacker and you can stream music from your iOS devices to the Raspberry Pi's audio out.

Mathematica and R

Just before Christmas Stephen Wolfram announced that Mathematica would be made freely available for the Raspberry Pi. I had used Mathematica a little at university and was curious to see it on the Pi. Alex Newman posted the installation instructions on the Wolfram community site. It is as simple as:

sudo apt-get update && sudo apt-get install wolfram-engine

As a little toy example I run Paul Nylander's Mathematica code for a Feigenbaum diagram:

Surprisingly, it took about 3.5 minutes to run. Curious to find out if my old R routine would run faster I installed R on my Pi:

sudo apt-get install r-base

and adapted the code to match the parameters of the Mathematica routine (well, so I think):

The R code finished after about 10 seconds. Surely, there must be ways to speed up the Mathematica code that I have to investigate.

I noticed that the monthly number of posts on R-bloggers stopped increasing over the last year. Indeed, the last couple of months saw a decline in posts compared to the previous year. Thus, has most been said and written about R already?


Who knows? Well, I took a stab at looking into the future. However, I can tell you already that I am not convinced by my predictions. But maybe someone else will be inspired to take this work forward.

First, I have to get the data - that's easy, I can scrape the monthly post counts from the R-bloggers homepage.


Looking at the incremental and cumulative plots, and believing that eventually the number of R posts will decrease, I thought that a logistic growth function would provide a nice fit to the data and also give an asymptotic view of the total number of posts on R-bloggers.

Although the fit, see below, looks reasonable at first glance, I don't believe it provides a sensible prediction of the future. The model would forecast only another 1,269 post by the end of 2016 with not much more to expect after that. Indeed the asymptotic total number of posts K is only 14,396. I don't believe this can be right, not even as a proxy, when the current count of monthly posts is well above 100.



I played around with data and the logistic growth function a little further, using annual instead of monthly data, changing the time horizon and fixing K, yet without much success.

Eventually I recalled a talk by Rob Hyndman's about his forecast package. After all, I have a time series here. So, applying the forecast function to the incremental data provides a somewhat more realistic prediction of 2,695 posts for the next 12 months, but with an increasing trend in monthly posts for 2014, which I find hard to believe given the observations over the last year.



Well, I presented two models here: One predicts a rapid decline in monthly posts on R-bloggers, while the other forecasts an increase. Neither feels right to me. Of course time will tell, but have you got any ideas or views?

Session Info

R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats  graphics  grDevices  utils  datasets  methods   base     

other attached packages:
[1] forecast_4.8 xts_0.9-7    zoo_1.7-10   XML_3.95-0.2

loaded via a namespace (and not attached):
 [1] colorspace_1.2-4        fracdiff_1.4-2         
 [3] grid_3.0.2              lattice_0.20-23        
 [5] nnet_7.3-7              parallel_3.0.2         
 [7] quadprog_1.5-5          Rcpp_0.10.6            
 [9] RcppArmadillo_0.3.920.1 tools_3.0.2            
[11] tseries_0.10-32

The Christmas and New Year's break is over, yet there is still time to return unwanted presents. Return to Santa was the title of an article in the Economist that highlighted the impact on online retailers, as return rates can be alarmingly high.

The article quotes a study by Christian Schulze of the Frankfurt School of Finance and Management, which analyses the return habits of customers who bought at least five items over a five year period from a large European online retailer. Although only a few figures are cited I will attempt to create a little model that replicates the customer behaviour and visualises the impact on overall profitability.

The study found that 5% of customers sent back more than 80% of the items they had bought; and that 1% of customers sent back at least 90% of their purchases. Or in other words 95% of customers send back less than 80% and 99% of customers send back less than 90%. To model this behaviour an S-shape curve seems appropriate, such as the logistic curve, as no-one can return more than they bought or less than nothing. With location and scale parameters m and s the logistic function can be fitted to the data, see the R code below.


The return rates do look quite high. However, if the products were shoes rather than books then I find them believable.

Additionally the article cites studies that suggest handling each returned item costs online sellers between $6 and $18, not to mention losses from items that are returned in unsaleable condition. Furthermore, without the cost of returns, online retailer's profits would be almost 50% higher.

Thus, to spin my toy model further, I assume 100 customers with revenues following an exponential distribution (λ=1/250), the cost ratio of sold goods to be lognormal (μ=-0.1, σ=0.1) and the cost of returns to follow a normal distribution with mean of $12 and standard deviation of $6.

In my simulation I could have made a profit of $1,979 instead of $1,441. Clearly the customers who return many items cause a real dent to my bottom line.

This situation is best visualised in what is often called a Whale Chart. Here I plot the cumulative profit against customers, with the most profitable customer on the left and the least profitable customer on the right. This chart shows me how much profit the first x number of customers generated. Often this graphs looks like a whale coming out of the water - hence its name.


In my little toy simulation I note that the first 20 most profitable customers would have generated more profit than the revenue of all customers. Indeed, profitability could have been 37% higher if it wasn't for loss making customers.

So, what shall I do? Manage my customers, know who I should reward and keep and whose loss wouldn't hurt at all. More customers are not the answer. I need more customers who return less.

R Code

Last week's Cologne R user group meeting was the best attended so far. Well, we had a great line up indeed. Matt Dowle came over from London to give an introduction to the data.table package. He was joined by his collaborator Arun Srinivasan, who is based in Cologne. Their talk was followed by Thomas Rahlf on Datendesign mit R (Data design with R).

data.table

Download slides

Matt's goal with the data.table package is to reduce times; time to write code and to execute code. His talk illustrated how the syntax of data.table, not unlike SQL, can produce shorter and more readable code that at the same time provides an efficient and fast way to analyse big in memory data sets with R. Arun presented on new developments in data.table 1.8.11, which not only fixes bugs but adds many new features such as melt/cast and further speed gains.

I said early that data.table rocks. For more details see the data.table home page.

Data design with R

Thomas Rahlf: Datendesign mit R

Thomas Rahlf talked about his forthcoming book Datendesign mit R (Data design with R). He shared with us his motivations and aims for the book. In his opinion there are many books that present beautiful charts and concepts (e.g. Tufte's books), but then don't show how they can be reproduced, as there are often done with software such as Adobe Illustrator. Or books explain the graphical functions of a software, yet fail to demonstrate how to create beautiful charts with them. Thus, Thomas' book will contain 100 examples demonstrating that desktop publishing quality charts can be produced with R and in some cases with the help of LaTeX. Indeed, all examples have about 40 lines of code and use the base R graphics system only and not grid or any add-ons such as lattice or ggplot2.

The book's accompanying web site gives you a taster already. The book itself will be published by Open Source Press next month.

The Schnitzel

Of course the evening ended with Schnitzel and Kölsch at the Lux.

The Luxus Schnitzel. Photo by Günter Faes

Next Kölner R meeting

The next meeting is scheduled for 26 February 2013 (Wednesday before Altweiber), with two talks by Diego de Castillo (Connecting R with databases) and Kim Kuen Tang (R and kdb+).

Please get in touch if you would like to present and share your experience, or indeed if you have a request for a topic you would like to hear more about. For more details see also our Meetup page.

Thanks again to Bernd Weiß for hosting the event and Revolution Analytics for their sponsorship.

Quick reminder: The next Cologne R user group meeting is scheduled for this Friday, 13 December 2013. We are delighted to welcome:

• Matt Dowle and Arun Srinivasan: Introduction to data.table
• Thomas Rahlf: Book presentation - Datendesign mit R

Further details and the agenda are available on our KölnRUG Meetup site.

Please sign up if you would like to come along. Notes from past meetings are available here.


The organisers, Bernd Weiß and Markus Gesmann, gratefully acknowledge the sponsorship of Revolution Analytics, who support the Cologne R user group as part of their vector programme.


View Larger Map

Following the very positive feedback that Andreas and I have received from delegates of the first R in Insurance conference in July of this year, we are planning to repeat the event next year. We have already reserved a bigger auditorium.

The second conference on R in Insurance will be held on Monday 14 July 2014 at Cass Business School in London, UK.

This one-day conference will focus again on applications in insurance and actuarial science that use R, the lingua franca for statistical computation. Topics covered may include actuarial statistics, capital modelling, pricing, reserving, reinsurance and extreme events, portfolio allocation, advanced risk tools, high-performance computing, econometrics and more. All topics will be discussed within the context of using R as a primary tool for insurance risk management, analysis and modelling.

The intended audience of the conference includes both academics and practitioners who are active or interested in the applications of R in insurance.

Invited talks will be given by:

• Arthur Charpentier, Département de mathématiques Université du Québec à Montréal
• Montserrat Guillen, Dept. Econometrics University of Barcelona together with Leo Guelman, Royal Bank of Canada (RBC Insurance division)

The members of the scientific committee are: Katrien Antonio (University of Amsterdam and KU Leuven), Christophe Dutang (Université du Maine, France), Jens Nielsen (Cass), Andreas Tsanakas (Cass) and Markus Gesmann (ChainLadder project).

Details about the registration and abstract submission process will be published soon on www.RinInsurance.com.

You can contact us via rinsuranceconference at gmail dot com.

The organisers, Andreas Tsanakas and Markus Gesmann, gratefully acknowledge the sponsorship of Mango Solutions, RStudio, Cybaea and PwC.