Keep on Fighting!

之前在群众大学给一些客官承诺过我会把我的硕士毕业论文发上来，其实早就在网上了，只不过没吱声而已，用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))

本文论点在于计算机自动化的可能性以及对长久以来的数学推导传统的疑惑，并非为了贬低经济学论文。

陈丽云这篇博客“真的是只大狐狸吗？对江西财经陈军昌博士的探究”让我想起一个长期以来我关心的一个话题。我对陈军昌这个人本身不如我对他的摘要的兴趣大，此君提到：

本文预言：在不久的未来，计算机技术将会借助非线性问题的进展彻底占据经济学的主流地位。这项技术不再是简单的用于经验数据的回归预测，而将成为主流形式化逻辑。 本文作者甚至计划在将来使用纯计算机程序的形式化逻辑写作一篇经济学论文。

我完全相信这段预言。倒不是说我觉得这种做法是对的，只是在目前可见的范围内，我严重怀疑堆积如山的经济学论文是否还需要人的脑子，我甚至想象，给一个Sweave模板，提供几个参数（如欲选择用什么模型、生成什么样式的图形），然后把数据读进来用R跑一遍，一篇论文就自动生成了，加上LaTeX本身就显得正式，这种论文一定人模狗样的，很能忽悠人。比如：

\Sexpr{names(dat)[1]} 的均值为 \Sexpr{mean(dat[ ,1])}，标准差为 xxx。图 \ref{...} 为 \Sexpr{figname[1]}：

……

\Sexpr{modname[1]}模型显示，斜率项为\Sexpr{coef(dat.lm)[2]}（t 检验结果为\Sexpr{ifelse{coef(summary(dat.lm))[2,4]<0.05, '显著', '不显著'}}）。

把目前主流杂志上的经济学论文遍历一下，总结一个八股规则，生成统计分析部分。至于结尾嘛，就把所有论文的结论部分分条存在一个数据库中，随机抽取5条就可以了，反正大家的结论都很NB，都对社会主义建设有重大意义。

调侃归调侃。如果这位陈博士的论文真的能被广为接受的话，我估摸着将来大多数期刊是不是要去喝西北风了。以上是经济学界的事情，与我没啥关系，暂不多说。还是回头说统计。

传统统计教科书大多会提及t检验中方差齐性这个问题，因为检验的假设条件是需要总体方差相等的。然而这个问题实际上可能并没有人们想象的那么重要，这里给两个简单的数值计算结果，看看方差不等对检验结果有什么影响。

par(mar = c(4, 4, 0.5, 0.5), mfrow = c(1, 2))
set.seed(123)
plot(pval <- t(replicate(1000, {
    x1 = rnorm(100, mean = 0, sd = runif(1, 0.5, 1))
    x2 = rnorm(100, mean = 1, sd = runif(1, 2, 5))
    c(t.test(x1, x2, var.equal = TRUE)$p.value, t.test(x1, x2,
        var.equal = FALSE)$p.value)
})), xlab = "P-value: equal variance", ylab = "P-value: unequal variance",
    pch = 20, asp = 1)
abline(0, 1)
plot(pval[, 1], pval[, 2] - pval[, 1], xlab = "P-value: equal variance",
    ylab = "Diff of p-values (unequal var - equal var)", pch = 20)

过程是：从两个正态总体中生成样本，第一个总体均值为0，标准差随机取自U(0.5, 1)，第二个总体均值为1，标准差取自U(2, 5)，显然两个总体标准差不相等，那么在t检验时设定和不设定方差相等的选项对结果有多大影响？把两种情况的P值都画出来：左图是原始P值，可见基本在对角线上，说明大致相等，若眼神儿不好，可看右图，即P值的差异，可见方差不等时P值偏大（原因很简单，因为Welch校正的自由度小于等于不校正的自由度，样本量相等的时候统计量的分母即标准误一样，因此统计量完全一样，自由度越小，P值越大嘛），但大多少呢？其实也没大多少。

很久以前，某同学问我，为什么R语言要叫R呢？我回答：因为它总是让人发出一声“啊！”。譬如，“啊，R语言真TM难学”，或“啊，R语言真神奇”，等等。啊，我今天在论坛上看见一个生态模拟的例子，便为发帖者编撰了如下场景：

