R 语言
强大的统计计算和作图语言
本来这篇流水账日志早就该发了,一直懒得写,放在草稿箱里已有一个月,还是拖出来简单写一写算鸟。正好开学前老板说,你娃儿去开了useR!会议,那开学第一周的组会上就你和Marie来说说你们的观察吧。于是乎,我又绞尽脑汁回忆了useR! 2010。
公元2010年7月19日凌晨3:45起床,搭租好的车从农村去机场,一小时后到达。路上看见天边电光闪闪,非常吓人,我心想完鸟,这飞机还能飞么?机场一路进去,木有任何延迟或取消航班的迹象,早上6点,按时起飞,我心里捏着一把汗。这机长也不是盖的,拖着飞机掉头就往没有闪电的方向飞,绕道去芝加哥,幸好菩萨及圣母保佑,平安到芝加哥ORD机场,大约40分钟后,搭下一班飞机去DC里根机场,算是故地重游。出机场又找预定好的车,开了一个小时才开到马里兰NIST,路上又经过那华盛顿纪念碑啊、林肯纪念堂啊、造币厂啊之类的地方,也不稀罕。
这次张翔张主席也从国内不远万里漂洋过海到NIST开会,我琢磨着他可能是唯一一个从中国大陆来的参会者,报名的人当中还有另一位兄弟似乎也要来,但后来不知为何掉链子了。张主席本来订好了租车,结果到了美帝,驾照上木有对应的英文翻译,于是也木有办法租车,而本小子至今还没学开车,那么就地铁公交加“11路”车,总之找到了旅馆。张主席到了NIST附近,没直接去旅馆,而是气定神闲找了家餐馆吃饭去了,我只好吭哧吭哧去找他吃饭。主席就是主席,也不是盖的,不拖箱子不拎包,而是整个登山包背着,连睡袋装备都有,那是相~昂~当的专业。
从农村去了城市,自然是不便,虽说是城里,NIST周围似乎也是鸟不拉屎地不长毛,找了家超市胡乱买些食物水果,回旅馆歇着。下面进入正题。
后来回到有网络的世界,git push一下,就同步上网了,供各位客官公开下载(见于本站“作品”页面)。过了一阵子,领导的领导们看见了,觉得这样做有点不靠谱。再过了一阵子,又有好几个人建议我不要放出来公开下载。我呢,也绝非来自石器时代、对“世风日下人心不古”一无所知,佛曰“我不下地狱那哪个来下地狱呢?”佛敢下地狱,那是因为他是佛,是佛又咋地?佛有核心竞争力啊。现在这个忙碌年代,我们(别人)想尽办法折腾别人(我们),折腾的途径就是评价。网站不看内容要看备案,活人不看要看KPI,百姓生活不看要看GDP……而混学界的人就得靠著作来评价,于是使出浑身解数去发表,自己努力也好,坑蒙拐骗也好。我这样赤果果把作品放在网上供任何人下载,说出来真是一件很恐怖的事情,娘亲呐,多少豺狼虎豹盯着呢。
本来想继续扯到李一、马云等人,似乎有点扯远了,我是在想像马云这样伟大的宗教领袖式人物为何需要道教思想的支撑。俺们这个时代的人,真的越来越没底气了?打住打住。
最后我接受了腾飞的意见,把原本公开的书稿撤下来了。倒不是怕盗版抄袭,主要是一本不完整的书稿,对读者似乎不是一件好事,等正式完成之后,再集中发给大家提意见,这样效用会大一些;再想一想,出版之前公开下载也会对出版社带来一些麻烦,于是撤销下载。当然,目前我仍然内部开放书稿,我熟悉放心的人以及答应我在COS论坛帮别人回答图形方面100帖的人(鼓励为别人做贡献的人),皆可获得我的书稿。这本书的目录部分可以公开浏览:
各位客官若有任何建议,请随时与我联系,我会在完稿前尽量多采纳读者的意见。
—————————————
注1:平民老百姓。
我第一天讲的统计图形讲稿后来已经更新传上来了,由于麦克风的不平稳性,估计很多人都没听清我在说什么,不过这也没啥,因为我讲的内容都已经出现在我的《现代统计图形》书稿中了,想要详细了解的读者可以去作品页面中下载不完整的书稿。第二天的演讲幻灯片和代码也都发布在统计之都上的第三届中国R语言会议纪要中了。
尽管中国R语言会议还未形成燎原之势,但八卦的R core们其实还是在偷偷关注的。我这次去马里兰NIST参加useR! 2010找机会又跟Martin Machaeler提了提我们的“山寨会议”,这位老大一边摆弄着手里的Emacs一边听我忽悠,到最后走出会场我才发现他其实没有真正想起来我是谁,因为他最后猛然想起来R Journal上貌似有一篇中国R会议的报道……老大们表示他们个人还是很有兴趣去一趟中国的,不过要把这三四百观众都带去恐怕有点困难,路费啊签证啊都麻烦。所以呢,我们下次可以考虑请一些有影响力的R core们去中国,办一届英文会议也未尝不可。
回头来看,这次会议的主要问题还是准备仓促(难道是我留下的病根?),收到的演讲不够覆盖R的功能;其次就是交流时间不够,没有有意制造足够的交流机会,当然这与会议时间长度有关,但会议日程安排上应该可以改进,用制度促进交流;最后,会议的产出也不明确,有个出版相关的小组就好了,即使不出版,发表在COS主站上也好,开完会就散伙,可惜了儿了。回头再线下商量。
会议T恤是这次会议的一个小亮点,下次要保持。顺便附图一张:某童鞋把这次R语言会议的T恤穿到了美帝农村钓鱼(其实一条鱼都没钓到)。
下期预告:《现代统计图形》书稿。这两天忙着搬家,估计要等几天了。
这论文呢,自我评价可以给个80分。最大的问题在于没有花时间去整理文章的结构,所以构架上稍微有点散乱(俗称“意识流”)。内容上熟悉我的博客的客官一眼就能看出来,其实都是些博客文章的汇总,只不过用LaTeX让它们变得“人模狗样”一些而已,好在本小子平时也积攒了这些鸡零狗碎的东西,动过自己的脑子。我觉得群众大学的毕业论文,很多都是一个套路:经济/金融数据套一个神奇的模型,直到最后整个世界一片和谐,读者在最后一章都能隐约看到上帝老爷子在朝你挥手。其实也没啥,找工作不容易,地球人也都知道写论文就是忽悠——漫漫人生路上一道工序。
由于本小子是个小人(小小的活人),所以总关心小人关心的事情(俗称“人本主义”)。这论文嘛,窃以为也没什么上下高低之分,说出你怎么想的就可以了,而不要总说“他们”怎么想怎么做。一定要有数学上的创新?一定要有人家看不懂的公式才是好论文?一定要有综述?一定要有长长的参考文献列表才是好论文?一定要板起脸?不能写八卦?不准幽默?……嗨,作茧自缚。几年前看到一篇好文章,颇具恶搞性质,建议各位客官收藏:How to write Consistently Boring Scientific Literature。
言归正传:本文是厌倦八股文和数学理论的产物,从理论角度来说,几乎没什么价值,不过这篇文章是用Sweave写的,完全具有可重复性和100%透明度,对文中结果有怀疑的客官可以自行运行代码;其次,统计模拟和图形的声音在界内太微弱,大家都很忙,有人在忙着推公式,有人在忙着编数据,有人在忙着把公式用到不知道是不是编出来的数据上,本小子跟着瞎掺和了点别的东西,仅此而已。甭管有用没用,敬请拍砖。
----------外一篇:坛霸是怎样练成的----------
曾经有童鞋称呼在下为“坛霸”,这个……有时候确实有那么点意思,无图无真相(两个多月没怎么回帖了,一鼓作气):
接下来我会陆续写第三届中国R语言会议、《现代统计图形》书稿和useR! 2010,若时间允许,我考虑一下电视剧《九阴真经》(93版)。
先说培训。若一切顺利,我们将在6月14日进行第三届中国R语言会议的会前培训,暂定由刘思喆和我来讲。本次培训计划上下午各3小时,培训费用=R的费用(即:free as in beer)。我的计划如下:
下载:幻灯片《现代统计图形》 培训时间:2010年6月14日下午2:00至5:00培训目标:1、了解图形的构成元素,初级用户可以知道作图的诸多可能性,高级用户可以任意自定义图形;2、了解统计图形的基本类型和适用情形,跳出单调的“饼图+条形图+折线图”的范围;3、了解R的四种图形系统,即基础图形系统、grid、 lattice和ggplot2;4、学会用图形辅助模型去探索和分析数据;5、了解其它靠谱或不靠谱的应用,如统计动画和交互式动态图形
培训内容:R语言是统计计算和统计作图的强有力工具,本次培训着重介绍后者,内容包括:对统计图形历史的简要回顾,说明统计图形的功能(约5分钟);介绍R语言的基础图形系统的基本构成,包括各种图形参数和基本图形元素(约40分钟);介绍R自身的 graphics包中的各种统计图形函数,包括直方图、等高线图、散点图矩阵等(约1小时);R附加包中的各种图形函数,包括地图、脸谱图、平行坐标图等(约20分钟);R的其它三种图形系统:grid、lattice和ggplot2(约20分钟);基于统计模型的图形应用,包括回归模型、主成分分析、光滑方法、分类与回归树等(约25分钟);其它图形应用,包括动画和交互式图形等(约10分钟)。本次培训的大纲主要遵循本人正在编写的《现代统计图形》一书,该书的不完整书稿可以从这里下载:http://yihui.name/cn/publication/
适合听众:推荐以下四类听众前来听课(1)英语阅读能力较好,对编程感兴趣(2)公司企业的数据分析人员,尤其是咨询公司(3)教数据分析相关课程的高校老师(4)领导的秘书
不适合听众:本培训为实用性质的培训,可能不适合以发表学术论文为目标的听众,统计图形看似过于浅显,在统计学术界非主流研究方向
近日来收到邮件少了,但各个问题都不太好直接回答。比如这则关于McNemar检验的问题:McNemar检验可以有两种形式的统计量,一为(b – c)2/(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,都对社会主义建设有重大意义。
调侃归调侃。如果这位陈博士的论文真的能被广为接受的话,我估摸着将来大多数期刊是不是要去喝西北风了。以上是经济学界的事情,与我没啥关系,暂不多说。还是回头说统计。
洒家真是牟想到啊,Dirk Eddelbuettel居然也会给我发邮件问animation问题。关注R与金融的人应该知道这个人,不过这一点可能只是他的副业,更多工作在于研究R的高性能计算以及维护Debian上的R以及其它Debian包。蓬荜生辉啊。本小子的动画包,自打一开始就没做任何宣传,后来中奖了也没在R里吱声。前些日子某波兰童鞋来信问动画问题,我忍不住回问了一下你们从哪儿获知的消息,结果该童鞋说,老师布置作业,班上一共30来人,每个人分50来个R包去研究,CRAN上所有的R包就差不多研究遍了,我晕。
前两天,中国群众大学校报小记者来信,要采访统计之都的“幕后黑手”,据说是校报办公室老师安排来采访的,我一看心里大概就明白咋回事了:八成可能跟群众大学的某大领导有关……其实COS早在办旧版的时候就已经惊动党中央,院领导对我提起过几次说上面对COS赞赏有加,你娃儿要好好办,本小子连声“嗻”。
宣传这事儿,我一直认为没太大必要,理由是自身不够牛的时候把别人引来会让别人失望(宣传提高了别人的预期),而采取小步渗透的方式则会让用户有足够的惊喜,并且以这种方式聚集的人群关系会比较紧密。
现在用Google搜“统计学网站”排名第一是统计之都,搜“统计学论坛”排名第一是COS论坛,而我们从未用任何手段(如邪恶的积分制度)引诱任何会员在任何地方宣传。说到底,踏实做事才是王道。不过既然台风来了猪都会飞,这次有了机会,飞它一把也无妨。
~~~~~~~~~~~~~~自相矛盾的分隔线~~~~~~~~~~~~~~
最近又写了个跟WinBUGS(或者OpenBUGS)有关的R包,取一很俗的名字,叫iBUGS,之所以加个i,是因为我想让它尽量智能。目前实现了以下几点智能:
- 智能找WinBUGS或OpenBUGS的安装路径,免得手工填写;
- 智能分析BUGS模型,找出参数列表,可用鼠标点选,免得手工填写;
- 分析R工作空间,提出数据变量列表,可用鼠标点选;
这个包主要是为了提供一个(聪明的)图形界面,它包括了R2WinBUGS包中bugs()函数的所有功能。可在图形界面中写BUGS模型,然后点“执行”按钮就可以跑MCMC了。向不做贝叶斯的童鞋们特别声明一下,BUGS和bug没有任何联系。
话说暑假里的R会议该正式进入议事日程了……
source()函数是可以输出干净整齐的R代码的,但从某个R版本开始,这项功能被去掉了,让我感到很不爽,于是乎,研究了一下前一个版本的R源文件,悟出了整理代码的本质,于是操刀写了一个清理代码的函数;很久以前,小邱想出了一个聪明的办法,使得这个函数可以保留某些注释语句(以前这个函数只能去掉所有注释);几天前,我在看gWidgets包的同时,为这个函数写了个界面。一个新的R包formatR就诞生了。
很多人写R代码都不管代码是否整齐可读,也许是没有受够读代码的折磨:本来读别人的代码就是很痛苦的事情,如果代码写得不整齐,就更要命了。幸好R可以自己处理自己的代码,以R攻R,就可以得到自动整理的整齐代码了。
这个包的界面很简略,就是一个文本框加上若干按钮。以下是效果图:
感兴趣的客官可以从svn获得最新源代码自行编译安装,注意这个包依赖于gWidgetsRGtk2包:
svn checkout svn://svn.r-forge.r-project.org/svnroot/animation/pkg/formatR R CMD INSTALL formatR
初级用户可以从CRAN安装(CRAN上的版本取决于我的提交,可能会有延迟):
## 启动R,然后
install.packages('formatR')
library(formatR)
## 将自动启动界面
第二个版本中加入了执行代码、字体设置以及其它选项。需要提醒的有两点:
- 这个代码清理功能只能保留整行的注释(也就是单独占一行的注释),行内注释会被删掉。
- 目前暂时还没搞清楚中文编码的问题,所以暂不支持中文,如果代码中有中文,将会被转换为乱码。有谁熟悉gWidgets的话可以帮帮我。






近期评论