Archive
A brainfuck interpreter for R
The deadline for my book on R is fast approaching, so naturally I’m in full procrastination mode. So much so that I’ve spent this evening creating a brainfuck interpreter for R. brainfuck is a very simple programming language: you get an array of 30000 bytes, an index, and just 8 eight commands. You move the index left or right along the array with < and >; increase or decrease the value at the current position with + and -; read and write characters using . and ,; and start and end loops with [ and ].
There seem to be two approaches to creating a brainfuck interpreter: directly execute the commands, or generate code in a sensible language and execute that. I’ve opted for the latter approach because it’s easier, at least in R. Generating R code and then calling eval is probably a little slower than directly executing commands, but that’s the least of your worries with brainfuck. Even writing a trivial page-long program will take you many million times longer than it takes to execute.
The fact that you have to mix data variables (that 30000 element raw vector and an index) with commands means that an object oriented approach is useful. The whole interpreter is stored in a single reference class, of type brainfuck. Rather than me showing you all the code here, I suggest that you take a look at it (or clone it) from its repository on bitbucket. (I’ll submit to CRAN soon.)
Here’s a Hello World example taken from Wikipedia. To use the brainfuck package, you just create/import your brainfuck program as a character vector (non-command characters are ignored, so you can comment your code). Call fuckbrain once to create the interpreter variable, then call its interpret method on each program that you want to run.
library(brainfuck)
hello_world <- "+++++ +++++ initialize counter (cell #0) to 10
[ use loop to set the next four cells to 70/100/30/10
> +++++ ++ add 7 to cell #1
> +++++ +++++ add 10 to cell #2
> +++ add 3 to cell #3
> + add 1 to cell #4
<<<< - decrement counter (cell #0)
]
> ++ . print 'H'
> + . print 'e'
+++++ ++ . print 'l'
. print 'l'
+++ . print 'o'
> ++ . print ' '
<< +++++ +++++ +++++ . print 'W'
> . print 'o'
+++ . print 'r'
----- - . print 'l'
----- --- . print 'd'
> + . print '!'
> . print '\n'"
bfi <- fuckbrain()
bfi$interpret()
A little Christmas Present for you
Here’s an excerpt from my chapter “Blood, sweat and urine” from The Bad Data Handbook. Have a lovely Christmas!
I spent six years working in the statistical modeling team at the UK’s Health and Safety
Laboratory. A large part of my job was working with the laboratory’s chemists, looking
at occupational exposure to various nasty substances to see if an industry was adhering
to safe limits. The laboratory gets sent tens of thousands of blood and urine samples
each year (and sometimes more exotic fluids like sweat or saliva), and has its own team
of occupational hygienists who visit companies and collect yet more samples.
The sample collection process is known as “biological monitoring.” This is because when
the occupational hygienists get home and their partners ask “How was your day?,” “I’ve
been biological monitoring, darling” is more respectable to say than “I spent all day
getting welders to wee into a vial.”
In 2010, I was lucky enough to be given a job swap with James, one of the chemists.
James’s parlour trick is that, after running many thousands of samples, he can tell the
level of creatinine in someone’s urine with uncanny accuracy, just by looking at it. This
skill was only revealed to me after we’d spent an hour playing “guess the creatinine level”
and James had suggested that “we make it more interesting.” I’d lost two packets of fig
rolls before I twigged that I was onto a loser.The principle of the job swap was that I would spend a week in the lab assisting with
the experiments, and then James would come to my office to help out generating the
statistics. In the process, we’d both learn about each other’s working practices and find
ways to make future projects more efficient.
In the laboratory, I learned how to pipette (harder than it looks), and about the methods
used to ensure that the numbers spat out of the mass spectrometer4 were correct. So as
well as testing urine samples, within each experiment you need to test blanks (distilled
water, used to clean out the pipes, and also to check that you are correctly measuring
zero), calibrators (samples of a known concentration for calibrating the instrument5),
and quality controllers (samples with a concentration in a known range, to make sure
the calibration hasn’t drifted). On top of this, each instrument needs regular maintaining
and recalibrating to ensure its accuracy.
Just knowing that these things have to be done to get sensible answers out of the ma?
chinery was a small revelation. Before I’d gone into the job swap, I didn’t really think
about where my data came from; that was someone else’s problem. From my point of
view, if the numbers looked wrong (extreme outliers, or otherwise dubious values) they
were a mistake; otherwise they were simply “right.” Afterwards, my view is more
nuanced. Now all the numbers look like, maybe not quite a guess, but certainly only an
approximation of the truth. This measurement error is important to remember, though
for health and safety purposes, there’s a nice feature. Values can be out by an order of
magnitude at the extreme low end for some tests, but we don’t need to worry so much
about that. It’s the high exposures that cause health problems, and measurement error
is much smaller at the top end.
Have my old job!
My old job at the Health & Safety Laboratory is being advertised, and at a higher pay grade to boot. (Though it is still civil service pay, and thus not going to make you rich.)
You’ll need to have solid mathematical modelling skills, particularly solving systems of ODEs, and be proficient at writing scientific code, preferably R or MATLAB or acslX. From chats with a few people at the lab, management are especially keen to get someone who can bring in money so grant writing and blagging skills are important too.
It’s a smashing place to work and the people are lovely. Also, you get flexitime and loads of holiday. If you are looking for a maths job in North West* England then I can heartily recommend applying.
*Buxton is sometimes North West England (when we get BBC local news) and sometimes in the East Midlands (like when we vote in European elections).
Indexing with factors
This is a silly problem that bit me again recently. It’s an elementary mistake that I’ve somehow repeatedly failed to learn to avoid in eight years of R coding. Here’s an example to demonstrate.
Suppose we create a data frame with a categorical column, in this case the heights of ten adults along with their gender.
(heights <- data.frame(
height_cm = c(153, 181, 150, 172, 165, 149, 174, 169, 198, 163),
gender = c("female", "male", "female", "male", "male", "female", "female", "male", "male", "female")
))
Using a factory fresh copy of R, the gender column will be assigned a factor with two levels: “female” and then “male”. This is all well and good, though the column can be kept as characters by setting stringsAsFactors = FALSE.
Now suppose that we want to assign a body weight to these people, based upon a gender average.
avg_body_weight_kg <- c(male = 78, female = 63)
Pop quiz: what does this next line of code give us?
avg_body_weight_kg[heights$gender]
Well, the first value of heights$gender is “female”, so the first value should be 63, and the second value of heights$gender is “male”, so the second value should be 78, and so on. Let’s try it.
avg_body_weight_kg[heights$gender] # male female male female female male male female female male # 78 63 78 63 63 78 78 63 63 78
Uh-oh, the values are reversed. So what really happened? When you use a factor as an index, R silently converts it to an integer vector. That means that the first index of “female” is converted to 1, giving a value of 78, and so on.
The fundamental problem is that there are two natural interpretations of a factor index – character indexing or integer indexing. Since these can give conflicting results, ideally R would provide a warning when you use a factor index. Until such a change gets implemented, I suggest that best practice is to always explicitly convert factors to integer or to character before you use them in an index.
avg_body_weight_kg[as.character(heights$gender)] avg_body_weight_kg[as.integer(heights$gender)]
Make your data famous!
I’m writing a book on R for O’Reilly, and I need interesting datasets for the examples. Any data that you provide will get you a mention in the book and in the publicity material, so it’s a great opportunity to publicise your work or your organisation.
Datasets from any area or industry are suitable; the only constraint is that it can be analysed with a few pages of R code to provide a result that a general reader might go “ooh”. There’s a chapter on data cleaning, so even dirty data is suitable!
All the data will be provided in an R package to accompany the book, so you need to be willing to make it publically available. I can help you anonymise the data, or strip out commercially sensitive parts if you require.
If you can provide anything, or you know someone who might be able to, then drop me an email at richierocks AT gmail DOT com. Thanks.
EDIT: There are some (quite) frequently asked questions already! Here are the answers; you can use your Jeopardy! skills to guess the questions.
1. The book is called “Learning R”, and it’s a fairly gentle introduction to the language, covering both how you program in R, and how you analyse data.
2. If you provide data, then yes, you can have an PDF of the pre-release version to make sure I haven’t done something silly with your dataset.
How long does it take to get pregnant?
My girlfriend’s biological clock is ticking, and so we’ve started trying to spawn. Since I’m impatient, that has naturally lead to questions like “how long will it take?”. If I were to believe everything on TV, the answer would be easy: have unprotected sex once and pregnancy is guaranteed.
A more cynical me suggests that this isn’t the case. Unfortunately, it is surpisingly difficult to find out the monthly chance of getting pregnant (technical jargon: the “monthly fecundity rate”, or MFR), given that you are having regular sex in the days leading up to ovulation. Everyone agrees that age has a big effect, with women’s peak fertility occuring somewhere around the age of 25. Beyond that point, the internet is filled with near-useless summary statistics like the chance of conceiving after one year. For example, the usually reliable NHS site says
Women become less fertile as they get older. For women aged 35, about 94 out of every 100
who have regular unprotected sex will get pregnant after three years of trying. However, for
women aged 38, only 77 out of every 100 will do so.
I found a couple of reasonably sciency links(George and Kamath, Socal Fertility) that suggest that the MFR is about 25% for a women aged 25, and 10% at age 35. The Scoal link also gives rates of 15% at age 30, 5% at age 40 and less than 1% at age 45. If the woman is too fat, too thin, a smoker, or has hormone problems, or is stressed, then the rate needs reducing.
Given the MFR, the probability of getting pregnant after a given number of months can be calculated with a negative binomial distribution.
months <- 0:60
p_preg_per_month <- c("25" = 0.25, "30" = 0.15, "35" = 0.1, "40" = 0.05, "45" = 0.01)
p_success <- unlist(lapply(
p_preg_per_month,
function(p) pnbinom(months, 1, p)
))
Now we just create a data frame suitable for passing to ggplot2 …
mfr_group <- paste( "MFR =", format(p_preg_per_month, digits = 2), "at age", names(p_preg_per_month) ) mfr_group <- factor(mfr_group, levels = mfr_group) preg_data <- data.frame( months = rep.int(months, length(mfr_group)) , mfr_group = rep(mfr_group, each = length(months)), p_success = p_success )
and draw the plot.
library(ggplot2)
(p <- ggplot(preg_data, aes(months, p_success, colour = mfr_group)) +
geom_point() +
scale_x_continuous(breaks = seq.int(0, 60, 12)) +
scale_y_continuous(breaks = seq.int(0, 1, 0.1), limits = c(0, 1)) +
scale_colour_discrete("Monthly fecundity rate") +
xlab("Months") +
ylab("Probability of conception") +
opts(panel.grid.major = theme_line(colour = "grey60"))
)
So almost half of the (healthy) 25 year olds get pregnant in the first monthtwo months, and after two years (the point when doctors start considering you to have fertility problems) more than 90% of 35 year olds should conceive. By contrast, just over 20% of 45 year old women will. In fact, even this statistic is over-optimistic: at this age, fertility is rapidly decreasing, and a 1% MFR at age 45 will mean a much lower MFR at age 47 and the negative binomial model breaks down.
Of course, from a male point of view, conception is an embarrassingly parallel problem: you can dramatically reduce the time to conceive a child by sleeping with lots of women at once. (DISCLAIMER: Janette, if you’re reading this, I’m not practising or advocating this technique!)
Be assertive!
assertive, my new package for writing robust code, is now on CRAN. It consists of lots of is functions for checking variables, and corresponding assert functions that throw an error if the condition doesn’t hold. For example, is_a_number checks that the input is numeric and scalar.
is_a_number(1) #TRUE
is_a_number("a") #FALSE
is_a_number(1:10) #FALSE
In the last two cases, the return value of FALSE has an attribute “cause” that indicates the cause of failure. When “a” is the input, the cause is “"a" is not of type 'numeric'.“, whereas for 1:10, the cause is “1:10 does not have length one.“. You can get or set the cause attribute with the cause function.
m <- lm(uptake ~ 1, CO2) ok <- is_empty_model(m) if(!ok) cause(ok)
The assert functions call an is function, and if the result is FALSE, they throw an error; otherwise they do nothing.
assert_is_a_number(1) #OK
assert_is_a_number("a") #Throws an error
There are also some has functions, primarily for checking the presence of attributes.
has_names(c(foo = 1, bar = 4, baz = 9)) has_dims(matrix(1:12, nrow = 3))
Some functions apply to properties of vectors. In this case, the assert functions can check that all the values conform to the condition, or any of the values conform.
x <- -2:2 is_positive(x) #The last two are TRUE assert_any_are_positive(x) #OK assert_all_are_positive(x) #Error
“Why would you want to use these functions?”, you may be asking. The dynamic typing and extreme flexibility of R means that it is very easy to have variables that are the wrong format. This is particularly true when you are dealing with user input. So while you know that the sales totals passed to your function should be a vector of non-negative numbers, or that the regular expression should be a single string rather than a character vector, your user may not. You need to check for these invalid conditions, and return an error message that the user can understand. assertive makes it easy to do all this.
Since this is the first public release of assertive, it hasn’t been widely tested. I’ve written a moderately comprehensive unit-test suite, but there are likely to be a few minor bugs here and there. In particular, I suspect there may be one or two typos in the documentation. Please give the package a try, and let me know if you find any errors, or if you want any other functions adding.
Benford’s Law and fraud in the Russian election
Earlier today Ben Goldacre posted about using Benford’s Law to try and detect fraud in the Russian elections. Read that now, or the rest of this post won’t make sense. This is a loose R translation of Ben’s Stata code.
The data is held in a Google doc. While it is possible to directly retrieve the contents with R, for a single document it is easier to save it a CSV, and load it from your own machine.
russian <- read.csv("Russian observed results - FullData.csv")
There are loads of ways of manipulating data and plotting it in R, and while you can do everything in the base R distribution, I’m going to use a few packages to make it easier.
library(reshape) library(stringr) library(ggplot2)
A little transformation is needed. We take only the columns containing the counts and manipulate the data into a “long” format with only one value per row.
russian <- melt(
russian[, c("Zhirinovsky", "Zyuganov", "Mironov", "Prokhorov", "Putin")],
variable_name = "candidate"
)
Now we add columns containing the first and last digits, extracted using regular expressions.
russian <- ddply(
russian,
.(candidate),
transform,
first.digit = str_extract(value, "[123456789]"),
last.digit = str_extract(value, "[[:digit:]]$"))
The table function gives us the counts of each number, and we compare these against the counts predicted by Benford’s Law.
first_digit_counts <- as.vector(table(russian$first.digit)) first_digit_actual_vs_expected <- data.frame( digit = 1:9, actual.count = first_digit_counts, actual.fraction = first_digit_counts / nrow(russian), benford.fraction = log10(1 + 1 / (1:9)) )
The counts of the last digit can be obtained in a similar way.
last_digit_counts <- as.vector(table(russian$last.digit))
last_digit_actual_vs_expected <- data.frame(
digit = 0:9,
count = last_digit_counts,
fraction = last_digit_counts / nrow(russian)
)
last_digit_actual_vs_expected$cumulative.fraction <- cumsum(last_digit_actual_vs_expected$fraction)
Here is the line graph…
a_vs_e <- melt(first_digit_actual_vs_expected[, c("digit", "actual.fraction", "benford.fraction")], id.var = "digit")
(fig1_lines <- ggplot(a_vs_e, aes(digit, value, colour = variable)) +
geom_line() +
scale_x_continuous(breaks = 1:9) +
scale_y_continuous(formatter = "percent") +
ylab("Counts with this first digit") +
opts(legend.position = "none")
)
and the histogram
(fig2_hist <- ggplot(russian, aes(value)) +
geom_histogram(binwidth = 20)
)
Radical Statistics was radical
Today I went to the Radical Statistics conference in London. RadStats was originally a sort of left wing revolutionary group for statisticians, but these days the emphasis is on exposing dubious statistics by companies and politicians.
Here’s a quick rundown of the day.
First up Roy Carr-Hill spoke about the problems with trying to collect demographic data and estimating soft measures of societal progress like wellbeing. (Household surveys exclude people not in households, like the homeless soldiers and old people in care homes; and English people claim to be 70% satisfied regardless of the question.)
Next was Val Saunders who started with a useful debunking of done methodological flaws in schizophrenia research, then blew it by detailing her own methodologically flaws research and making overly strong claims to have found the cause of that disease.
Aubrey Blunsohn and David Healy both talked about ways that the pharmaceutical industry fudges results. The list was impressively long, leading me to suspect that far to many people have spent far too long thinking of ways to game the system. The two main recommendations that resonated with me were to extend the trials register to phase 1 trials to avoid unfavourable studies being buried and for raw data to be made available for transparent analysis. Pipe dreams.
After lunch Prem Sikka pointed out that tax avoidance isn’t just shady companies trying to scam the system, but actually accountancy firms pay people to dream up new wheezes and sell them to those companies.
Ann Pettifor and final speaker Howard Reed had similar talks evangelising Keynesian stimulus (roughly, big government spending in times of recession) for the UK economy amongst some economic myth debunking. Thought provoking, though both speakers neglected to mention the limitations of such stimuli – you have to avoid spending in pork barrel nonsense (see Japan in the 90s, that buy-a-banger scheme in the UK in 2009) and you have to find a ways to turn of the taps w when recession is over.
The other speaker was Allyson Pollack who discussed debunking a dubious study by Zac Cooper claiming that patients being allowed to choose their surgeon improved success rated treating acute myocardial infarction. Such patients are generally unconscious while having their heart attack so out was inevitably nonsense.
Overall a great day.
Richie Cotton
Licensing
The non-code parts of 4D Pie Charts by Richard Cotton are licensed under a Creative Commons Attribution-NoDerivs 2.0 UK: England & Wales License. The code parts of the blog are licensed under the WTFPL v2.0.
