Keep on Fighting!

暑假回国，办完一堆事情，回家躲在小村子里悠哉悠哉，无丝竹之乱耳，无案牍之劳形，谈笑木有鸿儒，往来全是白丁^注1。那古人云“春困秋乏夏打盹，冬日漫漫睡不醒”。本小子虽然贪睡，但在闭关写书这件事情上不含糊。其实与其说在写，不如说在复制，因为以前已经用纯LaTeX写了很多内容，后来看到pgfSweave包和LyX之后，觉得纯LaTeX太麻烦了，于是把原来的内容都复制到LyX中接着写。到从村子里出关时，R的基础图形（graphics包）已经全部写完，包括常见的底层作图函数和几乎所有高层作图函数。

后来回到有网络的世界，git push一下，就同步上网了，供各位客官公开下载（见于本站“作品”页面）。过了一阵子，领导的领导们看见了，觉得这样做有点不靠谱。再过了一阵子，又有好几个人建议我不要放出来公开下载。我呢，也绝非来自石器时代、对“世风日下人心不古”一无所知，佛曰“我不下地狱那哪个来下地狱呢？”佛敢下地狱，那是因为他是佛，是佛又咋地？佛有核心竞争力啊。现在这个忙碌年代，我们（别人）想尽办法折腾别人（我们），折腾的途径就是评价。网站不看内容要看备案，活人不看要看KPI，百姓生活不看要看GDP……而混学界的人就得靠著作来评价，于是使出浑身解数去发表，自己努力也好，坑蒙拐骗也好。我这样赤果果把作品放在网上供任何人下载，说出来真是一件很恐怖的事情，娘亲呐，多少豺狼虎豹盯着呢。

本来想继续扯到李一、马云等人，似乎有点扯远了，我是在想像马云这样伟大的宗教领袖式人物为何需要道教思想的支撑。俺们这个时代的人，真的越来越没底气了？打住打住。

最后我接受了腾飞的意见，把原本公开的书稿撤下来了。倒不是怕盗版抄袭，主要是一本不完整的书稿，对读者似乎不是一件好事，等正式完成之后，再集中发给大家提意见，这样效用会大一些；再想一想，出版之前公开下载也会对出版社带来一些麻烦，于是撤销下载。当然，目前我仍然内部开放书稿，我熟悉放心的人以及答应我在COS论坛帮别人回答图形方面100帖的人（鼓励为别人做贡献的人），皆可获得我的书稿。这本书的目录部分可以公开浏览：

《现代统计图形》目录下载（PDF）

各位客官若有任何建议，请随时与我联系，我会在完稿前尽量多采纳读者的意见。

—————————————
注1：平民老百姓。

李丰同学建议在《现代统计图形》书稿中插入一些彩蛋，我想第一颗蛋就用这些杂乱无章的彩色点线吧。不过小的实在拿不定哪一幅图看着更“艺术”一些，因此请路人帮忙挑选一下。下图只是众多可能的随机图中的一幅，更多图形请下载查看后面的PDF文件。挑好之后把图中的起始数字告诉小的，小的会在书稿的致谢部分说明您的贡献。

下载全部40页PDF图形文件（2.7Mb）

之前在群众大学给一些客官承诺过我会把我的硕士毕业论文发上来，其实早就在网上了，只不过没吱声而已，用Git同步一下容易，但写几个字介绍一下就花时间了。感兴趣的客官可以去本站作品页面内下载：http://yihui.name/cn/publication/#GSM

这论文呢，自我评价可以给个80分。最大的问题在于没有花时间去整理文章的结构，所以构架上稍微有点散乱（俗称“意识流”）。内容上熟悉我的博客的客官一眼就能看出来，其实都是些博客文章的汇总，只不过用LaTeX让它们变得“人模狗样”一些而已，好在本小子平时也积攒了这些鸡零狗碎的东西，动过自己的脑子。我觉得群众大学的毕业论文，很多都是一个套路：经济/金融数据套一个神奇的模型，直到最后整个世界一片和谐，读者在最后一章都能隐约看到上帝老爷子在朝你挥手。其实也没啥，找工作不容易，地球人也都知道写论文就是忽悠——漫漫人生路上一道工序。

