用过底层语言做计算的人转入R语言的时候一般都会觉得R语言的运算太慢，这是一个常见的对R的误解或者对R的设计的不理解。在二三十年前Chambers等一群人在贝尔实验室设计S语言之前，统计计算主要是通过那些底层语言实现的，典型的如Fortran。当时有一个基于Fortran的统计分析库，用它的麻烦就在于无论做什么样的数据分析，都涉及到一大摞繁琐的底层代码，这让数据分析变得很没劲，统计学家有统计学家的事情，天天陷在计算机程序代码中也不是个办法，要摆脱底层语言，那就只能设计高层语言了。有所得必有所失，众所周知，高层语言一般来说比底层语言低效，但对用户来说更友好。举个简单的例子，用C语言计算均值时，我们得对一个向量（数组）从头循环到尾把每个值累加起来，得到总和，然后除以向量的长度，而均值在统计当中简直是再家常便饭不过了，要是每次都要来这么个循环，大家也都甭干活儿了，天天写循环好了。 前两天COS论坛上有个帖子提到“R语言的执行效率问题”，大意如下： 有120000行数据，用perl写成12万行R命令做t.test，然后执行，大概几分钟就算完了。如果只用R语言，把所有数据先读入，然后用循环，每一行做t.test，花了几个小时，到现在还没有算完。这说明一个问题，在R里执行单行命令要比用循环快，涉及到循环的问题，最好写成单行命令执行。 我不清楚作者是如何写这“12万行”R命令的，本文假设是把t.test(dat), i = 1, 2, ..., 120000写入一个文件，然后用source()执行之。面对这样一个问题，我们有若干种改进计算的可能性。首先我们看“硬”写入程序代码的方法： ## 为了使计算可重复，设定随机数种子 set.seed(123) ## 生成数据，随机数，120000行 x 100列矩阵 dat = matrix(rnorm(120000 * 100), 120000) nr = nrow(dat) nc = ncol(dat) ## 六种方法的P值向量 p1 = p2 = p3 = p4 = p5 = p6 = numeric(nr) ## via source() f = file("test.t") writeLines(sprintf("p1 = t.test(dat)$p.value", seq(nr), seq(nr)), f) system.time({ source(f) }) # user...
阅读全文 »

本文是对R中精算学专用包actuar使用的一个简单教程。actuar项目开始于2005年，在2006年2月首次提供公开下载，其目的就是将一些常用的精算函数引入R系统。目前，提供的函数主要涉及风险理论，损失分布和信度理论。 如题所示，本文是我在学习actuar包过程中的学习笔记，主要涉及这个包中一些函数的使用方法和细节，对一些方法的结论也有稍许探讨，因此能简略的地方简略，而讨论的地方可能讲的会比较详细。闲话少提，下面正式开始学习！
阅读全文 »

分层线性模型 (Hierarchical linear Model，简称 HLM，又称多层线性模型，Multilevel Linear Model)，HLM6.0 是分层线性模型软件，包含线性和非线性部分，可以读取大部份统计软件的数据如 SPSS, SAS, SYSTAT及STATA等等。HLM常用于社会科学和行为科学，因为它常有嵌套结构(Nested Structure)的数据，因此需用次模型(Sub-Model)或分层线性模型(Hierarchical Model)，HLM就是设计来专门解决此类问题的，HLM提供的模型包括2-level models、3-level models、Hierarchical Generalized Linear Models (HGLM)和Hierarchical Multivariate Linear Models (HMLM)等。 下面是中国人民大学统计学院分层线性模型的课件，主要内容为HLM6.0的操作 下载地址：http://cos.name/wp-content/uploads/2009/09/intro2HLM6.zip
阅读全文 »

