统计计算
R语言统计计算
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值越大嘛），但大多少呢？其实也没大多少。
方差齐与不齐时t检验的结果对照
很久以前，某同学问我，为什么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()
没想到魔鬼经济学（Freakonomics）在NY Times的Blog上也写起R来了。话说昨天有这么一篇“免费的超牛牛计算软件”（Free Super-Crunching Software），作者小小赞扬了一把Excel可以多么灵活多么实用，然后介绍了一下R。后面顿时狂风大作、飞沙走石、昏天黑地。
其中，第一位留言的同学给了一篇小文档，对统计分析的人来说值得警醒。
很多人对统计计算或数值计算不以为然，表现在他们充分相信计算机、在遇到计算问题时就是一句话，“交给软件去做就可以了”。不过我觉得稍微了解一下数字的常识还是很有必要的，一些表达式在数学上成立，在计算上未必成立。据R-FAQ介绍，在”The Elements of Programming Style”一书中有这么一句话：
10.0 times 0.1 is hardly ever 1.0. （10.0乘以0.1通常不是1.0。）
类似地，在R里面还有这样的例子：
> .17 - .6 + .43 [1] 5.551115e-17 > .17 + .43 - .6 [1] 0 > .6 - .17 - .43 [1] -5.551115e-17
奇怪么？不奇怪（而且注意这不是R的问题），想象与现实总是有差距。以前我见有人用R求极大似然估计，直接使用了密度函数值相乘再最大化，却不知成百上千个小数乘起来在计算机中会是怎样难以表达，我换作取对数求和的形式之后，估计就稳定多了。编程的人（尤其是统计编程的人）千万莫偷懒，写代码之前一定要想好把数学问题转换一下，让计算机能够适应你要用的数字。还有一个最简单的例子就是求组合数，理论上可以用排列数（函数factorial()）乘乘除除得出来，但是当n取好几百的大数时，想想n的排列是多大的数、计算机能否表达？这种情况下就应该用组合数本身的函数choose()——理论上等价，计算上不等价。万一choose()本身也很难计算了，那只能考虑先求对数再求和最后求幂，这样也许损失精度，但计算上可实现。
最近R-help里面出现过好几例类似的案子，甚至有人在优化时初始值取的是Inf（无穷大），让人哭笑不得、无可奈何。昨天又有人问1 - cumsum(rep(0.1, 10))的最后一个数字为什么不是0。
最近评论