ecol.simu = function(nr = 10, nc = 10, col.sp = c(1, 2),
    pch.sp = c(1, 2), col.die = 1, pch.die = 4, cex = 3,
    nmax = 50, interval = 1) {
    x = rep(1:nc, nr)
    y = rep(1:nr, each = nc)
    par(ann = FALSE)
    p = sample(rep(1:2, nr * nc/2), nr * nc)
    for (i in 1:nmax) {
        plot(1:nc, 1:nr, type = "n", xlim = c(0.5, nc + 0.5),
            ylim = c(0.5, nr + 0.5))
        abline(h = 1:nr, v = 1:nc, col = "lightgray", lty = 3)
        points(x, y, col = col.sp[p], pch = pch.sp[p], cex = cex)
        Sys.sleep(interval)
        idx = sample(nr * nc, 1)
        points(x[idx], y[idx], pch = pch.die, col = col.die,
            cex = cex, lwd = 3)
        tbl = as.vector(table(p))
        tbl = tbl + sign(p[idx] - 1.5) * c(1, -1)
        p[idx] = sample(1:2, 1, prob = tbl)
        Sys.sleep(interval)
    }
    p
}
par(mar = c(3, 3, 1, 1))
# 一步一步来
ecol.simu()
# 对于急性子和眼力好以及计算机速度快的同志，试试这个暴力模拟
ecol.simu(col.sp = c(8, 2), pch.sp = c(20, 17), nmax = 1000,
    interval = 0.05)

大意是：有两种物种，每次随机死一个，后继者的概率为生者的比例。

理论上应该是独立的啊，但今天偶然发现这么个现象：

不管怎么重复这幅散点图，两个累加变量之间总是有某种说不清的关系。怪哉，怪哉。

n = 1e+05
x = cumsum(rnorm(n, 0, 1))
y = cumsum(runif(n, -1, 1))
cl = gray(seq(0, 1, length = n))
plot(x, y, pch = ".", col = cl, xlab = 'CUM-Norm', ylab = 'CUM-Uniform')

数学公式推得快的看官请帮忙解释一下，谢谢！

---

下午2:00的补充：

我想明白这个道理了，因为累加之后数据的极差比较大，而每次点的位置与上一次点的位置之间的距离相对于这个极差来说就非常小，因此点的移动就是一小步一小步的，而整幅散点图看起来就会有聚类的现象，似乎两个累加变量之间有什么依随关系，实际上跟布朗运动差不多。

今天突然想起，开平方sqrt(x)和直接用x^0.5谁的速度快，本以为sqrt()会快一些，结果貌似不是这样：

> rm(list = ls(all = TRUE))
> set.seed(123)
> x = runif(1e+05)
> system.time(replicate(1000, {
+     sqrt(x)
+     NULL
+ }))
   user  system elapsed
  18.07    0.00   18.14
> system.time(replicate(1000, {
+     x^0.5
+     NULL
+ }))
   user  system elapsed
  10.47    0.00   10.50
> identical(sqrt(x), x^0.5)
 [1] TRUE

看来四则运算还是快于调用函数啊。（旁白：真煞风景，今天是圣诞前夜呀老大！我：啥时候见过美国人过五四青年节么？）

请各位客官谨记：R是一把锋利的刀，用得不好会割到自己（so sharp that you’ll cut yourself）。近日一位小同学的一个问题真是让我颇有些生气。

带着图形界面统计软件的思维来用R的话，十有八九会割到自己。在SPSS、SAS等工具中，虽然看着满屏幕的按钮，但大部分人可能也不害怕，因为不用管它们是什么意思，瞎选一通，按OK，下面就可以洋洋自得看着长篇大论的报表出来了，这种过程很是爽。

到了R的世界，满屏幕只有代码，后来好不容易明白了，原来R不用编程，调用现成函数就可以了，于是乎，开始把各式各样的数据、参数往函数里面扔，扔完了summary()一下，长篇大论的报表也出来了，甚爽。直到有一天，R向你报告说某地方出错了，于是傻了。