这个青年的经历，只代表他个人，没有任何群体的意义。 我想写下一段自白,这自白既是我个人的,也具有普遍意义,因为一个人经历过的事情所有的人都可以经历。 /*跟武汉博文视点合作，召集些身边的朋友，2009应届生，计算机背景，在毕业之前，讲讲自己求学、实习、找工作等的经历与感悟，文章将由电子工业出版社结集出版，在今天秋季学期开学之前出来。我是主编，也是作者之一，刚好经历跟大伙有重叠：经济学、软件工程、SAS、统计学、数据挖掘，文章贴出来，大伙多拍砖头。文中我提到COS，COS让我长见识，又结识不少好朋友。*/ 想启蒙至今，我的学习重心从文史哲转到经济学，最后以软件工程收场；在行业方面，在一家软件公司实习了一年多，找工作以药厂收场——用一个朋友的话讲，“背景是复杂得一塌糊涂”。虽然我常以“学一行，爱一行；做一行，爱一行”自勉，有时也不免回首，想想以前经历过的所有分岔口，玩味些那些没有走过的路，无论当时我多么坚决或者犹豫，一条路被选中，然后一路曲折至今。很多事情，很多选择，可能只有以后才能适合评说，现在我选择把它们记下来，在这么一个时点，我就要毕业，我有一个安身立命的去处，我有规划，但我不知道以后生活会给我什么样的分岔口，就跟我以前遭遇的一样。 ——————我的人生轨迹————— 项羽，“学书不成，去；学剑，又不成”。 —-《史记·项羽本纪 1.高中：文科生，喜文史 | |高考 ...
阅读全文 »

介结了在R中如何实现极大似然估计的一种实现方法。当然其在R中还有其它的实现方法。
阅读全文 »

R主要面向统计计算，似乎很少会用到面向对象的编程方法。但在统计计算中，在下列情形中使用面向对象的编程方法可以编程更有效率。 1）当需要用一种新的方式来表示数据，该方式与已有的数据类型有区别的时候。 2）当需要一个新的函数，该函数可以根据不同的参数类型做出不同的反应的时候。 在R中，经常需要定义一个新的函数，并且定义一个新的函数也是一项繁重的工作。相反，较少去定义一个新的类。但有时候定义一个类是一个很关键的步骤。一个类通常决定了如何对对象进行处理，决定了对象中应当包含什么样的信息。甚至有时候，类的定义决定你的项目的成败。(本文的PDF版本:R中的面象对象编程学习笔记)
阅读全文 »

鉴于最近有不少人在讨论用R软件绘制地图的问题，我也就跟着凑了凑热闹，对相应的方法学习了一番。下面的这篇文章是一个初步的介绍，还有很多内容仍在学习和探索中，如果大家有什么意见或建议，我将根据自己学习的情况对文章进行进一步的补充。
阅读全文 »

如何生成一个GeoBUGS格式的中国地图 第一节 导言 之前有些对GeoBUGS感兴趣的同学发邮件询问我有没有GeoBUGS的中国地图，以用于分析中国国内的一些空间数据。我想有必要将如何生成GeoBUGS格式的地图的方法分享给大家。这样的话， GeoBUGS就可以真正为我们所用，从而对于其他GeoBUGS没有自带的地图，我们也可以轻松生成了。本节不涉及统计分析，仅为GeoBUGS的研究使用者提供一个软件使用的技术参考。关于GeoBUGS的统计的书，国外实在是很多了，但关于这块地图定制的参考资料较少，故提出来供大家参阅。
阅读全文 »