由于本小子是个小人（小小的活人），所以总关心小人关心的事情（俗称“人本主义”）。这论文嘛，窃以为也没什么上下高低之分，说出你怎么想的就可以了，而不要总说“他们”怎么想怎么做。一定要有数学上的创新？一定要有人家看不懂的公式才是好论文？一定要有综述？一定要有长长的参考文献列表才是好论文？一定要板起脸？不能写八卦？不准幽默？……嗨，作茧自缚。几年前看到一篇好文章，颇具恶搞性质，建议各位客官收藏：How to write Consistently Boring Scientific Literature。

言归正传：本文是厌倦八股文和数学理论的产物，从理论角度来说，几乎没什么价值，不过这篇文章是用Sweave写的，完全具有可重复性和100%透明度，对文中结果有怀疑的客官可以自行运行代码；其次，统计模拟和图形的声音在界内太微弱，大家都很忙，有人在忙着推公式，有人在忙着编数据，有人在忙着把公式用到不知道是不是编出来的数据上，本小子跟着瞎掺和了点别的东西，仅此而已。甭管有用没用，敬请拍砖。

－－－－－－－－－－外一篇：坛霸是怎样练成的－－－－－－－－－－

曾经有童鞋称呼在下为“坛霸”，这个……有时候确实有那么点意思，无图无真相（两个多月没怎么回帖了，一鼓作气）：

接下来我会陆续写第三届中国R语言会议、《现代统计图形》书稿和useR! 2010，若时间允许，我考虑一下电视剧《九阴真经》（93版）。

如同经济学家不讲道德一样（学过经济学的人都知道这句话的意思），理论统计学家从某种程度上来说也不讲道德。我们常用的一些统计量通常都渐近服从某种分布（以卡方和正态为典型），看起来做理论的人对这些渐近理论都非常骄傲和自豪，我们在学习过程中也要一代一代传承下去。数学公式摆出来当然能唬人，也许唬到最后大家都为光着屁股的皇帝欢呼。坦白说，我对这些东西感到非常厌倦。

近日来收到邮件少了，但各个问题都不太好直接回答。比如这则关于McNemar检验的问题：McNemar检验可以有两种形式的统计量，一为(b – c)²/(b + c)，一为2b*log(2b/(b+c)) + 2c*log(2c/(b+c))，其中b和c是列联表非对角线上的频数。前者是McNemar检验本身的统计量，可以根据渐近正态分布得来（然后平方得到卡方），后者是似然比统计量（不带约束的似然除以带约束的，取对数，乘2）。McNemar检验看似复杂，实际上可以简化为检验b = c，或等价于检验一个n = b+c的二项分布中，是否p = 1/2（观察到X = b或c）。现在的问题是，这两种统计量有没有优劣之分？

作为一个懒得推公式的人，我向来喜欢用模拟回答问题，因为模拟的结果非常直截了当。我的考虑是，要看渐近统计量的优劣，那就看随着n增大，统计量和渐近分布有多接近好了。一个自然而然的想法当然是对若干统计量的观测值做分布检验了，比如KS检验。我们知道这两个统计量都是自由度为1的卡方分布，剩下的事情就是计算：

set.seed(123)
nmax = 1000
p = matrix(nrow = nmax, ncol = 2)
for (n in 2:nmax) {
    # 生成服从二项分布的随机数，分别计算两种统计量并作KS检验、记录P值
    b = rbinom(500, n, 0.5)
    x1 = (b - (n - b))^2/n
    x2 = 2 * b * log(2 * b/n) + 2 * (n - b) * log(2 * (n - b)/n)
    p[n, 1] = ks.test(x1, "pchisq", df = 1)$p.value
    p[n, 2] = ks.test(x2, "pchisq", df = 1)$p.value
}
# 调整一下数据格式，画图：随着n增大，P值如何变化？
library(ggplot2)
d = melt(p, varnames = c("n", "method"))
d$method = factor(d$method, labels = c("McNemar", "LRT"))
colnames(d)[3] = "p.value"
qplot(n, p.value, data = d, shape = method, geom = c("smooth", "point")) +
    scale_shape_manual(values = c(2, 3))

昨天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会议所致。剩下的联系人就分布在五湖四海了，结果请看：

PDF图形文件下载链接

频数文件：

CSV频数文件下载链接

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的基础图形元素什么图都可以画，从而不愿学习新的作图系统。