这里的案例是AdaBoost，这位同学用adabag包中的adaboost.M1()函数对树模型做boosting，却被告知无法进行。我看了一下数据，原来因变量是数值变量。于是火了，数值变量你咋用Adaboost.M1啊？它本身是对分类问题做的提升，对于一个回归问题非要驴唇对马嘴，这不净瞎扯么。

洒家满以为是个有趣的问题，结果饿着肚子回了邮件，真是亏大了。外专业的同仁也就罢了，俺不会说什么，关键是统计专业的。挥一挥衣袖，用膳去鸟。

今日在北大听史密斯商学院创业大赛报告，其中有一位参赛者幻灯片中提到了他们的用户数这几年呈几何级数增长，并拿用户数和时间作了一幅图，图中线条呈现出增长越来越快的趋势，其实这种做法有糊弄之嫌——增长越来越快的并不一定是几何级数增长方式。例如y=sin(x)+1在区间上增长也是越来越快，但它并非几何级数。

表达几何级数增长（或者指数增长）的方式一般是对y取对数，然后与x作图，看图形是否呈一条直线：若log(y)=a*x+b，那么显然是指数增长方式。人眼观察直线比观察曲线要容易得多，因此这种方法比用原始数据作图要更容易表达“几何级数增长”。下图左边为原始数据，右边为y轴取对数后的图形。R中处理起来非常简单，作图时添加参数log即可（可以对x轴或y轴或者同时取对数）。

文件下载链接

#png("exp_growth.png", width = 600, height = 500)
options(scipen = 5)
x = seq(1.51 * pi, 2 * pi, length = 100)
par(mfrow = c(2, 2), pch = 20, mar = c(5, 6, 1, 0.1),
    col = rgb(0, 0, 0, 0.5), las = 1, mgp = c(4, 1, 0))
plot(x, sin(x) + 1)
plot(x, sin(x) + 1, log = "y")
plot(x, exp(x))
plot(x, exp(x), log = "y")
#dev.off()

如果我没记错的话，Prof Lin在课堂上提到，Least Median of Squares是一个难题：它的目标是最小化误差平方的中位数。当时上课我没去考虑这个问题，晚上突然想起，这个问题不管怎么说至少可以有暴力解法，比如用BFGS等方法去优化就能算出系数了。

但是转念一想，LMS这个东西怎么看着眼熟，打开MASS一看，原来是自己以前用过的。MASS包的lqs()可以计算各种稳健和耐抗回归，LMS是其中一种。LMS是一种稳健回归方法（下第一幅图），然而MASS的文献指出它对大量堆在中间位置上的数值敏感（下第二幅图）。

图1 离群点对普通最小二乘有严重影响，LMS却没问题

图2 大量的中间值让LMS失效，普通最小二乘却没事

文件下载链接

#png("lms.png", width = 500, height = 400)
par(mar = c(4, 4, 0.1, 0.1))
library(MASS)
set.seed(1)
x = runif(50)
set.seed(1)
y = 2 + 3 * x + rnorm(50)
# an outlier
x = c(x, 2); y = c(y, 30)
plot(x, y, pch = 20)
abline(lqs(y ~ x, method = "lqs"), col = "blue", lty = 2)
abline(lm(y ~ x), col = "red")
legend("topleft", legend = c("OLS", "LMS"), col = c("red",
    "blue"), lty = 1:2, bty = "n")
#dev.off()

#png("lms-central.png", width = 500, height = 400)
par(mar = c(4, 4, 0.1, 0.1))
library(MASS)
set.seed(1)
x = runif(50)
set.seed(1)
y = 2 + 3 * x + rnorm(50)
# 500 central values
x = c(x, jitter(rep(mean(x), 500), 3))
y = c(y, jitter(rep(mean(y), 500), 3))
plot(x, y, pch = 20, cex = c(rep(1, 50), rep(0.1, 500)))
abline(lqs(y ~ x, method = "lqs"), col = "blue", lty = 2)
abline(lm(y ~ x), col = "red")
legend("topleft", legend = c("OLS", "LMS"), col = c("red",
    "blue"), lty = 1:2, bty = "n")
#dev.off()

应读者要求，将两幅图的数据分别附在这里。

第一幅图的数据下载

第二幅图的数据下载