统计图示
R语言作图
昨天Ripley教授向R提交了第50000次修改，Romain Francios对SVN的日志数据做了一些简单分析，我个人一直关心Ripley是不是整天不用睡觉（你看这老爷子一天到晚都在邮件列表中出没），这次正好验证一下他是不是24小时工作，数据和R代码参见50000 Revisions Committed to R。
Ripley从1999年加入R核心团队，从上图可以看出，他显然是不需要睡觉的——每个小时都可能有commit。时间分布呈双峰：早上7点到10点、下午3点到6点。看看其他人的工作时间，很容易发现Martin倾向于早起干活，而Peter倾向于每天晚上12点之后干活。
Romain的博客中有SVN的日志数据可以下载，感兴趣的同志们可以继续分析R core的工作行为。
上周上课，我们老爷子又说SAS is extremely powerful，SAS很靠谱，就差明说SAS没有Bug了。我正在整理课程笔记，SAS的事情，我改天要去找老爷子好好谈一谈。R从不说自己能担保什么，大家拼命找bug，拼命改进，这是开源软件的共同特征——没有人付钱，但就是有一群疯子半夜3点还在写代码。商业软件一向说自己能保证什么，可是一个bug二十年都没人会去修正（例），还好我们看不见源代码，要是能看见，后果不堪设想。
前些日子有位童鞋在R-help邮件列表里问如何生成那种单词大小与其出现频率成比例的图，这玩意儿也就是通常所说的标签云（Tag Cloud）。我琢磨了一下WordPress的插件wp-cumulus，发现其原理很简单，不过就是将标签信息以XML形式通过JavaScript传递给一个Flash文件，所以也很容易用R去实现这个传递过程，即：将文本、超级链接以及频数写成XML，然后嵌入到HTML文件中。整个过程参见Creating Tag Cloud Using R and Flash / JavaScript (SWFObject)这篇日志，函数源代码和示例数据都可以从那里下载。
效果是这样的：
（通过RSS阅读的童鞋们请打开原文链接在浏览器中观看，否则啥都看不到）
x区间上取一系列点，然后计算y = f(x)，把这些点描在笛卡尔坐标系中，然后用铅笔手工连起来。这个办法很笨很原始，但很有效，奇怪的是，到计算机如此发达的今天，仍然有人不断问如何用软件画曲线，仿佛中学时学的东西都还给数学老师了。例如，若干天前，有人问我如何用R画这个曲线：
我只好提醒他回忆中学的描点法。R里面的gamma()和curve()函数是现成的，把数学公式写成相应的代码就可以了：
curve((sqrt(2) * gamma(x/2) * 2)/(sqrt(x - 1) * gamma(x - 1)), 2, 50) # 当然，如果你非要按照x = seq(2, 50, ...)然后计算f(x)然后plot(x, f(x), type = "l") # 那也可以，只不过curve()就是干这事的，不必麻烦分两步走了
描点法画函数曲线
如同某些网文说人一辈子的守则在幼儿园都已经教完了，就看你长大还记不记得。所以每当面对一个貌似复杂的问题时，首先要想想幼儿园有没有学过。
“武功再高，也怕菜刀”绝对是真理，哈哈。
此时此刻，俺哆嗦着告诉各位客官，俺真的中奖了。
话说两个月前，我介绍了一下John Chambers软件奖，当时并没有指望获奖，只是打探打探情况，搞明白之后等我明年好好琢磨一个靠谱的软件之后再正式参加。结果，不知是美国人民都忙着对付金融危机了，还是一共就我一个人参加，或是怎么的，今日打开邮箱，发现我被通知得奖（中奖？）了。我愣了半天，觉得这是在做梦。
好吧，两个月前说了要请客，那么那篇日志底下的观众朋友们请赶紧拿出小本本记下，“谢某人欠俺一顿大餐”。
做动画都能获奖……敢问谱在何方啊……管它呢，反正是美国人的钱，俺遵旨，8月去华盛顿领钱，顺便领奖。
受Zhang H.启发，自个儿看了看Gmail邮箱的联系人，导出为Outlook的CSV格式，然后读进R，用正则表达式去掉*@部分，table()了一下。Hotmail邮箱的频数比我想象的要高，然后也没想到中科院植物所也排在前面，估计是那次R会议所致。剩下的联系人就分布在五湖四海了，结果请看：
频数文件：
x = read.csv("contacts.csv", stringsAsFactor = FALSE)
y = gsub("^.*@", "", x$E.mail.Address)
y = sort(table(y))
我们可以把R-help邮件列表当作一个练手的地方，比如今天有人想把地图上的指北箭头作为位图格式导入他的地图，而我则建议不如手工画一个：
north.arrow = function(x, y, h) {
polygon(c(x, x, x + h/2), c(y - h, y, y - (1 + sqrt(3)/2) *
h), col = "black", border = NA)
polygon(c(x, x + h/2, x, x - h/2), c(y - h, y - (1 + sqrt(3)/2) *
h, y, y - (1 + sqrt(3)/2) * h))
text(x, y, "N", adj = c(0.5, 0), cex = 4)
}
par(mar = rep(0, 4))
plot(1, type = "n", ylim = c(0, 1), axes = FALSE,
ann = FALSE)
north.arrow(1, 0.8, 0.4)
不过这种基础性练习带来的坏处就是，你总觉得用R的基础图形元素什么图都可以画，从而不愿学习新的作图系统。
用地图展示统计信息有一个明显的缺陷，就是各地区的面积大小可能会极大地影响你的视觉。如：北京上海面积小，但人口多，若是直接用地图的颜色或斜线密度表示人口多少，那么这些面积小的地区将会很不显眼，而事实上它们本应该是很“吸引眼球”的地方，这种情况下，我们可以将地区边界变形，使得面积和统计信息成比例。尽管你会看到一副很丑的图，但表达的信息却是被校正过的。
一个更明显的例子是美国的大选，若用红蓝标示各州的选举情况，似乎红方麦凯恩大胜了，而事实是麦凯恩获胜的州都是人口稀少的地方，总账算下来，仍然是奥巴马那小子赢了（比如他拿下了人口众多的加州）。看看以下两幅地图的对比：
Duncan Temple Lang几年前写了一个包叫Rcartogram，发在Omegahat上，不过不幸的是，Windows用户不能使用。它是基于一个C程序的，而这个程序里面的算法又涉及到另一个做快速傅里叶变换（FFT）的软件包，因此在Windows下编译起来很复杂（我花了很长时间也没搞清楚）。
Cartogram的一种算法可以参见”Diffusion-based method for producing density-equalizing maps“一文。
最近两天Friedrich在R-devel发了一个通知，告诉大家要开始准备Google Summer of Code了！然后一位叫Oleg Sklyar的同志迫不及待提出了两个主意，结果讨论很快就乱了套，因为他提的第一个主意就是关于交互式图形的，他说很关键的一点是R不能实现图形的缩放和平移（zooming and panning），我一看，不对劲呀，谁说不能缩放、平移呀，于是三下五除二给了一个例子：
##################################################################
# a demo for zooming and panning in R graphics
# by Yihui Xie, Feb 20, 2009
##################################################################
# a large number of points
plot(x <- rnorm(5000), y <- rnorm(5000), xlab = "x", ylab = "y")
xylim <- c(range(x), range(y))
zoom <- function(d, speed = 0.05) {
rx <- speed * (xylim[2] - xylim[1])
ry <- speed * (xylim[4] - xylim[3])
# global assignment '<<-' here!
xylim <<- xylim + d * c(rx, -rx, ry, -ry)
plot(x, y, xlim = xylim[1:2], ylim = xylim[3:4])
NULL
}
# Key `+`: zoom in; `-`: zoom out
# Left, Right, Up, Down: self-explaining
# `*`: reset
# Press other keys to quit
keybd <- function(key) {
switch(key, `+` = zoom(1), `-` = zoom(-1), Left = zoom(c(-1,
1, 0, 0)), Right = zoom(c(1, -1, 0, 0)), Up = zoom(c(0,
0, 1, -1)), Down = zoom(c(0, 0, -1, 1)), `*` = plot(x,
y), "Quit the program")
}
getGraphicsEvent(onKeybd = keybd)
##################################################################
上图仅限于Windows下操作，因为getGraphicsEvent()函数仅限于Windows图形设备使用。通过键盘上的加（放大）减（缩小）乘（恢复）和上下左右键可以很容易实现图形的动态缩放和平移。
不过后来Simon Urbanek老大发话了，Yihui你这不叫交互式图形；Simon老大说得对，不过俺只是为了说明R并不是像Oleg老大讲的那样不能实现缩放和平移。
前几天R-help邮件列表中有人问如何计算Alpha Shape（注：Alpha Shape可以看作是闭包Convex Hull的扩展，它可以通过调整Alpha参数计算更精细的“闭包”从而大致描述平面或空间上一群点的外形），我给了一个K-Means聚类的做法：
set.seed(1234)
devAskNewPage(ask = TRUE)
par(pch = 20)
dat = iris[, 1:2]
n = nrow(dat)
for (k in 2:30) {
ch = integer()
cl = kmeans(dat, k, 50)$cluster
plot(dat, main = paste("k =", k))
for (i in unique(cl)) {
idx = chull(tmp <- dat[cl == i, ])
ch = c(ch, as.integer(rownames(tmp[idx, ])))
polygon(tmp[idx, ], border = NA, col = rgb(0, 0, 0, 0.2))
}
plot(dat, main = paste("Polygon shape when k =", k))
polygon(dat[ch, ], col = rgb(0, 0, 0, 0.2)) # need to be ordered
}
通过选择不同的k，可以逐步描述出点的外形，不过这种方法太粗略，而且最后找出来的多边形的顶点也没有排序，因此不是太好的解决方案。Alpha Shape的算法之一是用某个固定半径的圆去“套”一对一对的点，当一对点都刚好落在圆上而且圆内不包含任何其它点的时候，这两个点就是形状的边界点。通过这样的方法找出所有的边界点，便描述出了Alpha Shape。有兴趣的看官可以把这几句话转化为R代码试试。注意圆的半径是可变的参数，不同的半径对形状的描述精确程度有不同，显然当半径很大时，算法找出的就是闭包。
写出代码的看官请不吝分享一下：）
一眨眼又到12点了，ft。简写一下：Liu M.发来邮件，写了一个中心极限定理的代码，我一看，不对劲（此处省去200字）。不过另一个现象倒是让我颇为吃惊：常数的核密度估计依然是钟形曲线！刚开始我还觉得纳闷，后来转念一想，哦，也对。
不信客官您试试plot(density(rep(0, 1000)))看密度估计图形是啥样的。英文版参见这里：When “Bell-shaped” is Far Far Away from Gaussian。
本例要说明的是小学里的道理，不能学小蝌蚪找妈妈，逮着钟形曲线就说正态。
对搞统计的人来说，务必谨慎对待那些经过处理的数据，若有可能，尽量画出原始数据，而不要仅仅给出综合信息。
最近评论