Keep on Fighting!

Here are some (trivial) R tips in the course Stat 511. I’ll update this post till the semester is over.

1. Formatting R Code

Reading code is pain, but the well-formatted code might alleviate the pain a little bit. The function tidy.source() in the animation package can help us format our R code automatically. By default it will read your code in the clipboard, parse it and return the well-formatted code. You have options to keep or remove the comments/blank lines and set the width of the code, etc. Spaces and indent will be added automatically. This can save us time typing spaces and paying attention to indent.

## install.packages('animation') if it is not installed yet
library(animation)
## copy some R code somewhere and type:
tidy.source()
## or specify the path of your code file
tidy.source(file.path(system.file(package = "graphics"), "demo", "image.R"))
## can also use a URL
tidy.source('http://www.public.iastate.edu/~dnett/S511/twofactor.R')
## remove blank lines
tidy.source('http://www.public.iastate.edu/~dnett/S511/twofactor.R',
           keep.blank.line = FALSE)
## remove comments
tidy.source('http://www.public.iastate.edu/~dnett/S511/twofactor.R',
           keep.comment = FALSE)

Continue reading »

I have to admit that the previous post on Christmas is actually not much fun. Today I received another pResent from Yixuan which is more interesting:

Life should be fun. I saw a post in R-help list saying Merry Christmas to other useRs, and I followed up by some R code which can produce a naive animation like this:

Here is the code to generate the above Flash animation with shining Christmas:

library(animation)
saveSWF({
    n = length(speed <- runif(angle <- runif(x <- strsplit("MERRY CHRISTMAS",
        "")[[1]], 0, 360), 0, 15))
    for (j in 1:300) {
        angle = angle + speed
        plot.new()
        plot.window(c(1, n), c(0, 1))
        for (i in 1:n) text(i, 0.5, x[i], srt = angle[i], cex = runif(1,
            1, 4), col = sample(colors(), 1))
        text(n, 0, "Yihui @ 2009-12-24 (http://yihui.name)",
            adj = c(1, 0), col = "white", cex = 0.8)
    }
}, interval = 0.04, dev = "pdf", outdir = getwd(), para = list(mar = rep(0,
    4), bg = "black"), width = 8, height = 1)
## in animation package (>=1.1-0), see demo('Xmas')

There are other animation formats in the R package animation:

1. use saveMovie() to get a GIF animation (need ImageMagick)
2. ani.start() and ani.stop() can produce an HTML page with the animation in it
3. saveLatex() can embed an animation into a PDF document

Since animation 1.0-9, we will be able to create a PDF document with an animation embedded in it; the function is saveLatex(), and its usage is similar to saveMovie() and saveSWF(): you pass an R expression for creating animations to this function, and this expression will be evaluated in the function; the image frames get recorded by a graphics device. In the end, a LaTeX document is written in a directory, and we can get a PDF document by running pdflatex on the document.

In fact, the key point is the LaTeX package named animate, which can be used to insert image frames into a PDF document to generate an animation. The interface of animations created by this package is quite similar to the HTML animation page by the R package animation, moreover, it also uses JavaScript (in PDF) to animate the image frames.

Today Romain Francois posted an interesting topic in the R-help list, and you can read his blog post for more details: celebrating R commit #50000. 50000 is certainly not a small number; we do owe R core members a big “thank you” for their great efforts in this fantastic statistical language in the 13 years. When I saw Romain’s data, I suddenly remembered a question I asked to one of Prof Ripley’s student a couple of years ago: does Prof Ripley ever sleep? And he answered “No!”. No wonder we can see Prof Ripley so frequently in the R-help/devel mailing list. If you have stayed on R-help list for enough long time, you’ll surely know several facts, e.g. Martin Maechler will arrive in less than 3 minutes if you dare call an R package “library”, and you will get “Ripleyed” if you are not careful enough in posting your R code.

> library(fortunes)
> fortune("Ripleyed")

And the fear of getting Ripleyed on the mailing list also makes me think, read,
and improve before submitting half baked questions to the list.
 -- Eric Kort
 R-help (January 2006)

Linlin Yan posted a cool (hot?) simulation of burning fire with R in the COS forum yesterday, which was indeed a warm welcome. I’m not sure whether our forum members will be scared by the “fire” under the title “Welcome to COS Forum”. The fire was mainly created by the function image() with carefully designed rows and columns in heated colors heat.colors(). Here is one of the pictures generated from his code:

Today Ruya Gokhan Kocer asked me how to use the R function identify() in off-screen graphics devices. Actually it’s pretty easy as long as we obtain the list returned by identify(pos = TRUE). For example,

Downdload the file here

# open a windows device
x11()
x = rnorm(20)
y = rnorm(20)
plot(x, y)
# identify 5 points
id = identify(x, y, n = 5, pos = TRUE)

# $ind
# [1]  2  6 10 14 16
#
# $pos
# [1] 1 1 4 4 1

# then open a bitmap device
png("identify.png")
plot(x, y)
# use the information from above mouse click
text(x[id$ind], y[id$ind], id$ind, pos = id$pos)
dev.off()

I was surprised to find the density estimation of a constant was also “bell-shaped” by default when a friend passed some R code to me to illustrate CLT, but I realized the reason soon.

Downdload the R code here

# png(width = 500, height = 300)
x = rep(0, 1000)
par(mfrow = c(1, 2), mar = c(4, 4, 0.1, 0.1))
plot(density(x), main = "")
plot(density(x), main = "")
rug(jitter(x))
# dev.off()

Note that I added a rug (jittered) to the right plot to tell you the true locations of the data points.

Yesterday I wrote some R code to simulate the quincunx, which looks like:

It will appear as a function in the R package animation soon.

Here is an example of dynamically selecting points using R (the function getGraphicsEvent()):

Downdload the R code here

par(bg = "black", mar = rep(0, 4), pch = 20)
xx = runif(100)
yy = runif(100)
plot(xx, yy, type = "n")
mousemove = function(buttons, x, y) {
    r = 0.2
    idx = (x - r < xx & xx < x + r) & (y - r < yy & yy < y + r)
    plot(xx, yy, type = "n")
    rect(x - r, y - r, x + r, y + r, border = "yellow", lty = 2)
    points(xx[idx], yy[idx], col = "yellow", cex = 2)
    points(xx[!idx], yy[!idx], col = "red")
    NULL
}
mousedown = function(buttons, x, y) {
    "Done"
}
getGraphicsEvent("Click mouse to exit", onMouseDown = mousedown,
    onMouseMove = mousemove)

We can adjust the parameter r as we wish.