新华网刚发布了一个关于学生冬季长跑的调查结果（于2009年4月27日13:52访问），一共调查了100人，结果中却出现了92.79%这样的比例数字，有常识的读者都知道，世上不存在0.79个人，因此这里面必然有某个地方是错的（姑且不妄言造假）。这则消息让我马上想起《统计陷阱》这本书，我们生活中有多少陷阱呢？ 从消息编辑人员的角度来说，他们可能觉得保留两位小数显得“精确”，而这种“精确精神”从数学的角度来说显得既可爱又可笑，如果小数位能表示精确，那何不保留100位小数呢？ 从统计人的角度来看，对这种调查报告中的比例数字应该有足够的警觉。很多调查报告并不会告诉我们究竟样本量多大（在这一点上新华网的调查网还比较诚实），这种情况下，我们应该弄清究竟调查了多少对象，当样本量很小的时候我们会怀疑这个调查的代表性。当我们看到比例66%的时候也许能想起来这是2/3（猜测样本量是3的倍数），但对29.1667%这个比例我们未必能很快反应出来分子和分母是多少，若报告公布方没有说明样本量，我们只能自己猜测；对于667这样的数字，我们很容易猜测这是6循环的四舍五入。最终大概思路就是拿比例去挨个乘以一系列整数，看看哪个结果接近整数，从而“还原”原来的分式n/N。以下是简单的R代码： > digit = ((1:100) * 0.29166666)%%1 # 整除1之后的“余数” > plot(digit, ylim = c(0, 1)) > idx = which((1 - digit) < 1e-05 | (digit - 0) < 1e-05) # 与0或1很靠近时 > points(idx, digit, pch = 20) > abline(v = idx, lty = 2) > axis(3, idx) > idx * 0.29166666 7 14 21 28 我们很容易发现分母（样本量）是24的倍数，因为上图中24的倍数乘以29.167%得到的结果很靠近整数；而具备某种特征的样本数量为7的倍数。根据具体的调查背景，我们可以自己猜测报告方究竟调查了多少人：24人？太少；960人？为什么不是1000人？…… 上面只是统计侦查的小游戏而已，当我们具备更多统计知识储备之后，便可以去考虑一些具体的统计模型输出是否存在造假嫌疑。我想，P值在0.05之下且很靠近0.05的时候，或P值一律接近于0的时候，我们不妨以小人之心揣测这个模型也许有问题。当因子分析中，50个变量能根据载荷被准确划分到作者预先设定的5个因子中时（5列因子载荷一律都是只在某个因子上取值极大），这个分析也许存在嫌疑。当然，所有的“小人之心”的前提假设都是：理想情况在现实中是不容易出现的（这是赤裸裸的假设检验逻辑）。
阅读全文 »

今天我来谈一点用R编程的经验吧。好像R的很多方面许多牛人都谈过了，比如R的打包啊，R的图形啊，下面我来谈谈R的测试。如果希望真正学到什么的话，还是要自己花时间实践的。 在很多情况下，你自己或和你的团队在一起开发一个R的项目，而伴随着开发的深入，测试就成了家常便饭。但是往往很多统计算法涉及到比较大的计算，比方说missing data的模型，比方说具有多层次结构的模型。测试这些R程序需要花费大量的时间。当然了如果你的程序规模很小，计算量也不大的情况下，大可不必杀鸡用牛刀，但是一般情况下，再小的一段统计算法，如果需要做一系列的simulation或是case study的话，测试都会花很多的时间。比如，如果你有几个实现了推广的ROC模型的R函数，或者是一个包含这些函数的R包，没有人能说这个程序真正管用，你要测试，那你就需要做simulation和case study。simulation简单的就是模拟出一系列预先设定模型参数的数据，让目标模型去fit，然后比较结果。case study则可以做各式各样的比较研究，特定数据的实例分析等等。好了，废话不多说，总之这篇文章就是告诉你怎么样方便地测试，随时随地想测就测。 其实核心就是调用Rterm，包含的主要文件如下: （1）创建文件：Rbatch.bat cd test date /T >> logout.txt time /T >> logout.txt cd.. C:\Progra~1\R\R-2.8.0\bin\Rterm.exe --vanilla <config.R> testout.txt （2）创建文件：config.R #如果你已经打好了R包，一下的代码可以减少，只要load你的包就好了。 #下面的版本是你还没有打包，但是调用了C++动态链接库，测试数据，和函数源代码需要source # PATH #--------------------------------------------------------- Path.autorun <- "C:\\autorun\\" Path.rdll <- "C:\\WORK\\Project\\src\\C++\\bin\\" Path.data <- "C:\\WORK\\Project\\Datasets\\Data\\R\\" Path.funs <- "C:\\WORK\\Project\\src\\R\\" Path.dump <- "C:\\autorun\\results\\" #--------------------------------------------------------- dyn.load(paste(Path.rdll,"RDLL.dll",sep="")) source(paste(Path.autorun,"source.data.R",sep="")) source(paste(Path.autorun,"source.funs.R",sep="")) source(paste(Path.autorun,"source.test.R",sep="")) source(paste(Path.autorun,"testProjectR.R",sep="")) source.data(Path.data) source.funs(Path.funs) testProjectR(paste(Path.autorun,"test\\",sep=""),TEST=0) （3）创建测试框架函数testProjectR in testProjectR.R 文件： testProjectR <- function(loc=getenv("STEST"), TEST=c(0,1,2)) { ## loc: the directory of the Project testing suite ## TEST: (0: the default loop test, 1: test buglist only, 2: loop test...
阅读全文 »