陈丽云这篇博客“真的是只大狐狸吗？对江西财经陈军昌博士的探究”让我想起一个长期以来我关心的一个话题。我对陈军昌这个人本身不如我对他的摘要的兴趣大，此君提到：
本文预言：在不久的未来，计算机技术将会借助非线性问题的进展彻底占据经济学的主流地位。这项技术不再是简单的用于经验数据的回归预测，而将成为主流形式化逻辑。 本文作者甚至计划在将来使用纯计算机程序的形式化逻辑写作一篇经济学论文。
我完全相信这段预言。倒不是说我觉得这种做法是对的，只是在目前可见的范围内，我严重怀疑堆积如山的经济学论文是否还需要人的脑子，我甚至想象，给一个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，都对社会主义建设有重大意义。
调侃归调侃。如果这位陈博士的论文真的能被广为接受的话，我估摸着将来大多数期刊是不是要去喝西北风了。以上是经济学界的事情，与我没啥关系，暂不多说。还是回头说统计。
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检验的结果对照
今天突然想起，开平方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
看来四则运算还是快于调用函数啊。（旁白：真煞风景，今天是圣诞前夜呀老大！我：啥时候见过美国人过五四青年节么？）
没想到魔鬼经济学（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。
近期评论