Example: Logistic distribution

I teach this algorithm in one of my classes and I’m always on the look-out for new examples. Something that escaped my notice is that it is easy to generate RN’s using this technique from the Logistic distribution. This distribution has CDF

and so we can generate a random number from the logistic distribution using the following formula:

November 27, 2010

I’ve just received my copy of Advanced Markov Chain Monte Carlo Methods, by Liang, Liu, & Carroll. Although my PhD didn’t really involve any Bayesian methodology (and my undergrad was devoid of any Bayesian influence), I’ve found that the sort of problems I’m now tackling in systems biology demand a Bayesian/MCMC approach. There are a number of introductory books on MCMC, but not that many on advanced techniques.

This book suggests that it could be used as a possible textbook or reference guide in a one-semester statistics graduate course. I’m not entirely convinced that it would be a good textbook, but as a reference it looks very promising. A word of warning, if you’re not familiar with MCMC then you should try the Dani Gamerman MCMC book first. Some later chapters look particularly interesting, including topics on adaptive proposals, population-based MCMC methods and dynamic weightings.

Anyway, I intend to work through the book (well at least a few of the chapters) and post my results/code as I go. Well that’s the plan anyway.

November 23, 2010

Each year I have the pleasure (actually it’s quite fun) of teaching R programming to first year mathematics and statistics students. The vast majority of these students have no experience of programming, yet think they are good with computers because they use facebook!

Debugging students' R scripts

The class has around 100 students, and there are eight practicals. In some of these practicals the students have to submit code. Although the code is “marked” by a script, this only detects if the code is correct. Therefore, I have to go through a lot of R functions by hand and find bugs.

First year the course ran, I had no style guide.

Result: spaghetti R code.

Second year: asked the students to indent their code.

In fact, during practicals I refused to debug in any R code that hadn’t been indented.

Result: nicer looking code and more correct code.

This year I intend to introduce a R style guide based loosely on Google’s and Hadley’s guides.

One point that’s in my guide and not (and shouldn’t be) in the above style guides, is that all functions must have one and only return statement. I tend to follow the single return rule for the majority of my R functions, but do, on occasions, break it. The bible of code styling, Code Complete, recommends that you use returns judiciously.

R Style Guide

This style guide is intended to be very light touch. It’s intended to give students the basis of good programming style, not be a guide for submitting to cran.

File names

File names should end in .R and, of course, be meaningful. Files should be stored in a meaningful directory – not your Desktop!
GOOD: predict_ad_revenue.R
BAD: foo.R

Variable & Function Names

Variable names should be lowercase. Use _ to separate words within a name. Strive for concise but meaningful names (this is not easy!)
GOOD: no_of_rolls
BAD: noOfRolls, free

November 16, 2010

In R, you can use both ‘=’ and ‘<-‘ as assignment operators. So what’s the difference between them and which one should you use?

What’s the difference?

The main difference between the two assignment operators is scope. It’s easiest to see the difference with an example:
##Delete x (if it exists)
> rm(x)
> mean(x=1:10) #[1] 5.5
> x #Error: object 'x' not found

Here x is declared within the function’s scope of the function, so it doesn’t exist in the user workspace. Now, let’s run the same piece of code with using the <- operator:
> mean(x <- 1:10)# [1] 5.5
> x # [1] 1 2 3 4 5 6 7 8 9 10

This time the x variable is declared within the user workspace.

When does the assignment take place?

In the code above, you may be tempted to thing that we “assign 1:10 to x, then calculate the mean.” This would be true for languages such as C, but it isn’t true in R. Consider the following function:
> a <- 1
> f <- function(a) return(TRUE)
> f <- f(a <- a + 1); a
[1] TRUE
[1] 1

Notice that the value of a hasn’t changed! In R, the value of a will only change if we need to evaluate the argument in the function. This can lead to unpredictable behaviour:
> f <- function(a) if(runif(1)>0.5) TRUE else a
> f(a <- a+1);a
[1] 2
> f(a <- a+1);a
[1] TRUE
[1] 2
> f(a <- a+1);a
[1] 3

November 13, 2010

Stochastic simulation is a very important tool for Systems Biology modelling. However, it is difficult to check the correctness of a stochastic simulator, since any two realizations from a single model will typically be different. Recognising this, we developed the discrete stochastic model test suite. The test suite contains four core models, and each model has a number of variants. We restricted our attention to models which admit direct analytic or numerical analysis. There are not many of these, and they are all fairly simple.

Today I updated the test suite with a new version. Only minor changes were made to the test suite (thanks to Sean Mauch for pointing them out). Here’s a list of changes:

November 6, 2010

Part of the reason R has become so popular is the vast array of packages available at the cran and bioconductor repositories. In the last few years, the number of packages has grown exponentially!

This is a short post giving steps on how to actually install R packages. Let’s suppose you want to install the ggplot2 package. Well nothing could be easier. We just fire up an R shell and type:
> install.packages("ggplot2")

In theory the package should just install, however:

if you are using Linux and don’t have root access, this command may not work. If you use Ubuntu then this isn’t a problem.

if you have limited space in /home/ then you may want to install the package in the non-default directory.

you will be asked to select your local mirror, i.e. which server should you use to download the package.

if you have a proxy, then you may run into problems.

Installing packages in a different directory

First, you need to designate a directory where you will store the downloaded packages. On my machine, I use the directory /data/Rpackages/ After creating a package directory, to install a package we use the command:
> install.packages("ggplot2", lib="/data/Rpackages/")
> library(ggplot2, lib.loc="/data/Rpackages/")

It’s a bit of a pain having to type /data/Rpackages/ all the time. To avoid this burden, we create a file .Renviron in our home area, and add the line R_LIBS=/data/Rpackages/ to it. This means that whenever you start R, the directory /data/Rpackages/ is added to the list of places to look for R packages and so:

> install.packages("ggplot2")
> library(ggplot2)

just works!

Setting the repository

Every time you install a R package, you are asked which repository R should use. To set the repository and avoid having to specify this at every package install, simply:

November 1, 2010

In the last six months or so, the behemoth of Q & A sites stackoverflow, decided to change tack and launch a number of other non-computing-language sites. To launch a site in the stackoverflow family, sites have to spend time gathering followers in Area51. Once a site has gained a critical mass, a new StackExchange (SE) site is born.

At present there are around twenty-one SE beta sites. Being rather bored this weekend, I decided to see how these sites are similar/different. For a first pass, I did something rather basic, but useful none the less.

First, we need to use the stackoverflow api to download some summary statistics for each site: