- R的极客理想:工具篇
- 张丹
- 15210字
- 2023-07-14 21:01:17
第1章
R语言基础包
本章主要介绍了为什么要学习R语言,R语言软件的安装,R语言的开发工具,以及R语言中常用的几个软件包,以帮助读者快速了解R语言的工具包,激发读者对R语言的学习兴趣。
1.1 R是最值得学习的编程语言
问题
为什么要学R语言?
引言
如果要问在Node、Lua、Python、Ruby和R这5种语言中,哪个语言在2014年的应用前景会更好,我会毫不犹豫地选择R,而且我认为R语言不仅是2014年,也是以后更长一段时间内的明星。在本书开篇,我们就谈谈为什么R语言是最值得学习的编程语言。
1.1.1 我的编程背景
本人是程序员、架构师,从编程入门到今天,一直深信着Java是改变世界的语言,Java已经做到了,而且一直很辉煌。但当Java越来越强大,覆盖领域越来越多,变得无所不能的时候,它反而不够专业,这就给了其他语言发展的机会。
我已使用Java语言11年,R语言3年,Node 1年,对于这个问题 “哪个语言在2014年的应用前景会更好”,我选择R语言。
1.1.2 为什么我会选择R
从下面的9个方面来说明我选择R的原因。
R的基因
R的发展
R的社区和资源
R的哲学
R的使用者
R的语法
R的思维模式
R解决的问题
R的不足
1.R的基因
1992年,新西兰奥克兰大学的Ross Ihaka和Robert Gentleman 两位统计学家,为了方便教授初等统计课程,发明了一种编程语言,因为他们名字的首字母都是R,于是R便成为这门语言的名称。
从开始学习R语言,我就开始了知识的跨界思考。统计学基于概率论,概率论又基于数学,用计算机的方式编程,同时解决某个领域的实际问题。多种学科知识的交集,决定着我们解决问题的能力。统计的基因,让R语言与众不同!
2.R的发展
R一直在小众领域成长着,最早也只有统计学家在用,主要用来代替SAS做统计计算。然而时代在进步,随着大数据的爆发,R终于在这一波浪潮中被工业界所发现。然后,越来越多有工程背景的人加入到这个圈子,对R的计算引擎、性能以及各种程序包进行改进和升级,让R获得了新生。
我们现在用到的R语言软件,已经越来越接近工业软件的标准了。由工程师推动的R的发展,其速度远远地超过了由统计学家推动的发展。随着人们对数据分析需求的进一步增加,R会以更快的速度继续发展,将成为免费的、开源的数据分析软件的代名词。
3.R的社区和资源
R的发展离不开R的各种社区支持,尤其是R的官方社区支持。在R的官方网站中,我们可以下载到R语言软件、R的第三方软件包和R的其他支持软件。当然,我不得不承认R的官方网站(http://www.r-project.org/)从Web页上看起来太简陋了,稍微调整一下CSS样式表,都会比现在好看很多。也许这种简单、无修饰也是统计学家的基因吧。R语言的社区资源同其他语言一样丰富,除了官方社区,还有开发者论坛(http://r.789695.n4.nabble.com/)、R-Journal列表(http://journal.r-project.org/)、软件包列表、R语言图书列表以及R用户组等。
R是自由软件,因此开发者可以开发自己的软件包,封装自己的功能,然后在CRAN(http://cran.rstudio.com/)上面发布。截止到2014年2月,共有5236个R包在CRAN上面发布。可能很多人会说只有5236个包,数量太少了。这是因为CRAN是需要提交申请的,每个包都必须经过R语言小组审核、检查后才会发布出来,而且审核非常严格。高质量是发布一个新的R包的基本要求。由于CRAN过于严格的审查,让很多开发者选择在RForge( https://r-forge.r-project.org/ )上发布R包,还有些R包是基于Github发布的,我也在Github上面发布了自己的R包:https://github.com/bsspirit/chinaWeather 。
下面列出与R语言相关的主要社区和资源。
R官方网站:http://www.r-project.org/
R开发者论坛:http://r.789695.n4.nabble.com/
CRAN:http://cran.rstudio.com/
RForge:https://r-forge.r-project.org/
R新闻和博客:http://www.r-bloggers.com/
统计之都:http://cos.name/
4.R的哲学
每种语言都有自己的设计理念和哲学,而我体会的R的哲学就是“静下心做事情”。R不需要很长的代码,也不需要设计模式。一个函数调用,传几个参数,就能实现一个复杂的统计模型。我们需要思考的是用什么模型、传什么参数,而不是怎么进行程序设计。我们可能会用R实现“从一个数学公式,变成一个统计模型”的过程,我们也可能会考虑“如何让一个分类器结果更准确”,但我们不必思考一个算法的“时间复杂度是多少,空间复杂度是多少”。
R的哲学,可以让你把数学和统计学的知识,变成计算模型,这也是R的基因所决定的。
5.R的使用者
R语言早期主要是学术界统计学家在用,后来逐渐被其他很多领域的学者所用。R语言有各种不同的应用领域,包括统计分析、应用数学、计量经济、金融分析、财经分析、人文科学、数据挖掘、人工智能、生物信息学、生物制药、全球地理科学、数据可视化等。
近几年,由互联网引发的大数据革命让工业界的人开始认识R,加入R。当越来越多的有工程背景的人加入到R语言使用者的队伍后,R才开始向着全领域发展,逐步实现工业化的要求。现在,R已不仅仅是学术界的语言,它还是工业界必备的语言。
下面列出一些推动R语言在工业界发展的R包。
RevolutionAnalytics公司的RHadoop产品,让R可以直接调用Hadoop集群资源。
RStudio公司的RStudio产品,给了我们对编辑软件新的认识。
RMySQL、ROracle、RJDBC打通了R和数据库之间的访问通道。
rmongodb、rredis、RHive、rHBase、RCassandra打通了R和NoSQL数据库之间的访问通道。
Rmpi、snow打通了单机多核并行计算的通道。
Rserve、rwebsocket 打通了R语言的跨平台通信的通道。
6.R的语法
R是面向对象语言,语法如同Python。但R的语法很自由,很多函数的名字看起来都是那么随意,这也是R的哲学的一部分吧!例如,看到如下这样的赋值语法,有其他语言基础的程序员,肯定会崩溃的。
> a<-c(1,2,3,4)->b > a [1] 1 2 3 4 > b [1] 1 2 3 4
随机取正态分布N(0,1)的10个数,又是这么的简单。
> rnorm(10) [1] -0.694541401 1.877780959 -0.178608091 0.004362026 [5] 0.836891967 1.794961298 0.115284187 0.155175219 [9] 0.464028612 -0.842569561
用R画鸢尾花的数据集的散点图,有非常好的可视化效果。
> data(iris) #加载数据集 > head(iris) #查看前6行数据集 Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa > plot(iris) # 画图
输出结果见图1-1。
图1-1 鸢尾花的数据集的散点图
正是因为R自由的哲学,让R的语法独特而简洁,我已经喜欢上这种哲学了。
7.R的思维模式
R语言让我跳出了原有的思维定式。使用R语言,我们应该像统计学家那样思考问题,而不是拘泥于程序员的思维模式。统计学家的思维模式,是先考虑为什么,再考虑做什么。而程序员的思维模式,是直接考虑怎么做,等有了结果再考虑为什么。
R语言是直接面向数据的语言。在我们的日常生活中,无论做什么事情都会产生数据,上网有浏览数据,买东西有消费数据,就算什么都不干,也会受大气PM2.5的影响,有空气污染指数数据。利用R语言,我可以直接分析这些数据。面向什么业务,就分析什么数据,不需要从产品经理向程序员的角色转换,不需要考虑有什么功能,更不需要考虑程序设计的事。跳出程序员的思维模式,我们所能认知的东西会更多,于是也能找到更适合自己的定位。
8.R解决的问题
当数据成为生产资料的时候,R就是为人们能运用生产资料创造价值的生产工具,R语言主要解决的是数据的问题。整个人类文明所获得的全部数据中,有90%以上是自互联网诞生以来产生的;当Hadoop帮助人们解决了大数据存储的问题后,如何发现数据的价值,则成为当前最火的话题。R语言具有强大的统计分析能力,这就让它成为数据分析最好的工具。所以,R要解决的问题,就是如何挖掘数据价值的问题。
9.R的不足
尽管前面说了R的各种优点,但我们依然不能说R就是完美无缺的,因为R也有很多不足。具体来说,R的缺点有下面5个。
R软件是统计学家编写的,并不如软件工程师编写的软件那么健壮。
R软件的性能,存在一些问题。
R语言很自由,语法命名不太规范,需要花时间熟悉。
R语言的内核编程,要比普通的R包使用,难度大得多。
R语言结合了很多数学、概率、统计的基础知识,学起来有一定门槛。
R的这些不足,都是可以克服的。当有更多有工程背景的人加入的时候,R语言会比现在更强大,会帮助使用者创造更多的价值。
1.1.3 R的应用前景
R可以做所有SAS能做的事情。SAS系统全称为Statistics Analysis System,是国际上最知名的商业分析软件工具之一。SAS用于决策支持的大型集成信息系统,其重要组成部分和核心功能是统计分析功能。在数据处理和统计分析领域,SAS系统被誉为国际上的标准软件系统,堪称统计软件界的巨无霸。
R和SAS处于完全的竞争的关系中,R的免费和开放,让R有着更广阔的应用前景。下面给出当今R应用最热门的领域。
统计分析:统计分布、假设检验、统计建模。
金融分析:量化策略、投资组合、风险控制、时间序列、波动率。
数据挖掘:数据挖掘算法、数据建模、机器学习。
互联网:推荐系统、消费预测、社交网络。
生物信息学:DNA分析、物种分析。
生物制药:生存分析、制药过程管理。
全球地理科学:天气、气候、遥感数据。
数据可视化:静态图、可交互的动态图、社交图、地图、热图、与各种JavaScript库的集成。
本书会介绍R语言在统计分析、金融分析、数据挖掘、推荐系统、社交网络等领域的应用。R有着非常广阔的应用前景,而且R也将成为新一代的最有能力创造价值的工具。
1.1.4 时代赋予R的任务
R语言是在大数据时代被工业界了解和认识的语言,R语言被时代赋予了挖掘数据价值、发现数据规律以及创造数据财富的任务。R语言也是帮助人们发挥智慧和创造力的最好的生产工具,因此我们不仅要学好R语言,还要用好R语言,为社会注入更多的创新的生产力。
总而言之,在这5种语言中,R是最特殊的,R被赋予了与其他语言不同的使命。R的基因决定了R将成为2014年,也可能是以后更长一段时间的明星。因此我认为“R是最值得学习的编程语言”。不论你正在读书,还是已经工作,掌握R语言这个工具并找最适合自己的位置将会前途无量。
1.2 R的历史版本安装
问题
在Linux Ubuntu上,如何安装不同版本的R?
引言
R语言已进入到了3.0的时代,但有些第三方的R包还处于2.15的状态,没有升级,如RHadoop等。我们要用这些R包的时候,就需要指定版本的R软件。对于Windows来说,这是很简单的操作,只要安装不同的(.exe)文件就行了;对于Linux系统来说,就不那么容易了,需要我们手动进行配置。不熟悉Linux系统的同学,在这里就很容易卡住。所以,本节就讲一下如何在Linux Ubuntu系统中安装R语言软件包的指定版本。
1.2.1 R在Windows中安装
通过R的官方网站(http://cran.r-project.org/) ,我们可以下载Linux、MacOS、Windows系统的R语言安装包。R在Windows系统中安装非常简单,下载可执行文件(.exe),双击即可进行安装。安装后就能运行R语言的界面,如图1-2所示。
图1-2 R在Windows系统中的安装界面
1.2.2 R在Linux Ubuntu中安装
本书使用的Linux是Ubuntu 12.04.2 LTS 64bit的系统,安装R语言软件包可以通过Ubuntu的apt-get工具进行安装。下面就介绍在Linux Ubuntu中安装R语言的过程。
~ R # 检查R是否已安装 The program 'R' is currently not installed.You can install it by typing: sudo apt-get install r-base-core ~ sudo apt-get install r-base-core # 根据提示安装R语言软件包 ~ R --version # 检查R的版本 R version 2.14.1 (2011-12-22) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-pc-linux-gnu (64-bit)
前面的检查结果表明,我们安装的是R的默认版本,即2.14.1版,这与本书中R的版本是不符的,接下来我们希望安装最新版本R的软件包。
1.2.3 R的最新版本安装
首先,删除Linux Ubuntu系统中原有的R软件包,代码如下:
~ sudo apt-get autoremove r-base-core # 删除系统中原有的R软件包
接下来,找到一个Ubuntu的软件源镜像(http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu/ ),Linux Ubuntu 12.04对应的名字是precise,进入到precise/目录,找到r-base-core相关的文件,发现有多个R的版本。把这个软件源,增加到apt的sources.list文件中,代码如下:
~ sudo sh -c "echo deb http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu precise/ >>/etc/apt/sources.list" # 在sources.list文件最下面,新加一行 ~ sudo apt-get update # 更新源 ~ sudo apt-get install r-base-core # 再次安装R语言软件包 ~ R –version # 检查R的版本 R version 3.0.3 (2014-03-06) -- "Warm Puppy" Copyright (C) 2014 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit)
这时我们就安装了最新的R语言版本——3.0.3版。
1.2.4 R的指定版本安装
由于本书中的例子,大部分是基于3.0.1完成的,而RHadoop的例子是基于2.15.3完成的,因此我们还需要指定R的安装版本。
1.安装R的2.15.3版本
~ sudo apt-get autoremove r-base-core # 删除系统中原有的R软件包 ~ sudo apt-get install r-base-core=2.15.3-1precise0precise1 # 安装R的2.15.3版本 ~ R –version # 检查R语言软件包版本 R version 2.15.3 (2013-03-01) -- "Security Blanket" Copyright (C) 2013 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: x86_64-pc-linux-gnu (64-bit)
2.安装R的3.0.1版本
~ sudo apt-get autoremove r-base-core # 删除系统中原有的R软件包 ~ sudo apt-get install r-base-core=3.0.1-6precise0 # 安装R的3.0.1版本 ~ R –version # 检查R语言软件包版本 R version 3.0.1 (2013-05-16) -- "Good Sport" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64-pc-linux-gnu (64-bit)
这样我们就可以很方便地指定安装不同版本的R的语言软件包,满足不同的应用需求!
1.3 fortunes 记录R语言的大智慧
问题
如何才能更深入地了解R, 它的起源、成长、经历是怎样的?
引言
R语言是在大数据“火”起来后,映入人们眼球的。但R语言的社区已经存在很多年,我们并不知道,R语言在很长的历史时期,有着什么样的智慧。不过,正有人悄悄地记录着R语言的大智慧。
1.3.1 fortunes介绍
fortunes库是一个R语言的语录集,截至2013年12月14日,一共总结了360条R-help的留言。这些都是R语言智慧的精华,让R语言的后辈使用者,可以更了解R语言的本身,了解R的精神。
1.3.2 fortunes安装
本节使用的系统环境是:
Linux: Ubuntu 12.04.2 LTS 64bit
R: 3.0.1 x86_64-pc-linux-gnu
注 fortunes同时支持Windows 7环境和Linux环境。
Fortunes的安装过程如下。
~ R # 启动R程序 > install.packages("fortunes") # 安装fortunes包 > library(fortunes) # 加载fortunes包 > ?fortunes # 查看帮助
1.3.3 fortunes包的使用
fortunes包的使用非常简单,只有一个函数fortune()。
> fortune() # 随机查看一条语录 Barry Rowlingson: Your grid above has 8*6 = 42 points. (That was a subtle Hitchhikers Guide To The Galaxy reference there, honest, and not a stupid dumb multiplication mistake on my part after working four 18-hour days on the trot...) Peter Dalgaard: [...] Don't panic, just throw yourself at the ground and miss. -- Barry Rowlingson and Peter Dalgaard R-help (March 2004) > fortune(108) # 指定查看一条语录 Actually, I see it as part of my job to inflict R on people who are perfectly happy to have never heard of it.Happiness doesn't equal proficient and efficient.In some cases the proficiency of a person serves a greater good than their momentary happiness. -- Patrick Burns R-help (April 2005)
完整的语录下载地址是http://cran.r-project.org/web/packages/fortunes/vignettes/fortunes.pdf。静下心来阅读这些智慧精华就能更了解R语言本身。想用好一门语言,就需要更深入地了解它。
1.4 formatR 代码自动化排版
问题
如何写出让别人看得懂,且符合规范的代码呢?
引言
新手写的代码,大都不注重代码规范,以为实现功能就算完成了。这种代码不仅别人不愿意读,过几个月后再看自己都会觉得很烂。不仅仅新手如此,大多数程序员写的代码都没有考虑如何让别人看着更方便。程序员最痛苦的事情,不是每天加班写程序,而是每天加班读懂别人写的程序。最后,有人实在忍受不了其他人的丑陋代码,便开始制定代码编程规范,又有人去实现代码的自动化排版工具。formatR就是这样的一个R语言自动化排版的工具。
1.4.1 formatR介绍
formatR包是一个实用的包,提供了R代码格式化功能,可以自动设置空格、缩进、换行等代码格式,让代码看起来更友好。formatR包中的API中主要有下面5个函数。
tidy.source: 对代码进行格式化
tidy.eval: 输出格式化后的R代码和运行结果
usage: 格式化函数定义,并按指定宽度输出
tidy.gui: 一个GUI工具,支持编辑并格式化R代码
tidy.dir: 对某个目录下,所有R脚本进行格式化
1.4.2 formatR安装
本节使用的系统环境是:
Win7 64bit
R: 3.0.1 x86_64-w64-mingw32/x64 b4bit
注 formatR同时支持Windows 7环境和Linux环境。formatR的安装过程如下:
~ R # 启动R程序 > install.packages("formatR") # 安装formatR包 library(formatR) # formatR加载
1.4.3 formatR的使用
1.字符串格式化
tidy.source()函数,以字符串作为输入参数,对代码格式化。
> tidy.source(text = c("{if(TRUE)1 else 2; if(FALSE){1+1", "## comments", "} else 2}")) { if (TRUE) 1 else 2 if (FALSE) { 1 + 1 ## comments } else 2 }
通过执行tidy.source()函数,把代码进行了重新格式化,让我们一眼就可以看得懂。
2.文件格式化
messy.R是一个不太规范的R程序文件。我们读入这个文件,然后通过tidy.source()函数,以文件对象作为输入参数,进行代码格式化。
> messy = system.file("format", "messy.R", package = "formatR") > messy [1] "C:/Program Files/R/R-3.0.1/library/formatR/format/messy.R"
messy.R 的原始代码输出:
> src = readLines(messy) > cat(src,sep="\n") # a single line of comments is preserved 1+1 if(TRUE){ x=1 # inline comments }else{ x=2;print('Oh no...ask the right bracket to go away!')} 1*3 # one space before this comment will become two! 2+2+2 # 'short comments' lm(y~x1+x2, data=data.frame(y=rnorm(100),x1=rnorm(100),x2=rnorm(100))) ### only'single quotes' are allowed in comments 1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1+1 ## comments after a long line'a character string with \t in it' ## here is a long long long long long long long long long long long long long long long long long long long long comment
格式化后的代码输出:
> tidy.source(messy) # a single line of comments is preserved 1 + 1 if (TRUE) { x = 1 # inline comments } else { x = 2 print("Oh no...ask the right bracket to go away!") } 1 * 3 # one space before this comment will become two! 2 + 2 + 2 # 'short comments' lm(y ~ x1 + x2, data = data.frame(y = rnorm(100), x1 = rnorm(100), x2 =rnorm(100))) ### only 'single quotes' are allowed in comments 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 +1 + 1 + 1 ## comments after a long line "a character string with \t in it" ## here is a long long long long long long long long long long long long long long long long ##long long long long comment
可以看出,格式化后的输出,经过了空格、缩进、换行、注释等处理,代码可读性就增强了。
3.格式化并输出R脚本文件
新建R脚本文件demo.r。
~ vi demo.r a<-1+1;a;matrix(rnorm(10),5); if(a>2) { b=c('11',832);"#a>2";} else print('a is invalid!!')
格式化demo.r。
> x = "demo.r" > tidy.source(x) a <- 1 + 1 a matrix(rnorm(10), 5) if (a > 2) { b = c("11", 832) "#a>2" } else print("a is invalid!!")
输出格式化结果到文件demo2.r,如图1-3所示。
> f="demo2.r" > tidy.source(x, keep.blank.line = TRUE, file = f) > file.show(f)
图1-3 输出格式化结果到文件
4.输出格式化代码和运行结果
使用tidy.eval()函数,以字符串形式,执行R脚本:
> tidy.eval(text = c("a<-1+1;a", "matrix(rnorm(10),5)")) a <- 1 + 1 a ## [1] 2 matrix(rnorm(10), 5) ## [,1] [,2] ## [1,] 0.65050729 0.1725221 ## [2,] 0.05174598 0.3434398 ## [3,] -0.91056310 0.1138733 ## [4,] 0.18131010 -0.7286614 ## [5,] 0.40811952 1.8288346
这样直接在当前的运行环境中,就输出了代码和运行结果。
5.格式化函数定义
通过usage()函数可以只打印出函数定义,跳过函数细节。以var()函数为例,输入var,默认会打印出一个函数细节。
> var function (x, y = NULL, na.rm = FALSE, use) { if (missing(use)) use <- if (na.rm) "na.or.complete" else "everything" na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs","everything", "na.or.complete")) if (is.na(na.method)) stop("invalid 'use' argument") if (is.data.frame(x)) x <- as.matrix(x) else stopifnot(is.atomic(x)) if (is.data.frame(y)) y <- as.matrix(y) else stopifnot(is.atomic(y)) .Call(C_cov, x, y, na.method, FALSE) } <bytecode: 0x0000000008fad030> <environment: namespace:stats> > usage(var) # 通过usage,只打印函数定义 var(x, y = NULL, na.rm = FALSE, use)
有时候函数定义也很长,比如lm()函数,通过usage的width参数可以控制函数的显示宽度。
> usage(lm) lm(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...) > usage(lm,width=30) # usage的width参数,控制函数的显示宽度 lm(formula, data, subset, weights, na.action, method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts = NULL, offset, ...)
6.GUI工具
tidy.gui()函数是一个GUI的工具,可以在界面上编辑并格式化R代码。首先安装gWidgetsRGtk2库:
> install.packages("gWidgetsRGtk2") also installing the dependencies 'RGtk2', 'gWidgets'
打开GUI控制台:
> library("gWidgetsRGtk2") > g = tidy.gui()
我们输入一段不太好看的代码,如图1-4所示。
图1-4 tidy的GUI
点击“转换”,结果如图1-5所示,可以看到,在GUI的编辑器中,R语言的代码被格式化了!
图1-5 格式化后的代码
7.格式化目录中的文件
tidy.dir()函数可以批量格式化文件,对目录中的所有文件进行格式化操作。下面新建目录:dir,并在目录dir中新建两个R脚本文件:dir.r, dir2.r。
~ mkdir dir # 新建目录 dir ~ cd dir ~ vi dir.r # 用vi新建文件dir.r a<-1+1;a;matrix(rnorm(10),5); ~ vi dir2.r if(a>2) { b=c('11',832);"#a>2";} else print('a is invalid!!')
执行tidy.dir:
> tidy.dir(path="dir") tidying dir/dir.r tidying dir/dir2.r
分别查看dir.r和dir2.r:
~ vi dir.r a <- 1 + 1 a matrix(rnorm(10), 5) ~ vi dir2.r if (a > 2) { b = c("11", 832) "#a>2" } else print("a is invalid!!")
我们发现不规则的代码,已经被格式化了!
1.4.4 formatR的源代码解析
通过上面的使用,我们不难发现,formatR包的核心函数就是tidy.source()函数,从Github上面找到源代码:https://github.com/yihui/formatR/blob/master/R/tidy.R。我将在代码中增加注释:
tidy.source = function( source = 'clipboard', keep.comment = getOption('keep.comment', TRUE), keep.blank.line = getOption('keep.blank.line', TRUE), replace.assign = getOption('replace.assign', FALSE), left.brace.newline = getOption('left.brace.newline', FALSE), reindent.spaces = getOption('reindent.spaces', 4), output = TRUE, text = NULL, width.cutoff = getOption('width'), ... ) { if (is.null(text)) { # 判断输入来源为剪贴板 if (source == 'clipboard' && Sys.info()['sysname'] == 'Darwin') { source = pipe('pbpaste') } } else { # 判断输入来源为字符串 source = textConnection(text); on.exit(close(source)) } text = readLines(source, warn = FALSE) # 按行读取来源数据 if (length(text) == 0L || all(grepl('^\\s*$', text))) { # 文件行数判断 if (output) cat('\n', ...) return(list(text.tidy = text, text.mask = text)) } if (keep.blank.line && R3) { # 空行处理 one = paste(text, collapse = '/n') # record how many line breaks before/after n1 = attr(regexpr('^/n*', one), 'match.length') n2 = attr(regexpr('\n*$', one), 'match.length') } if (keep.comment) text = mask_comments(text, width.cutoff, keep.blank.line) # 注释处理 text.mask = tidy_block(text, width.cutoff, replace.assign && length(grep('=', text))) # 把输入的R代码,先转成表达式,再转回字符串。用来实现对每个语句的截取 text.tidy = if (keep.comment) unmask.source(text.mask) else text.mask # 对注释排版 text.tidy = reindent_lines(text.tidy, reindent.spaces) # 重新定位缩进 if (left.brace.newline) text.tidy = move_leftbrace(text.tidy) # 扩号换行 if (keep.blank.line && R3) text.tidy = c(rep('', n1), text.tidy, rep('', n2)) # 增加首尾空行 if (output) cat(paste(text.tidy, collapse = '\n'), '\n', ...) # 在console打印格式化后的结果 invisible(list(text.tidy = text.tidy, text.mask = text.mask)) # 返回,但不打印结果 }
1.4.5 源代码中的Bug
在读源代码的过程中,我发现有一个小问题,即在R 3.0.1版本,没有对向右赋值操作(->)进行处理。我已经就这个问题给作者提Bug了,参见https://github.com/yihui/formatR/ issues/31。Bug测试代码如下:
> c('11',832)->x2 > x2 [1] "11" "832" > tidy.source(text="c('11',832)->x2") # 格式化代码 c("11", 832) <- x2 > tidy.eval(text="c('11',832)->x2") c("11", 832) <- x2 Error in eval(expr, envir, enclos) : object 'x2' not found
Bug已修复。作者回复:“这个问题已经在R 3.0.2中修正了。”
> formatR::tidy.source(text="c('11',832)->x2") # 格式化代码 x2 <- c("11", 832) > sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] formatR_0.10.3
formatR包提供的功能非常实用,特别是读别人写的不规范的代码的时候。建议各IDE厂商能把formatR作为标准的格式化工具直接嵌入编辑器的工具中。让我们把阅读别人的代码,也变成一件快乐的事情吧。
1.5 多人在线协作R开发RStudio Server
问题
R语言开发,哪个工具最好用?
引言
RStudio是R语言开发中的利器,是最好用的R语言IDE集成环境。RStudio Server更是利器中的神器。不仅提供了Web的功能,可以安装到远程服务器上,通过Web进行访问,还支持多用户的协作开发。如此神器,快来动手试一下吧。
1.5.1 RStudio和RStudio Server
RStudio 是一个强大的、免费的、开源的R语言集成开发环境的应用软件,可以安装在Windows、Linux和Mac 不同操作系统上。RStudio Server 是一个基于Web访问的RStudio云端开发环境,需要安装在Linux服务器上面,支持多用户远程访问使用。
1.5.2 RStudio Server安装
本文使用的系统环境是:
Linux: Ubuntu Server 12.04.2 LTS 64bit
R: 3.0.1 x86_64-pc-linux-gnu
IP: 192.168.1.13
注 RStudio Server只支持Linux系统环境。下载最新版本RStudio Server的地址是http://www.rstudio.com/ide/download/server.html。
在Linux Ubuntu环境中,下载并安装64位的Rstudio Server:
~ sudo apt-get install gdebi-core ~ sudo apt-get install libapparmor1 # Required only for Ubuntu, not Debian ~ wget http://download2.rstudio.org/rstudio-server-0.97.551-amd64.deb ~ sudo gdebi rstudio-server-0.97.551-amd64.deb
安装后,RStudio Server会自动启动运行。
~ ps -aux|grep rstudio-server # 查看RStudio Server运行进程 998 2914 0.0 0.1 192884 2568 ? Ssl 10:40 0:00 /usr/lib/ rstudio-server/bin/rserver
可以看到,RStudio Server的服务已启动,8787端口被打开。
1.5.3 RStudio Server使用
通过浏览器,我们访问RStudio Server: http://192.168.1.13:8787,IP地址为RStudio Server服务器的地址,如图1-6所示。
图1-6 RStudio Server登录界面
RStudio Server 登录需要用Linux 系统的用户账号。如果想增加或减少用户,直接对Linux系统用户进行操作就可以了。我的环境中用户登录,用户名是conan,密码是conan111。登录之后看到的界面如图1-7所示。
图1-7 RStudio Server的Web界面
1.RStudio Server 的系统配置
RStudio Server主要有两个配置文件,默认文件不存在。
/etc/rstudio/rserver.conf /etc/rstudio/rsession.conf
设置端口和ip控制:
~ vi /etc/rstudio/rserver.conf www-port=8080 # 监听端口 www-address=127.0.0.1 # 允许访问的IP地址,默认为0.0.0.0
重启Rstudio Server服务器,配置生效:
~ sudo rstudio-server restart
会话配置管理:
~ vi /etc/rstudio/rsession.conf session-timeout-minutes=30 # 会话超时时间 r-cran-repos=http://ftp.ctex.org/mirrors/CRAN/ # 设置CRAN资源库
2.RStudio Server 的系统管理
启动、停止、重启 RStudio Server服务器的命令如下:
~ sudo rstudio-server start # 启动 ~ sudo rstudio-server stop # 停止 ~ sudo rstudio-server restart # 重启
查看运行中的R进程:
~ sudo rstudio-server active-sessions PID TIME COMMAND 6817 00:00:03 /usr/lib/rstudio-server/bin/rsession -u zd
指定PID, 停止运行中的R进程:
~ sudo rstudio-server suspend-session 6817 ~ sudo rstudio-server active-sessions # 再次查看进程 PID TIME COMMAND
停止所有运行中的R进程:
~ sudo rstudio-server suspend-all
强制停止运行中的R进程,此操作优先级最高,立刻执行。
~ sudo rstudio-server force-suspend-session <pid> ~ sudo rstudio-server force-suspend-all
RStudio Server临时下线,不允许Web访问,并给用户友好的错误提示:
~ sudo rstudio-server offline rstudio-server start/running, process 6880
RStudio Server上线:
~ sudo rstudio-server online rstudio-server start/running, process 6908
RStudio Server的其他操作和单机版的RStudio一样。
1.5.4 RStudio Server多人协作
1.增加新用户和新用户组
~ sudo groupadd hadoop # 创建Hadoop用户组 ~ sudo useradd hadoop -g hadoop # 创建Hadoop用户并加入到Hadoop用户组 ~ sudo passwd hadoop # 设置Hadoop用户的密码 ~ sudo adduser hadoop sudo # 增加Hadoop用户到sudo组 ~ sudo mkdir /home/hadoop # 创建Hadoop用户的home目录 ~ sudo chown -R hadoop:hadoop /home/hadoop # 给/home/hadoop目录及子目录,设置用户权限
以Hadoop用户登录,检查用户是否设置成功
~ ssh hadoop@localhost # 通过 ssh 远程登录 ~ bash ~ pwd # 查看登录后的访问目录 /home/hadoop
新打开浏览器窗口通过Hadoop账号登录 ,如图1-8所示。
图1-8 以Hadoop账号登录
2.Git代码共享
首先安装Git:
~ sudo apt-get install git ~ ssh-keygen -t rsa #生成rsa密钥对 ~ cat /home/conan/.ssh/id_rsa.pub #查看公钥 ssh-rsa AAAAB3NzaC1yc2EAAAADAQABAAABAQDMmnFyZe2RHpXaGmENdH9kSyDyVzRas4GtRwMNx+qQ 4QsB8xVTrIbFayG2ilt+P8UUkVYO0qtUJIaLRjGy/SvQzzL7JKX12+VyYoKTfKvZZnANJ414d6oZpbDw sC0Z7JARcWsFyTW1KxOMyesmzNNdB+F3bYN9sYNiTkOeVNVYmEQ8aXywn4kcljBhVpT8PbuHl5eadSLt 5zpN6bcX7tlquuTlRpLi1e4K+8jQo67H54FuDyrPLUYtVaiTNT/xWN6IU+DQ9CbfykJ0hrfDU1d1LiLQ 4K2Fdg+vcKtB7Wxez2wKjsxb4Cb8TLSbXdIKEwSOFooINw25g/Aamv/nVvW1 conan@conan-deskop
接下来,我们需要把本地项目上传到Github。首先在Github上创建一个新的项目rstudio-demo,地址为https://github.com/bsspirit/rstudio-demo,通过下面的操作上传本地目录到rstudio-demo项目。
~ mkdir /home/conan/R/github # 创建rstudio-demo项目目录 ~ cd /home/conan/R/github ~ git init # 初始化Git ~ git add # 增加当前目录及子目录到本地Git库 ~ git commit -m 'first comment' # 在本地Git库提交 ~ git remote add origin git@github.com:bsspirit/rstudio-demo.git # 绑定当前目录和github的项目 ~ git push -u origin master # 上传本地Git库中的代码到Github
打开RStudio设置到/home/conan/R/github目录,tools–>version control –> project setup,如图1-9所示。
图1-9 在RStudio中,配置Github地址
在RStudio中修改sayHello.r的代码:
sayHello<-function(name){ print(paste("hello",name)) } sayHello("Conan") sayHello("World")
点击tools–>version control–> commit提交,如图1-10所示。
图1-10 通过RStudio进行Git操作
上传到Github,只需要点击tools–>version control–> push,如图1-11所示。
图1-11 在Github上查看提交的操作
RStudio有如此强大的功能,极大地降低了编程的门槛。还没有用过的同学,赶紧去体验一把极客的感觉吧!
1.6 R和JSON的傻瓜式编程
问题
如何让R语言的数据类型转换成JSON数据类型?
引言
JSON作为一种轻量级数据格式,被大量地应用在各种程序环境中。JSON(JavaScript Object Notation)是JavaScript的内嵌的标准对象,同时也是MongoDB的表结构存储类型。JSON是半结构化的,可以表达出丰富的文档含义。JSON文档比XML文档要少很多,更适合于网络传输。早期R语言编程很少会用到JSON,但随着R语言的壮大,R也在伸向各种领域,JSON就是与其他领域的一个交点。如何让R语言傻瓜式转型JSON呢?请看下文介绍。
1.6.1 rjson包介绍
rjson是一个R语言与JSON进行转换的包,非常简单,支持R自身语言转型和基于C类库的转型两种方式。rjson包提供的函数只有3个,fromJSON(), newJSONParser(), toJSON()。 后面我们将介绍如何使用rjson包。本节使用的系统环境是:
Windows 7: x86_64-w64-mingw32/x64 (64-bit)
R: version 3.0.1
注 rjson同时支持Windows 7环境和Linux环境。
1.安装并加载rjson
> install.packages("rjson") # 安装rjson > library(rjson) # 加载rjson
接下来,我们进行测试。新建JSON文件fin0.json:
~ vi fin0.json { "table1": { "time": "130911", "data": { "code": [ "TF1312", "TF1403", "TF1406" ], "rt_time": [ 130911, 130911, 130911 ] } }, "table2": { "time": "130911", "data": { "contract": [ "TF1312", "TF1312", "TF1403" ], "jtid": [ 99, 65, 21 ] } } }
2.调用函数fromJSON(): 从JSON到R
从fin0.json文件中,读取JSON并解析成R语言对象。我们通常把字节或者文字格式转型为程序对象的过程叫反序列化过程,与之相反的过程,叫做序列化过程。
> json_data <- fromJSON(paste(readLines("fin0.json"), collapse="")) > json_data $table1 $table1$time [1] "130911" $table1$data $table1$data$code [1] "TF1312" "TF1403" "TF1406" $table1$data$rt_time [1] 130911 130911 130911 $table2 $table2$time [1] "130911" $table2$data $table2$data$contract [1] "TF1312" "TF1312" "TF1403" $table2$data$jtid [1] 99 65 21
检查转型后,对应R语言的数据类型:
> class(json_data) [1] "list" > class(json_data$table2) [1] "list" > class(json_data$table2$data) [1] "list" > class(json_data$table2$data$jtid) [1] "numeric" > class(json_data$table1$data$code) [1] "character"
我们看到原JSON对象转型后,除最内层外,其他都被解析为R的列表(list)类型,最内层则是基本(numeric,character)类型。在R对象结构中取JSON数据的一个叶子节点,JSON的索引路径为 json.table1.data.code[0]。
> json_data$table1$data$code [1] "TF1312" "TF1403" "TF1406" > json_data$table1$data$code[1] [1] "TF1312"
3.toJSON() 从R到JSON
把R对象转型为JSON串,这个过程叫做序列化过程。还以刚才的json_data为例。
> json_str<-toJSON(json_data) > print(json_str) [1] "{\"table1\":{\"time\":\"130911\",\"data\":{\"code\":[\"TF1312\",\"TF1403\",\"TF1406\"],\"rt_time\":[130911,130911,130911]}},\"table2\":{\"time\":\"130911\",\"data\":{\"contract\":[\"TF1312\",\"TF1312\",\"TF1403\"],\"jt id\":[99,65,21]}}}" > cat(json_str) {"table1":{"time":"130911","data":{"code":["TF1312","TF1403","TF1406"],"rt_time":[130911,130911,130911]}},"table2":{"time":"130911","data":{"contract":["TF1312","TF1312","TF1403"],"jtid":[99,65,21]}}}
我们只要使用toJSON()函数,就可以实现R对象向JSON的转型。如果用print()函数输出,结果就是带转义的输出(\");如果直接用cat函数输出,结果就是标准的JSON串格式。把JSON输出到文件fin0_out.json, 有2种方法,即writeLines()和sink()。
> writeLines(json_str, "fin0_out1.json") # writeLines方法 > sink("fin0_out2.json") # sink方法 > cat(json_str) > sink()
虽然写法不同,但是输出结果是一个样的,writeLines最后新建一个空行。
{"table1":{"time":"130911","data":{"code":["TF1312","TF1403","TF1406"],"rt_time": [130911,130911,130911]}},"table2":{"time":"130911","data":{"contract":["TF1312", "TF1312","TF1403"],"jtid":[99,65,21]}}}
4.C语言库和R语言库转型,并进行性能测试
我们对fromJSON进行性能测试:
> system.time( y <- fromJSON(json_str,method="C") ) 用户 系统 流逝 0 0 0 > system.time( y2 <- fromJSON(json_str,method = "R") ) 用户 系统 流逝 0.02 0.00 0.02 > system.time( y3 <- fromJSON(json_str) ) 用户 系统 流逝 0 0 0
我们可以看到,基于C语言库的操作比基于R语言库的要快,因为数据量很小,所以0.02并不明显。当JSON串很大的时候,这个差距就会变得相当大了。fromJSON默认使用C语言库的方法,所以我们平时处理不用加method='C'的参数。下面做toJSON的性能测试。
> system.time( y <- toJSON(json_data,method="C") ) 用户 系统 流逝 0 0 0 > system.time( y2 <- toJSON(json_data,method = "R") ) 用户 系统 流逝 0.02 0.00 0.01 > system.time( y3 <- toJSON(json_data) ) 用户 系统 流逝 0 0 0
解释同前。
1.6.2 RJSONIO包介绍
RJSONIO包,提供了2个主要的操作,把JSON串反序列化成R对象,把R对象序列化成JSON串。RJSONIO包的两个主要函数是fromJSON(), toJSON(),它还包括几个辅助函数,即asJSVars(), basicJSONHandler(), Bob(), isValidJSON(), readJSONStream()。RJSONIO包解决了rjson包序列化大对象慢的问题。RJSONIO依赖于底层的C语言类库libjson。
1.安装并加载RJSONIO
> install.packages("RJSONIO") > library(RJSONIO)
2.fromJSON() 从JSON到R
同rjson一样,测试fromJSON()函数。
> json_data <- fromJSON(paste(readLines("fin0.json"), collapse="")) > json_data $table1 $table1$time [1] "130911" $table1$data $table1$data$code [1] "TF1312" "TF1403" "TF1406" $table1$data$rt_time [1] 130911 130911 130911 $table2 $table2$time [1] "130911" $table2$data $table2$data$contract [1] "TF1312" "TF1312" "TF1403"$table2$data$jtid [1] 99 65 21
我们发现与 rjson的结果是一样,R对象除最内层外,其他都是列表(list)类型。下面取叶子节点:
> json_data$table1$data$code [1] "TF1312" "TF1403" "TF1406" > json_data$table1$data$code[1] [1] "TF1312"
3.toJSON() 从R到JSON
做toJSON的性能测试:
> json_str<-toJSON(json_data) > print(json_str) [1] "{\n\"table1\": {\n\"time\": \"130911\",\n\"data\": {\n\"code\": [\"TF1312\", \"TF1403\", \"TF1406\" ],\n\"rt_time\": [ 1.3091e+05, 1.3091e+05,1.3091e+05 ]\n}\n},\n\"table2\": {\n\"time\": \"130911\",\n\"data\": {\n\"contract\": [ \"TF1312\", \"TF1312\", \"TF1403\" ],\n\"jtid\": [99,65,21 ]\n}\n}\n}" > cat(json_str) { "table1": { "time": "130911", "data": { "code": [ "TF1312", "TF1403", "TF1406" ], "rt_time": [ 1.3091e+05, 1.3091e+05, 1.3091e+05 ] } }, "table2": { "time": "130911", "data": { "contract": [ "TF1312", "TF1312", "TF1403" ], "jtid": [99,65,21 ] } } }
toJSON的函数输出与rjson是不一样的,这个输出是格式化的。下面输出到文件:
> writeLines(json_str, "fin0_io.json")
文件结果:
{ "table1": { "time": "130911", "data": { "code": [ "TF1312", "TF1403", "TF1406" ], "rt_time": [ 1.3091e+05, 1.3091e+05, 1.3091e+05 ] } }, "table2": { "time": "130911", "data": { "contract": [ "TF1312", "TF1312", "TF1403" ], "jtid": [ 99, 65, 21 ] } } }
4.isValidJSON() 验证JSON是否合法
验证JSON的格式是否合法。
> isValidJSON(json_str) Error in file(con, "r") : cannot open the connection > isValidJSON(json_str,TRUE) #合法JSON [1] TRUE > isValidJSON(I('{"foo" : "bar"}')) #合法JSON [1] TRUE > isValidJSON(I('{foo : "bar"}')) #不合法JSON [1] FALSE
5.asJSVars() 转换为JavaScript变量格式
> cat(asJSVars( a = 1:10, myMatrix = matrix(1:15, 3, 5))) a = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ] ; myMatrix = [ [ 1, 4, 7, 10, 13 ], [ 2, 5, 8, 11, 14 ], [ 3, 6, 9, 12, 15 ] ] ;
得到两个JavaScript变量,即数组a和二维数组myMatrix.
1.6.3 自定义JSON的实现
在R语言中,我们最常用的类型是data.frame,接下来我们自己实现data.frame类型和JSON的转换。首先,把R的data.frame对象转成我们定义的JSON格式。
定义输出JSON的格式:
[ { "code": "TF1312", "rt_time": "152929", "rt_latest": 93.76, "rt_bid1": 93.76, "rt_ask1": 90.76, "rt_bsize1": 2, "rt_asize1": 100, "optionValue": -0.4, "diffValue": 0.6 } ]
定义R语言data.frame类型对象:
> df<-data.frame( code=c('TF1312','TF1310','TF1313'), rt_time=c("152929","152929","152929"), rt_latest=c(93.76,93.76,93.76), rt_bid1=c(93.76,93.76,93.76), rt_ask1=c(90.76,90.76,90.76), rt_bsize1=c(2,3,1), rt_asize1=c(100,1,11), optionValue=c(-0.4,0.2,-0.1), diffValue=c(0.6,0.6,0.5) ) > df code rt_time rt_latest rt_bid1 rt_ask1 rt_bsize1 rt_asize1 optionValue diffValue 1 TF1312 152929 93.76 93.76 90.76 2 100 -0.4 0.6 2 TF1310 152929 93.76 93.76 90.76 3 1 0.2 0.6 3 TF1313 152929 93.76 93.76 90.76 1 11 -0.1 0.5
直接使用toJSON,输出的JSON串按列组合成了数组,并不是我们想要的。
> cat(toJSON(df)) { "code": ["TF1312", "TF1310", "TF1313"], "rt_time": ["152929", "152929", "152929"], "rt_latest": [93.76, 93.76, 93.76], "rt_bid1": [93.76, 93.76, 93.76], "rt_ask1": [90.76, 90.76, 90.76], "rt_bsize1": [2, 3, 1], "rt_asize1": [100, 1, 11], "optionValue": [-0.4, 0.2, -0.1], "diffValue": [0.6, 0.6, 0.5] }
我们要对data.frame类型进行数据处理:
> library(plyr) #用plyr进行数据转换 > cat(toJSON(unname(alply(df, 1, identity)))) [ { "code": "TF1312", "rt_time": "152929", "rt_latest": 93.76, "rt_bid1": 93.76, "rt_ask1": 90.76, "rt_bsize1": 2, "rt_asize1": 100, "optionValue": -0.4, "diffValue": 0.6 }, { "code": "TF1310", "rt_time": "152929", "rt_latest": 93.76, "rt_bid1": 93.76, "rt_ask1": 90.76, "rt_bsize1": 3, "rt_asize1": 1, "optionValue": 0.2, "diffValue": 0.6 }, { "code": "TF1313", "rt_time": "152929", "rt_latest": 93.76, "rt_bid1": 93.76, "rt_ask1": 90.76, "rt_bsize1": 1, "rt_asize1": 11, "optionValue":-0.1, "diffValue": 0.5 } ]
输出的结果已经通过alply()函数做了数据变换,这正是我希望的按行输出的结果。
1.6.4 JSON性能比较
性能比较我们做2组测试,一组是rjson和RJSONIO 对大对象进行序列化(toJSON)测试,另一组是RJSONIO包序列化(toJSON) 列式输出和行式输出的测试。
1.rjson和RJSONIO 对大对象进行序列化(toJSON)测试
创建一个rjson测试脚本,在命令行运行。
> library(rjson) > df<-data.frame( + a=rep(letters,10000), + b=rnorm(260000), + c=as.factor(Sys.Date()-rep(1:260000)) + ) > system.time(rjson::toJSON(df)) 1.01 0.02 1.03 > system.time(rjson::toJSON(df)) 1.01 0.03 1.04 > system.time(rjson::toJSON(df)) 0.98 0.05 1.03
同样,再创建一个RJSONIO测试脚本,在命令行运行。
> library(RJSONIO) > df<-data.frame( + a=rep(letters,10000), + b=rnorm(260000), + c=as.factor(Sys.Date()-rep(1:260000)) + ) > system.time(RJSONIO::toJSON(df)) 2.23 0.02 2.24 > system.time(RJSONIO::toJSON(df)) 2.30 0.00 2.29 > system.time(RJSONIO::toJSON(df)) 2.25 0.01 2.26
对比结果发现,rjson的性能优于RJSONIO。
2.rjson和RJSONIO,序列化(toJSON)列式输出和行式输出的测试
创建一个rjson测试脚本,在命令行运行。
> library(rjson) > library(plyr) > df<-data.frame( + a=rep(letters,100), + b=rnorm(2600), + c=as.factor(Sys.Date()-rep(1:2600)) + ) > system.time(rjson::toJSON(df)) 0.01 0.00 0.02 > system.time(rjson::toJSON(df)) 0.01 0.00 0.02 > system.time(rjson::toJSON(unname(alply(df, 1, identity)))) 1.55 0.02 1.56 > system.time(rjson::toJSON(unname(alply(df, 1, identity)))) 0.83 0.00 0.83
同样,再创建一个RJSONIO测试脚本,在命令行运行。
> library(RJSONIO) > library(plyr) > df<-data.frame( + a=rep(letters,100), + b=rnorm(2600), + c=as.factor(Sys.Date()-rep(1:2600)) + ) > system.time(RJSONIO::toJSON(df)) 0.03 0.00 0.03 > system.time(RJSONIO::toJSON(df)) 0.04 0.00 0.03 > system.time(RJSONIO::toJSON(unname(alply(df, 1, identity)))) 2.82 0.02 2.84 > system.time(RJSONIO::toJSON(unname(alply(df, 1, identity)))) 2.06 0.00 2.06
通过测试我们发现,用toJSON直接列式输出,比行式输出要高效得多,这是因为行输出需要多一步数据处理的过程。这里说rjson比RJSONIO高效,而前面说RJSONIO包解决了rjson包序列化大对象慢的问题,似乎有些矛盾。当然,我的一次的测试不能证明这个结论的,也希望大家在有条件、有精力的情况下,自己做一些测试。
1.7 R语言的高质量图形渲染库Cairo
问题
如何让R语言画出无锯齿的高清图?
引言
R语言不仅在统计分析和数据挖掘领域计算能力强大,它在数据可视化领域也不逊于昂贵的商业软件。当然,R在可视化上强大,其背后离不开各种开源软件包的支持,Cairo就是这样一个用于矢量图形处理的类库。Cairo可以创建高质量的矢量图形(GIF、SVG、PDF、PostScript) 和位图(PNG、JPEG、TIFF),同时支持在后台程序中高质量渲染!本节将介绍Cairo在R语言中的使用。
1.7.1 Cairo介绍
Cairo 是一个用于图形绘图和渲染的免费库,支持复杂的 2D 的绘图功能,支持硬件加速。虽然,Cairo 是用C语言编写的,但提供多种语言的接口,允许其他语言直接调用,包括有 C++、C#、Java、Python、Perl、Ruby、Scheme、Smalltalk 等语言。Cairo发布的许可协议为 GNU Lesser General Public License version 2.1(LGPL) 或 Mozilla Public License 1.1(MPL)。R语言Cairo接口的官方发布页是http://www.rforge.net/Cairo/。
1.7.2 Cairo包安装
本节使用的系统环境是:
Linux: Ubuntu 12.04.2 LTS 64bit
R: 3.0.1 x86_64-pc-linux-gnu
注 caTools同时支持Windows 7环境和Linux环境。
Cairo包在Linux Ubuntu系统中的安装过程如下:
~ sudo apt-get install libcairo2-dev # Cairo的底层依赖库 ~ sudo apt-get install libxt-dev ~ R # 启动R程序 > install.packages("Cairo") # 安装Cairo包
1.7.3 Cairo使用
Cairo使用起来非常简单,和基础包grDevices中的函数对应。
CairoPNG: 对应grDevices:png()。
CairoJPEG: 对应grDevices:jpeg()。
CairoTIFF: 对应grDevices:tiff()。
CairoSVG: 对应grDevices:svg()。
CairoPDF: 对应grDevices:pdf()。
我常用的图形输出,就是png和svg。下面检查Cairo的兼容性:
> library(Cairo) # 加载Cairo包 > Cairo.capabilities() # 检查Cairo包支持的图片格式 png jpeg tiff pdf svg ps x11 win raster TRUE TRUE FALSE TRUE TRUE TRUE TRUE FALSE TRUE
从兼容性的检查结果,我们可以查看Cairo支持的图形输出格式:
支持:png、jpeg、pdf、svg、ps、x11(Linux桌面)、raster
不支持:tiff、win(win桌面)
注 如果是Windows系统,则x11为FALSE, win为TRUE。
下面比较 CairoPNG() 和 png() 的输出效果。
1.散点图
首先我们来画一个6000个点的散点图。
> x<-rnorm(6000) # 随机取6000个点坐标 > y<-rnorm(6000) > png(file="plot4.png",width=640,height=480) # png函数 > plot(x,y,col="#ff000018",pch=19,cex=2,main = "plot") > dev.off() > CairoPNG(file="Cairo4.png",width=640,height=480) # CairoPNG函数 > plot(x,y,col="#ff000018",pch=19,cex=2,main = "Cairo") > dev.off()
在当前目录,会生成2个png文件,即plot4.png和Cairo4.png,见图1-12和图1-13。
图1-12 以png()函数生成的图
图1-13 以CairoPNG()函数生成的图
SVG图形输出代码为:
> svg(file="plot-svg4.svg",width=6,height=6) > plot(x,y,col="#ff000018",pch=19,cex=2,main = "plot-svg") > dev.off() > CairoSVG(file="Cairo-svg4.svg",width=6,height=6) > plot(x,y,col="#ff000018",pch=19,cex=2,main = "Cairo-svg") > dev.off()
在当前目录,会生成2个svg文件,即 plot-svg4.svg和Cairo-svg4.svg,可以把svg的图片拖拽到浏览器中显示。
2.三维截面图
接下来,我们再画一个三维截面图,比较函数png()和CairoPNG()输出的效果。
> x <- seq(-10, 10, length= 30) > y <- x > f <- function(x,y) { r <- sqrt(x^2+y^2); 10 * sin(r)/r } > z <- outer(x, y, f) > z[is.na(z)] <- 1 > png(file="plot2.png",width=640,height=480) # PNG图 > op <- par(bg = "white", mar=c(0,2,3,0)+.1) > persp(x, y, z,theta = 30, phi = 30,expand = 0.5,col = "lightblue",ltheta =120,shade = 0.75,ticktype = "detailed",xlab = "X", ylab = "Y", zlab ="Sinc(r)",main = "Plot") > par(op) > dev.off() > CairoPNG(file="Cairo2.png",width=640,height=480) > op <- par(bg = "white", mar=c(0,2,3,0)+.1) > persp(x, y, z,theta = 30, phi = 30,expand = 0.5,col = "lightblue",ltheta =120,shade = 0.75,ticktype = "detailed",xlab = "X", ylab = "Y", zlab ="Sinc(r)",main = "Cairo") > par(op) > dev.off()
在当前目录,会生成2个png文件,即plot2.png和Cairo2.png,见图1-14和图1-15。
3.大量文字的图片
最后,我们针对包含大量文字的图片比较png()和CairoPNG()函数输出的效果。
> library(MASS) # 加载MASS包 > data(HairEyeColor) # 加载HairEyeColor数据集 > x <- HairEyeColor[,,1]+HairEyeColor[,,2] > n <- 100 > m <- matrix(sample(c(T,F),n^2,replace=T), nr=n, nc=n) > png(file="plot5.png",width=640,height=480) # PNG图 > biplot(corresp(m, nf=2), main="Plot") > dev.off() > CairoPNG(file="Cairo5.png",width=640,height=480) > biplot(corresp(m, nf=2), main="Cairo") > dev.off()
图1-14 以png()函数生成的图
图1-15 以CairoPNG()函数生成的图
在当前目录,会生成2个png文件,即plot5.png和Cairo5.png,见图1-16和图1-17。
图1-16 以png()函数生成的图
图1-17 以CairoPNG()函数生成的图
我们查看两个文件的属性,发现以png直接生成的图54KB,以CairoPNG生成的图43.8KB,如图1-18所示。
图1-18 对比两幅图的文件属性
SVG图形输出代码为:
> svg(file="plot-svg5.svg",width=6,height=6) > biplot(corresp(m, nf=2), main="Plot-svg") > dev.off() > CairoSVG(file="Cairo-svg5.svg",width=6,height=6) > biplot(corresp(m, nf=2), main="Cairo-svg") > dev.off()
在当前目录,会生成2个svg文件,即plot-svg5.svg和Cairo-svg5.svg。
就以上的3个例子来说,我们分辨不出CairoPNG() 和 png()之间有太大区别,只是Cairo感觉更淡、更柔和一些。关于这一点谢益辉补充说,看不出区别是因为现在R的png()设备本来就是默认用cairo,几年前png()不用cairo,所以有时候PNG图的质量很差。多数情况下是几乎没有区别,只有少数情况下抗锯齿的表现不一样。
1.8 caTools:一个奇特的工具集
问题
R除了统计,还能干什么?
引言
R语言生来就是自由的,不像Java和PHP等有统一的规范约束。R语言不仅命名、语法各包各异,就连功能也是各种混搭。caTools库就是这种混搭库,包括了不相关的几组函数工具集,有图片处理的,有编解码的,有分类器的,有向量计算的,有科学计算的,而且都很好用!以至于我都不知道如何用简短的语言去描述这个包了!唯有用“奇特”来概括它的特点。
1.8.1 caTools介绍
caTools是一个基础的工具包,包括移动窗口(Moving Window)统计,二进制图片读写,快速计算曲线的面积AUC,LogitBoost分类器,base64的编码器/解码器,快速计算舍入、误差、求和、累计求和等函数。下面分别介绍caTools包中的API。
(1)二进制图片读写
read.gif & write.gif: gif格式图片的读写。
read.ENVI & write.ENVI: ENVI格式图片的读写,如GIS图片。
(2)base64的编码器/解码器:
base64encode: 编码器。
base64decode: 解码器。
注 Base64是一种基于64个可打印字符来表示二进制数据的表示方法。由于2的6次方等于64,所以每6个位元为一个单元,对应某个可打印字符。
(3)快速计算曲线的面积AUC
colAUC: 计算ROC曲线的面积(AUC)。
combs: 向量元素的无序组合。
trapz: 数值积分形梯形法则。
(4)LogitBoost分类器
LogitBoost: logitBoost分类算法。
predict.LogitBoost: 预测logitBoost分类。
sample.split: 把数据切分成训练集和测试集。
(5)快速计算工具
runmad: 计算向量中位数。
runmean: 计算向量均值。
runmin & runmax: 计算向量最小值和最大值。
runquantile: 计算向量位数。
runsd: 计算向量标准差。
sumexact, cumsumexact: 无误差求和,针对编程语言关于double类型的精度优化。
1.8.2 caTools安装
本节使用的系统环境是:
Linux: Ubuntu Server 12.04.2 LTS 64bit
R: 3.0.1 x86_64-pc-linux-gnu
注 caTools同时支持Windows 7环境和Linux环境。
安装caTools包非常简单,只需要一条命令就可以完成。
~ R # 启动R程序 > install.packages("caTools") # 执行caTools安装命令
1.8.3 caTools使用
在R环境中,加载caTools类库:
> library(caTools)
1.二进制图片读写gif
(1)写一个gif图片
> write.gif(volcano, "volcano.gif", col=terrain.colors, flip=TRUE, scale="always", comment="Maunga Whau Volcano") #取datasets::volcano数据集,写入volcano.gif
(2)读一个gif图片到内存,再从内存输出
> y = read.gif("volcano.gif", verbose=TRUE, flip=TRUE) #读入图片到变量y GIF image header Global colormap with 256 colors Comment Extension Image [61 x 87]: 3585 bytes GIF Terminator > attributes(y) #查看变量的属性 $names [1] "image" "col" "transparent" "comment" > class(y$image) #查看图片存储矩阵 [1] "matrix" > ncol(y$image) #列数 [1] 61 > nrow(y$image) #行数 [1] 87 > head(y$col,10) #颜色表,取前10个 [1] "#00A600FF" "#01A600FF" "#03A700FF" "#04A700FF" "#05A800FF" "#07A800FF" "#08A900FF" [8] "#09A900FF" "#0BAA00FF" "#0CAA00FF" > y$comment #查看图片备注 [1] "Maunga Whau Volcano" > Figure-(y$image, col=y$col, main=y$comment, asp=1) #通过y变量画图,如图1-19所示(注:读者可自己运行代码,查看生成的彩色图片。——编辑注")
(3)创建一个Git动画
> x <- y <- seq(-4*pi, 4*pi, len=200) > r <- sqrt(outer(x^2, y^2, "+")) > image = array(0, c(200, 200, 10)) > for(i in 1:10) Figure-[,,i] = cos(r-(2*pi*i/10))/(r^.25) > write.gif(image, "wave.gif", col="rainbow") > y = read.gif("wave.gif") > for(i in 1:10) Figure-(y$image[,,i], col=y$col, breaks=(0:256)-00.5, asp=1)
在当前的操作目录,会生成一个 wave.gif 的文件,如图1-20所示。动画演示效果参见本书在线资源中的文件wave.gif。
图1-19 火山地形图
图1-20 同心圆变化动画
我们看到,caTools与谢益辉的animation包有一样的GIF输出动画功能。但使用caTools做gif动画,不需要依赖其他软件库。而animation的saveGIF函数需要依赖于ImageMagick或者GraphicsMagick等第三方软件。
2.Base64的编码器/解码器
Base64常用于处理文本数据的场合,表示、传输、存储一些二进制数据,包括MIME的email、email via MIME、在XML中存储复杂数据。Base64的编码解码,只支持向量(vector)类型,不支持data.frame和list类型。把一个boolean向量编解码,代码如下:
> size=1 # 设置每个元素占用的字节数 > x = (10*runif(10)>5) > y = base64encode(x, size=size) > z = base64decode(y, typeof(x), size=size) > x # 原始数据 [1] FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE > y # 编码后的密文 [1] "AAAAAAEBAAAAAA==" > z #解码后的明文 [1] FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE
把一个字符串编解码,代码如下:
> x = "Hello R!!" # character > y = base64encode(x) > z = base64decode(y, typeof(x)) > x # 原始数据 [1] "Hello R!!" > y # 编码后的密文 [1] "SGVsbG8gUiEh" > z # 解码后的明文 [1] "Hello R!!"
错误测试:把一个数据框编解码。
> data(iris) > class(iris) [1] "data.frame" > head(x) #原始数据 Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa > base64encode(x) #对数据框进行编码 Error in writeBin(x, raw(), size = size, endian = endian) : can only write vector objects
3.ROC曲线
ROC曲线就是接收者操作特征曲线(receiver operating characteristic curve)的简称,这是一种坐标图式的分析工具,主要用于选择最佳的信号侦测模型、舍弃次佳的模型和在同一模型中设定最佳阈值。ROC曲线首先由二战中的电子工程师和雷达工程师发明,用来侦测战场上的敌军载具(飞机、船舰等),即信号检测,之后很快就被引进心理学来进行信号的知觉检测。数十年来,ROC分析被用于医学、无线电、生物学、犯罪心理学等领域,而且最近在机器学习(machine learning)和数据挖掘(data mining)领域也得到了很好的发展。
取MASS::cats数据集,3列分别是Sex(性别)、Bwt(体重)、Hwt(心脏重量)。
> library(MASS) > data(cats) #加载数据集 > head(cats) #打印前6行数据 Sex Bwt Hwt 1 F 2.0 7.0 2 F 2.0 7.4 3 F 2.0 9.5 4 F 2.1 7.2 5 F 2.1 7.3 6 F 2.1 7.6 > colAUC(cats[,2:3], cats[,1], plotROC=TRUE) #计算ROC的曲线面积AUC,输出图片,如图1-21所示 Bwt Hwt F vs.M 0.8338451 0.759048
从AUC判断分类器(预测模型)优劣的标准:
AUC = 1,是完美分类器,采用这个预测模型时,不管设定什么阈值都能得出完美预测。在绝大多数预测的场合,不存在完美分类器。
0.5 < AUC < 1,优于随机猜测。这个分类器(模型)如果妥善设定阈值的话,能有预测价值。
AUC = 0.5,跟随机猜测一样(例如丢硬币),模型没有预测价值。
AUC < 0.5,比随机猜测还差;但只要总是反预测而行,就优于随机猜测。
从图1-21我们看到Bwt和Hwt都在(0.5,1)之间,因此,cats的数据集是一个真实有效数据集。如果cats的数据集,是一个通过分类预测的数据集,用AUC对数据集的评分,就可以检验分类器的好坏了。
图1-21 ROC曲线
4.向量元素的无序组合
combs(v,k)函数,用于创建无序的组合的矩阵,矩阵的列表示以几种元素进行组合,矩阵的行表示每种不同的组合。参数v是向量;k是数值,小于等于v的长度[1:length(v)]。
> combs(2:5, 3) [,1] [,2] [,3] [1,] 2 3 4 [2,] 2 3 5 [3,] 2 4 5 [4,] 3 4 5 > combs(c("cats", "dogs", "mice"), 2) [,1] [,2] [1,] "cats" "dogs" [2,] "cats" "mice" [3,] "dogs" "mice" > a = combs(1:4, 2) #快速组合构建矩阵 > a [,1] [,2] [1,] 1 2 [2,] 1 3 [3,] 1 4 [4,] 2 3 [5,] 2 4 [6,] 3 4 > b = matrix( c(1,1,1,2,2,3,2,3,4,3,4,4), 6, 2) > b [,1] [,2] [1,] 1 2 [2,] 1 3 [3,] 1 4 [4,] 2 3 [5,] 2 4 [6,] 3 4
5.数值积分梯形法则
梯形法则原理:将被积函数近似为直线函数,被积的部分近似为梯形。
> x = (1:10)*pi/10 > trapz(x, sin(x)) [1] 1.934983 > x = (1:1000)*pi/1000 > trapz(x, sin(x)) [1] 1.999993
6.LogitBoost分类器
取datasets::iris数据集,5列分别是Sepal.Length(花萼长)、Sepal.Width((花萼宽)、Petal.Length(花瓣长)、Petal.Width(花瓣宽)、Species(种属)。
> head(iris) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 4.9 3.0 1.4 0.2 setosa 3 4.7 3.2 1.3 0.2 setosa 4 4.6 3.1 1.5 0.2 setosa 5 5.0 3.6 1.4 0.2 setosa 6 5.4 3.9 1.7 0.4 setosa > Data = iris[,-5] > Label = iris[, 5] > model = LogitBoost(Data, Label, nIter=20) #训练模型 > model # 模型数据 $Stump feature threshhold sign feature threshhold sign feature threshhold sign [1,] 3 1.9 -1 2 2.9 -1 4 1.6 1 [2,] 4 0.6 -1 3 4.7 -1 3 4.8 1 [3,] 3 1.9 -1 2 2.0 -1 4 1.7 1 [4,] 4 0.6 -1 3 1.9 1 3 4.9 1 [5,] 3 1.9 -1 4 1.6 -1 4 1.3 1 [6,] 4 0.6 -1 1 6.5 1 2 2.6 -1 [7,] 3 1.9 -1 3 1.9 1 4 1.7 1 [8,] 4 0.6 -1 2 2.0 -1 2 3.0 -1 [9,] 3 1.9 -1 3 5.0 -1 3 5.0 1 [10,] 4 0.6 -1 2 2.9 1 1 4.9 -1 [11,] 3 1.9 -1 3 1.9 1 3 4.4 1 [12,] 4 0.6 -1 2 2.0 -1 4 1.7 1 [13,] 3 1.9 -1 3 5.1 -1 2 3.1 -1 [14,] 4 0.6 -1 2 2.0 -1 3 5.1 1 [15,] 3 1.9 -1 3 1.9 1 1 6.5 -1 [16,] 4 0.6 -1 4 1.6 -1 3 5.1 1 [17,] 3 1.9 -1 2 3.1 1 2 3.1 -1 [18,] 4 0.6 -1 3 1.9 1 1 4.9 -1 [19,] 3 1.9 -1 2 2.0 -1 4 1.4 1 [20,] 4 0.6 -1 3 5.1 -1 2 2.2 -1 $lablist [1] setosa versicolor virginica Levels: setosa versicolor virginica attr(,"class") [1] "LogitBoost" > Lab = predict(model, Data) # 分类预测,Lab只显示分类的结果, Prob显示各分类的概率 > Prob = predict(model, Data, type="raw") > t = cbind(Lab, Prob) # 结果合并,打印前6列 > head(t) Lab setosa versicolor virginica [1,] 1 1 0.017986210 1.522998e-08 [2,] 1 1 0.002472623 3.353501e-04 [3,] 1 1 0.017986210 8.315280e-07 [4,] 1 1 0.002472623 4.539787e-05 [5,] 1 1 0.017986210 1.522998e-08 [6,] 1 1 0.017986210 1.522998e-08
前6条数据,Lab列表示数据属于分类1,即setosa。其他3列setosa、versicolor、virginica,分类表示属于该分类的概率是多少。下面设置迭代次数,比较分类结果与实际数据结果。
> table(predict(model, Data, nIter= 2), Label) #设置迭代次数为2 Label setosa versicolor virginica setosa 48 0 0 versicolor 0 45 1 virginica 0 3 45 > table(predict(model, Data, nIter=10), Label) #设置迭代次数为10 Label setosa versicolor virginica setosa 50 0 0 versicolor 0 47 0 virginica 0 1 47 > table(predict(model, Data), Label) #默认迭代次数, 训练时LogitBoost的nIter值 Label setosa versicolor virginica setosa 50 0 0 versicolor 0 49 0 virginica 0 0 48
从上面3次测试结果可以看出,迭代次数越多模型分类越准确。下面随机划分训练集和测试集,并分类预测。
> mask = sample.split(Label) # 随机取训练集 > length(which(mask)) # 训练集99条记录 [1] 99 > length(which(!mask)) # 测试集51条记录 [1] 51 > model = LogitBoost(Data[mask,], Label[mask], nIter=10) #训练模型 > table(predict(model, Data[!mask,], nIter=2), Label[!mask]) #分类预测 setosa versicolor virginica setosa 16 0 0 versicolor 0 15 3 virginica 0 1 12 > table(predict(model, Data[!mask,]), Label[!mask]) setosa versicolor virginica setosa 17 0 0 versicolor 0 16 4 virginica 0 1 13
7.快速计算工具runmean
均线在股票交易上非常流行,是一种简单、实用的看盘指标。下面我们对时间序列数据取均线,输出结果如图1-22所示。
> BJsales # 取datasets::BJsales数据集 Time Series: Start = 1 End = 150 Frequency = 1 [1] 200.1 199.5 199.4 198.9 199.0 200.2 198.6 200.0 200.3 201.2 201.6 201.5 [13] 201.5 203.5 204.9 207.1 210.5 210.5 209.8 208.8 209.5 213.2 213.7 215.1 [25] 218.7 219.8 220.5 223.8 222.8 223.8 221.7 222.3 220.8 219.4 220.1 220.6 > plot(BJsales,col="black", lty=1,lwd=1, main = "Moving Window Means") > lines(runmean(BJsales, 3), col="red", lty=2, lwd=2) > lines(runmean(BJsales, 8), col="green", lty=3, lwd=2) > lines(runmean(BJsales,15), col="blue", lty=4, lwd=2) > lines(runmean(BJsales,24), col="magenta", lty=5, lwd=2) > lines(runmean(BJsales,50), col="cyan", lty=6, lwd=2)
图1-22 时间序列数据的均线图
我们从图1-22中看到6条线, 黑色是原始数据,其他各种颜色线代表不同单位的均线。比如,红色(red)线表示3个点的平均,绿色(green)线表示8个点的平均。
8.快速计算工具组合
对数据集取最大路径(runmax)、最小路径(runmin)、均值路径(runmean)和中位数路径(runmed),输出结果如图1-23所示。
> n=200;k=25 > set.seed(100) > x = rnorm(n,sd=30) + abs(seq(n)-n/4) > plot(x, main = "Moving Window Analysis Functions (window size=25)") > lines(runmin (x,k), col="red",lty=1, lwd=1) > lines(runmax (x,k), col="cyan",lty=1, lwd=1) > lines(runmean(x,k), col="blue",lty=1, lwd=1) > lines(runmed (x,k), col="green",lty=2, lwd=2)
图1-23 路径图
9.无误差求和
各种编程语言,用计算机做小数计算的时候,都会出现计算误差的情况,R语言也有这个问题。首先用sum()函数求和。
> x = c(1, 1e20, 1e40, -1e40, -1e20, -1) #计算误差求和 > a = sum(x); print(a) [1] -1e+20 > b = sumexact(x); print(b) # 无计算误差求和 [1] 0
我们对向量x求和,各组值加起来正好为0。但是sum()函数,结果是-1e+20,这是由于编程精度的问题造成的计算误差。通过sumexact()函数修正后,就没有误差了。下面用cumsum()函数累积求和。
> a = cumsum(x); print(a) # 计算误差累积求和 [1] 1e+00 1e+20 1e+40 0e+00 -1e+20 -1e+20 > b = cumsumexact(x); print(b) # 无计算误差累积求和 [1] 1e+00 1e+20 1e+40 1e+20 1e+00 0e+00
cumsum()函数同样出现了精度的误差,需要用cumsumexact()函数来修正。最后还是要用“奇特”来概括这个工具集,相信你也发现了它的奇特。