下载掌阅APP，畅读海量书库
立即打开
本章主要介绍了时间序列数据处理的3个工具包，帮助读者掌握时间序列在R语言中的数据结构和基本使用。
问题
R语言怎么处理时间序列数据？
引言
时间序列分析是一种动态数据处理的统计方法，通过对时间序列数据的分析，我们可以感觉到世界正改变着什么！R语言作为统计分析的利器，对时间序列处理有着强大的支持。在R语言中，单独为时间序列数据定义了一种数据类型zoo，zoo是时间序列的基础，也是股票分析的基础。本节将介绍zoo库在R语言中的结构和使用。
2.1.1 zoo包介绍
zoo是一个R语言类库，zoo类库中定义了一个名为zoo的S3类型对象，用于描述规则的和不规则的有序的时间序列数据。zoo对象是一个独立的对象，包括索引、日期、时间，只依赖于基础的R环境。zooreg对象继承了zoo对象，只能用于规则的时间序列数据。R语言中很多其他的程序包，都是以zoo和zooreg作为时间序列数据的基础的！zoo包的API主要有6类，下面一一介绍。
（1）基础对象
·zoo:有序的时间序列对象。
·zooreg:规则的时间序列对象，继承zoo对象。与zoo相比，不同之处在于zooreg要求数据是连续的。
（2）类型转换
·as.zoo:把一个对象转型为zoo类型。
·plot.zoo:为plot函数提供zoo的接口。
·xyplot.zoo:为lattice的xyplot函数提供zoo的接口。
·ggplot2.zoo:为ggplot2包提供zoo的接口。
（3）数据操作
·coredata:查看或编辑zoo的数据部分。
·index:查看或编辑zoo的索引部分。
·window.zoo:按时间过滤数据。
·merge.zoo:合并多个zoo对象。
·read.zoo:从文件读写zoo序列。
·aggregate.zoo:计算zoo数据。
·rollapply:对zoo数据的滚动处理。
·rollmean:对zoo数据的滚动计算均值。
（4）NA值处理
·na.fill:NA值的填充。
·na.locf:替换NA值。
·na.aggregate:计算统计值替换NA值。
·na.approx:计算插值替换NA值。
·na.StructTS:计算季节Kalman滤波替换NA值。
·na.trim:过滤有NA的记录。
（5）辅助工具
·is.regular:检查是否是规则的序列。
·lag.zoo:计算步长和差分。
·MATCH:取交集。
·ORDER:值排序，输出索引。
（6）显示控制
·yearqtr:以年季度显示时间。
·yearmon:以年月显示时间。
·xblocks:作图沿x轴分割图形。
·make.par.list:用于给plot.zoo和xyplot.zoo数据格式转换。
2.1.2 zoo安装
本节使用的系统环境是：
·Win764bit
·R:3.0.1x86_64-w64-mingw32/x64b4bit
注 zoo同时支持Windows 7环境和Linux环境。
zoo包的安装过程如下：
~ R #
启动R
程序
> install.packages("zoo") #
安装zoo
包
> library(zoo) #
加载zoo
包
2.1.3 zoo包的使用
1.zoo对象
zoo对象包括两部分，即数据部分和索引部分。首先是函数定义：
zoo(x = NULL, order.by = index(x), frequency = NULL)
其中x是数据部分，允许类型为向量、矩阵、因子；order.by是索引部分，字段唯一性要求，用于排序；frequency是每个时间单元显示的数量。
以下代码构建一个zoo对象，以时间为索引，产生的结果是图2-1。特别要注意的一点是，zoo对象可以接受不连续的时间序列数据。
> x.Date <- as.Date("2003-02-01") + c(1, 3, 7, 9, 14)
– 1 #
定义一个不连续的日期的向量
> x.Date
[1] "2003-02-01" "2003-02-03" "2003-02-07" "2003-02-09" "2003-02-14"
> class(x.Date)
[1] "Date"
> x <- zoo(rnorm(5), x.Date) #
定义不连续的zoo
对象
> x
2003-02-01 2003-02-03 2003-02-07 2003-02-09 2003-02-14
0.01964254 0.03122887 0.64721059 1.47397924 1.29109889
> class(x)
[1] "zoo"
> plot(x) #
画图显示
图2-1 时间序列图
接下来，我们以数字为索引创建多组时间序列。用如下代码，生成一个有12个元素的4行3列的矩阵，以数字0:10为索引，创建一个zoo类型对象y，并以图形输出y，产生的结果如图2-2所示。
> y <- zoo(matrix(1:12, 4, 3),0:10) > y 0 1 5 9 1 2 6 10 2 3 7 11 3 4 8 12 4 1 5 9 5 2 6 10 6 3 7 11 7 4 8 12 8 1 5 9 9 2 6 10 10 3 7 11 > plot(y) # 矩阵的每一列为一组时间序列图
图2-2 分组的时间序列图
2.zooreg对象
首先是函数定义：
zooreg(data, start = 1, end = numeric(), frequency = 1,
deltat = 1, ts.eps = getOption("ts.eps"), order.by = NULL)
下面是参数说明。
·data:数据部分，允许类型为向量、矩阵、因子。
·start:时间部分，开始时间。
·end:时间部分，结束时间。
·frequency:每个时间单元显示的数量。
·deltat:连续观测的采样周期，不能与frequency同时出现，例如，取每月的数据，为1/12。
·ts.eps:时间序列间隔，当数据时间间隔小于ts.eps时，使用ts.eps作为时间间隔。通过getOption(“ts.eps”)设置，默认是1e-05。
·order.by:索引部分，字段唯一性要求，用于排序,继承zoo的order.by。
以下代码构建一个zooreg对象，以连续的年（季度）时间为索引，产生的结果是图2-3。
> zooreg(1:10, frequency = 4, start = c(1959, 2))
1959(2) 1959(3) 1959(4) 1960(1) 1960(2) 1960(3) 1960(4) 1961(1) 1961(2) 1961(3)
1 2 3 4 5 6 7 8 9 10
> as.zoo(ts(1:10, frequency = 4, start = c(1959, 2)))
1959(2) 1959(3) 1959(4) 1960(1) 1960(2) 1960(3) 1960(4) 1961(1) 1961(2) 1961(3)
1 2 3 4 5 6 7 8 9 10
> zr<-zooreg(rnorm(10), frequency = 4, start = c(1959, 2))
> plot(zr)
图2-3 zooreg对象的时间序列图
3.zoo对象与zooreg对象的区别
zoo对象与zooreg对象的区别体现在计算步长和差分方面。
·lag（步长）：zoo根据索引计算，zooreg根据值计算。
·diff（差分）：zoo根据索引计算，zooreg根据值计算。
例如，对同一组值不连续的数据（1,2,3,6,7,8），二者的计算结果如下。
> x <- c(1, 2, 3, 6, 7, 8) > zz <- zoo(x, x) > zr <- as.zooreg(zz) > lag(zz, k = -1) # 计算步长 2 3 6 7 8 1 2 3 6 7 > lag(zr, k = -1) 2 3 4 7 8 9 1 2 3 6 7 8 > diff(zz) # 计算差分 2 3 6 7 8 1 1 3 1 1 > diff(zr) 2 3 7 8 1 1 1 1
4.zoo对象的类型转换
首先，把对象从其他类型转型到zoo类型。
> as.zoo(rnorm(5)) #
把一个基本类型的向量转型到zoo
类型
1 2 3 4 5
-0.4892119 0.5740950 0.7128003 0.6282868 1.0289573
> as.zoo(ts(rnorm(5), start = 1981, freq = 12))
1981(1) 1981(2) 1981(3) 1981(4) 1981(5)
2.3198504 0.5934895 -1.9375893 -1.9888237 1.0944444
> x <- as.zoo(ts(rnorm(5), start = 1981, freq = 12)); x #
把一个ts
类型转型到zoo
类型
1981(1) 1981(2) 1981(3) 1981(4) 1981(5)
1.8822996 1.6436364 0.1260436 -2.0360960 -0.1387474
其次，把对象从zoo类型转型到其他类型。
> as.matrix(x) #
把zoo
类型，转型到矩阵
x
1981(1) 1.8822996
1981(2) 1.6436364
1981(3) 0.1260436
1981(4) -2.0360960
1981(5) -0.1387474
> as.vector(x) #
把zoo
类型，转型到数字向量
[1] 1.8822996 1.6436364 0.1260436 -2.0360960 -0.1387474
> as.data.frame(x) #
把zoo
类型，转型到数据框
x
1981(1) 1.8822996
1981(2) 1.6436364
1981(3) 0.1260436
1981(4) -2.0360960
1981(5) -0.1387474
> as.list(x) #
把zoo
类型，转型到列表
[[1]]
1981(1) 1981(2) 1981(3) 1981(4) 1981(5)
1.8822996 1.6436364 0.1260436 -2.0360960 -0.1387474
5.用ggplot2画时间序列
由于ggplot2不支持zoo类型的数据，因此需要通过ggplot2::fortify()函数，调用zoo::fortify.zoo()函数，把zoo类型转换成ggplot2可识别的类型后，ggplot2才可以对zoo类型数据的画图。以下代码用ggplot2画zoo类型的时间序列图，产生的结果是图2-4。
> library(scales)
> x.Date <- as.Date(paste(2003, 02, c(1, 3, 7, 9, 14), sep = "-")) #
构建数据对象
> x <- zoo(rnorm(5), x.Date)
> xlow <- x - runif(5)
> xhigh <- x + runif(5)
> z <- cbind(x, xlow, xhigh)
> z #
显示数据集
x xlow xhigh
2003-02-01 -0.36006612 -0.88751958 0.006247816
2003-02-03 1.35216617 0.97892538 2.076360524
2003-02-07 0.61920828 0.23746410 1.156569424
2003-02-09 0.27516116 0.09978789 0.777878867
2003-02-14 0.02510778 -0.80107410 0.541592929
#
对zoo
类型的数据，用fortify()
转换成data.frame
类型
> g<-ggplot(aes(x = Index, y = Value), data = fortify(x, melt = TRUE))
> g<-g+geom_line()
> g<-g+geom_line(aes(x = Index, y = xlow), colour = "red", data = fortify(xlow))
> g<-g+geom_ribbon(aes(x = Index, y = x, ymin = xlow, ymax = xhigh), data =fortify(x), fill = "darkgray")
> g<-g+geom_line()
> g<-g+xlab("Index") + ylab("x")
> g
图2-4 用ggplot2画的时间序列图
6.zoo对象的数据操作
使用coredata()函数修改zoo类型的数据部分。
> x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
> x <- zoo(matrix(rnorm(20), ncol = 2), x.date)
> coredata(x) #
查看数据部分
[,1] [,2]
[1,] -1.04571765 0.92606273
[2,] -0.89621126 0.03693769
[3,] 1.26938716 -1.06620017
[4,] 0.59384095 -0.23845635
[5,] 0.77563432 1.49522344
[6,] 1.55737038 1.17215855
[7,] -0.36540180 -1.45770721
[8,] 0.81655645 0.09505623
[9,] -0.06063478 0.84766496
[10,] -0.50137832 -1.62436453
> coredata(x) <- matrix(1:20, ncol = 2) #
修改数据部分
> x #
查看修改后的数据集
2003-01-01 1 11
2003-01-03 2 12
2003-01-05 3 13
2003-01-07 4 14
2003-02-09 5 15
2003-02-11 6 16
2003-02-13 7 17
2003-03-15 8 18
2003-03-17 9 19
2003-04-19 10 20
使用index()函数修改zoo类型的索引部分。
> x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-")) > x <- zoo(matrix(rnorm(20), ncol = 2), x.date) > index(x) # 查看索引部分 [1] "2003-01-01" "2003-01-03" "2003-01-05" "2003-01-07" "2003-02-09" [6] "2003-02-11" "2003-02-13" "2003-03-15" "2003-03-17" "2003-04-19" > index(x) <- 1:nrow(x) # 修改索引部分 > index(x) # 查看修改后的索引部分 [1] 1 2 3 4 5 6 7 8 9 10
使用window.zoo()函数按时间过滤数据。
> x.date <- as.Date(paste(2003, rep(1:4, 4:1), seq(1,20,2), sep = "-"))
> x <- zoo(matrix(rnorm(20), ncol = 2), x.date)
> window(x, start = as.Date("2003-02-01"), end = as.Date("2003-03-01"))
#
取日期从2003-02-01
到2003-03-01
之间的数据
2003-02-09 0.7021167 -0.3073809
2003-02-11 2.5071111 0.6210542
2003-02-13 -1.8900271 0.1819022
> window(x, index = x.date[1:6], start = as.Date("2003-02-01"))
#
取日期从2003-02-01
开始的，且索引日期在x.date[1:6]
中的数据
2003-02-09 0.7021167 -0.3073809
2003-02-11 2.5071111 0.6210542
> window(x, index = x.date[c(4, 8, 10)]) #
取索引日期在x.date[c(4, 8, 10)]
中的数据
2003-01-07 1.4623515 -1.198597
2003-03-15 -0.5898128 1.318401
2003-04-19 -0.4209979 -1.648222
使用merge.zoo()合并多个zoo对象。
> y1 <- zoo(matrix(1:10, ncol = 2), 1:5);y1 #
创建2
个zoo
数据
1 1 6
2 2 7
3 3 8
4 4 9
5 5 10
> y2 <- zoo(matrix(rnorm(10), ncol = 2), 3:7);y2
3 1.4810127 0.13575871
4 -0.3914258 0.06404148
5 0.6018237 1.85017952
6 1.2964150 -0.12927481
7 0.2211769 0.32381709
> merge(y1, y2, all = FALSE) #
以相同的索引值合并数据
y1.1 y1.2 y2.1 y2.2
3 3 8 0.9514985 1.7238941
4 4 9 -1.1131230 -0.2061446
5 5 10 0.6169665 -1.3141951
> merge(y1, y2, all = FALSE, suffixes = c("a", "b")) #
自定义数据列的名字
a.1 a.2 b.1 b.2
3 3 8 0.9514985 1.7238941
4 4 9 -1.1131230 -0.2061446
5 5 10 0.6169665 -1.3141951
> merge(y1, y2, all = TRUE) #
合并完整的数据集，空数据默认以NA
填充
y1.1 y1.2 y2.1 y2.2
1 1 6 NA NA
2 2 7 NA NA
3 3 8 0.9514985 1.7238941
4 4 9 -1.1131230 -0.2061446
5 5 10 0.6169665 -1.3141951
6 NA NA 0.5134937 0.0634741
7 NA NA 0.3694591 -0.2319775
> merge(y1, y2, all = TRUE, fill = 0) #
合并完整的数据集，空数据以0
填充
y1.1 y1.2 y2.1 y2.2
1 1 6 0.0000000 0.0000000
2 2 7 0.0000000 0.0000000
3 3 8 0.9514985 1.7238941
4 4 9 -1.1131230 -0.2061446
5 5 10 0.6169665 -1.3141951
6 0 0 0.5134937 0.0634741
7 0 0 0.3694591 -0.2319775
使用aggregate.zoo()函数对zoo数据进行计算。
> x.date <- as.Date(paste(2004, rep(1:4, 4:1), seq(1,20,2), sep = "-")) # 创建zoo 类型数据集x > x <- zoo(rnorm(12), x.date); x 2004-01-01 2004-01-03 2004-01-05 2004-01-07 2004-02-09 2004-02-11 0.67392868 1.95642526 -0.26904101 -1.24455152 -0.39570292 0.09739665 2004-02-13 2004-03-15 2004-03-17 2004-04-19 -0.23838695 -0.41182796 -1.57721805 -0.79727610 > x.date2 <- as.Date(paste(2004, rep(1:4, 4:1), 1, sep = "-")); x.date2 # 创建时间向量x.date2 [1] "2004-01-01" "2004-01-01" "2004-01-01" "2004-01-01" "2004-02-01" [6] "2004-02-01" "2004-02-01" "2004-03-01" "2004-03-01" "2004-04-01" > x2 <- aggregate(x, x.date2, mean); x2 # 计算x 以x.date2 为时间分割规则的均值 2004-01-01 2004-02-01 2004-03-01 2004-04-01 0.2791904 -0.1788977 -0.9945230 -0.7972761
7.zoo对象数据函数化处理
使用rollapply()函数对zoo数据进行函数化处理。
> z <- zoo(11:15, as.Date(31:35))
> rollapply(z, 2, mean) #
从起始日开始，计算连续2
日的均值
1970-02-01 1970-02-02 1970-02-03 1970-02-04
11.5 12.5 13.5 14.5
> rollapply(z, 3, mean) #
从起始日开始，计算连续3
日的均值
1970-02-02 1970-02-03 1970-02-04
12 13 14
等价操作变换：用rollapply()实现aggregate()的操作。
> z2 <- zoo(rnorm(6))
> rollapply(z2, 3, mean, by = 3) # means of nonoverlapping groups of 3
2 5
-0.3065197 0.6350963
> aggregate(z2, c(3,3,3,6,6,6), mean) # same
3 6
-0.3065197 0.6350963
等价操作变换：用rollapply()实现rollmean()的操作。
> rollapply(z2, 3, mean) # uses rollmean which is optimized for mean
2 3 4 5
-0.3065197 -0.7035811 -0.1672344 0.6350963
> rollmean(z2, 3) # same
2 3 4 5
-0.3065197 -0.7035811 -0.1672344 0.6350963
8.NA值处理
使用na.fill()函数进行NA填充。
> z <- zoo(c(NA, 2, NA, 3, 4, 5, 9, NA));z #
创建有NA
值的zoo
对象
1 2 3 4 5 6 7 8
NA 2 NA 3 4 5 9 NA
> na.fill(z, "extend") #
用extend
的方法，填充NA
值，即NA
前后项的均值填充
1 2 3 4 5 6 7 8
2.0 2.0 2.5 3.0 4.0 5.0 9.0 9.0
> na.fill(z, -(1:3)) #
自定义填充NA
值，即-(1:3)
循环填充
1 2 3 4 5 6 7 8
-1 2 -2 3 4 5 9 -3
> na.fill(z, c("extend", NA)) #
用extend
，配合自定义的方法，即extend
和自定义规则循环填充
1 2 3 4 5 6 7 8
2 2 NA 3 4 5 9 9
使用na.locf()函数进行NA替换。
> z <- zoo(c(NA, 2, NA, 3, 4, 5, 9, NA, 11));z 1 2 3 4 5 6 7 8 9 NA 2 NA 3 4 5 9 NA 11 > na.locf(z) # 用NA 的前一项的值，替换NA 值 2 3 4 5 6 7 8 9 2 2 3 4 5 9 9 11 > na.locf(z, fromLast = TRUE) # 用NA 的后一项的值，替换NA 值 1 2 3 4 5 6 7 8 9 2 2 3 3 4 5 9 11 11
使用na.aggregate()函数的统计计算的值替换NA值。
> z <- zoo(c(1, NA, 3:9),c(as.Date("2010-01-01") + 0:2,as.Date("2010-02-01") +
0:2,as.Date("2011-01-01") + 0:2));z
2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01
1 NA 3 4 5 6 7
2011-01-02 2011-01-03
8 9
> na.aggregate(z) #
计算排除NA
的其他项的均值，替换NA
值
2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01
1.000 5.375 3.000 4.000 5.000 6.000 7.000
2011-01-02 2011-01-03
8.000 9.000
> na.aggregate(z, as.yearmon) #
以索引的年月分组的均值，替换NA
值
2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01
1 2 3 4 5 6 7
2011-01-02 2011-01-03
8 9
> na.aggregate(z, months) #
以索引的月份分组的均值，替换NA
值
2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01
1.0 5.6 3.0 4.0 5.0 6.0 7.0
2011-01-02 2011-01-03
8.0 9.0
> na.aggregate(z, format, "%Y") #
以正则表示的索引的年份分组的均值，替换NA
值
2010-01-01 2010-01-02 2010-01-03 2010-02-01 2010-02-02 2010-02-03 2011-01-01
1.0 3.8 3.0 4.0 5.0 6.0 7.0
2011-01-02 2011-01-03
8.0 9.0
使用na.approx()函数计算插值替换NA值。
> z <- zoo(c(2, NA, 1, 4, 5, 2), c(1, 3, 4, 6, 7, 8));z
1 3 4 6 7 8
2 NA 1 4 5 2
> na.approx(z)
1 3 4 6 7 8
2.000000 1.333333 1.000000 4.000000 5.000000 2.000000
> na.approx(z, 1:6)
1 3 4 6 7 8
2.0 1.5 1.0 4.0 5.0 2.0
使用na.StructTS()函数计算季节Kalman滤波替换NA值，产生的结果是图2-5。
> z <- zooreg(rep(10 * seq(4), each = 4) + rep(c(3, 1, 2, 4), times = 4),
start = as.yearqtr(2000), freq = 4)
> z[10] <- NA
> zout <- na.StructTS(z);zout
> plot(cbind(z, zout), screen = 1, col = 1:2, type = c("l", "p"), pch = 20)
图2-5 替换NA值的时间序列
使用na.trim()函数，去掉有NA的行。
> xx <- zoo(matrix(c(1, 4, 6, NA, NA, 7), 3), c(2, 4, 6));xx 2 1 NA 4 4 NA 6 6 7 > na.trim(xx) 6 6 7
9.数据显示格式
以“年+季度”格式输出
> x <- as.yearqtr(2000 + seq(0, 7)/4);x #
以年季默认格式输出
[1] "2000 Q1" "2000 Q2" "2000 Q3" "2000 Q4" "2001 Q1" "2001 Q2" "2001 Q3"
[8] "2001 Q4"
> format(x, "%Y Quarter %q") #
以年季自定义格式输出
[1] "2000 Quarter 1" "2000 Quarter 2" "2000 Quarter 3" "2000 Quarter 4"
[5] "2001 Quarter 1" "2001 Quarter 2" "2001 Quarter 3" "2001 Quarter 4"
> as.yearqtr("2001 Q2")
[1] "2001 Q2"
> as.yearqtr("2001 q2")
[1] "2001 Q2"
> as.yearqtr("2001-2")
[1] "2001 Q2"
以“年+月份”格式输出
> x <- as.yearmon(2000 + seq(0, 23)/12) ;x #
以年月默认格式输出
[1] "
一月 2000" "
二月 2000" "
三月 2000" "
四月 2000" "
五月 2000"
[6] "
六月 2000" "
七月 2000" "
八月 2000" "
九月 2000" "
十月 2000"
[11] "
十一月 2000" "
十二月 2000" "
一月 2001" "
二月 2001" "
三月 2001"
[16] "
四月 2001" "
五月 2001" "
六月 2001" "
七月 2001" "
八月 2001"
[21] "
九月 2001" "
十月 2001" "
十一月 2001" "
十二月 2001"
> as.yearmon("mar07", "%b%y")
[1] NA
> as.yearmon("2007-03-01")
[1] "
三月 2007"
> as.yearmon("2007-12")
[1] "
十二月 2007"
10.区间分割
使用xblock()函数，以不同的颜色划分3个区间，即（-Inf，15）、[15，30]和（30，Inf），产生的是图2-6。
图2-6 有分割线的时间序列
> set.seed(0) > flow <- ts(filter(rlnorm(200, mean = 1), 0.8, method = "r")) > rgb <- hcl(c(0, 0, 260), c = c(100, 0, 100), l = c(50, 90, 50), alpha = 0.3) > plot(flow) > xblocks(flow > 30, col = rgb[1]) ## high values red > xblocks(flow < 15, col = rgb[3]) ## low value blue > xblocks(flow >= 15 & flow <= 30, col = rgb[2]) ## the rest gray
11.从文件读入时间序列数据创建zoo对象
我们首先创建一个文件，并将其命名为read.csv，代码如下。
~ vi read.csv 2003-01-01,1.0073644,0.05579711 2003-01-03,-0.2731580,0.06797239 2003-01-05,-1.3096795,-0.20196174 2003-01-07,0.2225738,-1.15801525 2003-02-09,1.1134332,-0.59274327 2003-02-11,0.8373944,0.76606538 2003-02-13,0.3145168,0.03892812 2003-03-15,0.2222181,0.01464681 2003-03-17,-0.8436154,-0.18631697 2003-04-19,0.4438053,1.40059083
然后读入文件并生成zoo序列。
> r <- read.zoo(file="read.csv",sep = ",", format = "%Y-%m-%d") #
以zoo
格式读入数据
> r #
查看数据
V2 V3
2003-01-01 1.0073644 0.05579711
2003-01-03 -0.2731580 0.06797239
2003-01-05 -1.3096795 -0.20196174
2003-01-07 0.2225738 -1.15801525
2003-02-09 1.1134332 -0.59274327
2003-02-11 0.8373944 0.76606538
2003-02-13 0.3145168 0.03892812
2003-03-15 0.2222181 0.01464681
2003-03-17 -0.8436154 -0.18631697
2003-04-19 0.4438053 1.40059083
> class(r) #
查看数据类型
[1] "zoo"
我们已经完全掌握了zoo库及zoo对象的使用，接下来就可以放手去用R处理时间序列数据了！
问题
如何进行复杂的时间序列数据处理？
引言
本节将继续2.1节，介绍zoo包的扩展实现。看上去简单的时间序列，却内含复杂的规律。zoo作为时间序列的基础库，是面向通用的设计，可以用来定义股票数据，也可以分析天气数据。但由于业务行为的不同，我们需要更多的辅助函数帮助我们更高效地完成任务。xts扩展了zoo，提供更多的数据处理和数据变换的函数。
2.2.1 xts介绍
xts是对时间序列数据（zoo）的一种扩展实现，目标是为了统一时间序列的操作接口。实际上，xts类型继承了zoo类型，丰富了时间序列数据处理的函数，API定义更贴近使用者，更实用，更简单！
1.xts数据结构
xts扩展zoo的基础结构，由3部分组成，如图2-7所示。
·索引部分：时间类型向量。
·数据部分：以矩阵为基础类型，支持可以与矩阵相互转换的任何类型。
·属性部分：附件信息，包括时区和索引时间类型的格式等。
图2-7 xts数据结构
2.xts的API介绍
（1）xts基础
·xts:定义xts数据类型，继承zoo类型。
·coredata.xts:查看或编辑xts对象的数据部分。
·xtsAttributes:查看或编辑xts对象的属性部分。
·[.xts]:用[]语法，取数据子集。
·dimnames.xts:查看或编辑xts维度名。
·sample_matrix:测试数据集，包括180条xts对象的记录，matrix类型。
·xtsAPI:C语言API接口。
（2）类型转换
·as.xts:转换对象到xts（zoo）类型。
·as.xts.methods:转换对象到xts函数。
·plot.xts:为plot函数提供xts的接口作图。
·.parseISO8601:把字符串（ISO8601格式）输出为POSIXct类型的，包括开始时间和结束时间的list对象。
·firstof:创建一个开始时间，POSIXct类型。
·lastof:创建一个结束时间，POSIXct类型。
·indexClass:取索引类型。
·.indexDate:索引的日期。
·.indexday:索引的日期，同.indexDate。
·.indexyday:索引的年（日）值。
·.indexmday:索引的月（日）值。
·.indexwday:索引的周（日）值。
·.indexweek:索引的周值。
·.indexmon:索引的月值。
·.indexyear:索引的年值。
·.indexhour:索引的时值。
·.indexmin:索引的分值。
·.indexsec:索引的秒值。
（3）数据处理
·align.time:以下一个时间对齐数据，秒，分钟，小时。
·endpoints:按时间单元提取索引数据。
·merge.xts:合并多个xts对象，重写zoo::merge.zoo函数。
·rbind.xts:数据按行合并，为rbind函数提供xts的接口。
·split.xts:数据分割，为split函数，提供xts的接口。
·na.locf.xts:替换NA值，重写zoo:na.locf函数。
（4）数据统计
·apply.daily:按日分割数据，执行函数。
·apply.weekly:按周分割数据，执行函数。
·apply.monthly:按月分割数据，执行函数。
·apply.quarterly:按季分割数据，执行函数。
·apply.yearly:按年分割数据，执行函数。
·to.period:按期间分割数据。
·period.apply:按期间执行自定义函数。
·period.max:按期间计算最大值。
·period.min:按期间计算最小值。
·period.prod:按期间计算指数。
·period.sum:按期间求和。
·nseconds:计算数据集包括多少秒。
·nminutes:计算数据集包括多少分。
·nhours:计算数据集包括多少时。
·ndays:计算数据集包括多少日。
·nweeks:计算数据集包括多少周。
·nmonths:计算数据集包括多少月。
·nquarters:计算数据集包括多少季。
·nyears:计算数据集包括多少年。
·periodicity:查看时间序列的期间。
（5）辅助工具
·first:从开始到结束设置条件取子集。
·last:从结束到开始设置条件取子集。
·timeBased:判断是否是时间类型。
·timeBasedSeq:创建时间的序列。
·diff.xts:计算步长和差分。
·isOrdered:检查向量是否是顺序的。
·make.index.unique:强制时间唯一，增加毫秒随机数。
·axTicksByTime:计算X轴刻度标记位置按时间描述。
·indexTZ:查询xts对象的时区。
2.2.2 xts包的安装
本节使用的系统环境是：
·Win764bit
·R:3.0.1x86_64-w64-mingw32/x64b4bit
注 xts同时支持Windows 7环境和Linux环境。
xts的安装过程如下：
~ R #
启动R
程序
> install.packages("xts") #
安装xts
包
also installing the dependency 'zoo'
> library(xts) #
加载xts
2.2.3 xts包的使用
1.xts对象的基本操作
查看xts包中的测试数据集sample_matrix。
> data(sample_matrix) #
加载sample_matrix
数据集
> head(sample_matrix) #
查看sample_matrix
数据集的前6
条数据
Open High Low Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
接下来，定义一个xts类型对象。
> sample.xts <- as.xts(sample_matrix, descr='my new xts object')
#
创建一个xts
对象，并设置属性descr
> class(sample.xts) # xts
是继承zoo
类型的对象
[1] "xts" "zoo"
> str(sample.xts) #
打印对象结构
An 'xts' object on 2007-01-02/2007-06-30 containing:
Data: num [1:180, 1:4] 50 50.2 50.4 50.4 50.2 ...
- attr(*, "dimnames")=List of 2
..$ : NULL
..$ : chr [1:4] "Open" "High" "Low" "Close"
Indexed by objects of class: [POSIXct,POSIXt] TZ:
xts Attributes:
List of 1
$ descr: chr "my new xts object"
> attr(sample.xts,'descr') #
查看对象的属性descr
[1] "my new xts object"
在[]中，通过字符串匹配进行xts数据查询。
> head(sample.xts['2007']) #
选出2007
年的数据
Open High Low Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
> head(sample.xts['2007-03/']) #
选出2007
年03
月的数据
Open High Low Close
2007-03-01 50.81620 50.81620 50.56451 50.57075
2007-03-02 50.60980 50.72061 50.50808 50.61559
2007-03-03 50.73241 50.73241 50.40929 50.41033
2007-03-04 50.39273 50.40881 50.24922 50.32636
2007-03-05 50.26501 50.34050 50.26501 50.29567
2007-03-06 50.27464 50.32019 50.16380 50.16380
> head(sample.xts['2007-03-06/2007']) #
选出2007
年03
月06
日到2007
年的数据
Open High Low Close
2007-03-06 50.27464 50.32019 50.16380 50.16380
2007-03-07 50.14458 50.20278 49.91381 49.91381
2007-03-08 49.93149 50.00364 49.84893 49.91839
2007-03-09 49.92377 49.92377 49.74242 49.80712
2007-03-10 49.79370 49.88984 49.70385 49.88698
2007-03-11 49.83062 49.88295 49.76031 49.78806
> sample.xts['2007-01-03'] #
选出2007
年01
月03
日的数据
Open High Low Close
2007-01-03 50.2305 50.42188 50.2305 50.39767
2.用xts对象画图
用xts对象可以画曲线图（图2-8）和K线图（图2-9），下面是产生这两种图的代码，首先是曲线图：
> data(sample_matrix) > plot(as.xts(sample_matrix)) Warning message: In plot.xts(as.xts(sample_matrix)) : only the univariate series will be plotted
警告信息提示，只有单变量序列将被绘制，即只画出第一列数据sample_matrix[,1]的曲线。
图2-8 曲线图
然后是K线图：
> plot(as.xts(sample_matrix), type='candles') # 画K 线图
图2-9 K线图
3.xts对象的类型转换
创建首尾时间函数firstof()和lastof()。
> firstof(2000) # 2000 年的第一天，时分秒显示省略 [1] "2000-01-01 CST" > firstof(2005,01,01) [1] "2005-01-01 CST" > lastof(2007) # 2007 年的最后一天，最后一秒 [1] "2007-12-31 23:59:59.99998 CST" > lastof(2007,10) [1] "2007-10-31 23:59:59.99998 CST"
创建首尾时间。
> .parseISO8601('2000') #
以ISO8601
格式，创建2000
年首尾时间
$first.time
[1] "2000-01-01 CST"
$last.time
[1] "2000-12-31 23:59:59.99998 CST"
> .parseISO8601('2000-05/2001-02')
#
以ISO8601
格式，创建2000
年05
月开始，2001
年02
月结束的时间
$first.time
[1] "2000-05-01 CST"
$last.time
[1] "2001-02-28 23:59:59.99998 CST"
> .parseISO8601('2000-01/02')
$first.time
[1] "2000-01-01 CST"
$last.time
[1] "2000-02-29 23:59:59.99998 CST"
> .parseISO8601('T08:30/T15:00')
$first.time
[1] "1970-01-01 08:30:00 CST"
$last.time
[1] "1970-12-31 15:00:59.99999 CST"
创建以时间类型为索引的xts对象。
> x <- timeBasedSeq('2010-01-01/2010-01-02 12:00') #
创建POSIXt
类型时间
> head(x)
[1] "2010-01-01 00:00:00 CST"
[2] "2010-01-01 00:01:00 CST"
[3] "2010-01-01 00:02:00 CST"
[4] "2010-01-01 00:03:00 CST"
[5] "2010-01-01 00:04:00 CST"
[6] "2010-01-01 00:05:00 CST"
> class(x)
[1] "POSIXt" "POSIXct"
> x <- xts(1:length(x), x) #
以时间为索引创建xts
对象
> head(x)
[,1]
2010-01-01 00:00:00 1
2010-01-01 00:01:00 2
2010-01-01 00:02:00 3
2010-01-01 00:03:00 4
2010-01-01 00:04:00 5
2010-01-01 00:05:00 6
> indexClass(x)
[1] "POSIXt" "POSIXct"
格式化索引时间的显示。
> indexFormat(x) <- "%Y-%b-%d %H:%M:%OS3" #
通过正则格式化索引的时间显示
> head(x)
[,1]
2010-
一月-01 00:00:00.000 1
2010-
一月-01 00:01:00.000 2
2010-
一月-01 00:02:00.000 3
2010-
一月-01 00:03:00.000 4
2010-
一月-01 00:04:00.000 5
2010-
一月-01 00:05:00.000 6
查看索引时间。
> .indexhour(head(x)) # 按小时取索引时间 [1] 0 0 0 0 0 0 > .indexmin(head(x)) # 按分钟取索引时间 [1] 0 1 2 3 4 5
4.xts对象的数据处理
数据对齐。
> x <- Sys.time() + 1:30 > align.time(x, 10) # 整10 秒对齐，秒位为10 的整数倍 [1] "2013-11-18 15:42:30 CST" "2013-11-18 15:42:30 CST" [3] "2013-11-18 15:42:30 CST" "2013-11-18 15:42:40 CST" [5] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST" [7] "2013-11-18 15:42:40 CST" "2013-11-18 15:42:40 CST" [29] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" > align.time(x, 60) # 整60 秒对齐，秒位为0 ，分位为整数 [1] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" [3] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" [5] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" [7] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" [9] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST" [11] "2013-11-18 15:43:00 CST" "2013-11-18 15:43:00 CST"
按时间分割数据，并计算。
> xts.ts <- xts(rnorm(231),as.Date(13514:13744,origin="1970-01-01"))
> apply.monthly(xts.ts,mean) #
按月计算均值，以每月的最后一日显示
[,1]
2007-01-31 0.17699984
2007-02-28 0.30734220
2007-03-31 -0.08757189
2007-04-30 0.18734688
2007-05-31 0.04496954
2007-06-30 0.06884836
2007-07-31 0.25081814
2007-08-19 -0.28845938
> apply.monthly(xts.ts,function(x) var(x))
#
按月计算自定义函数(
方差)
，以每月的最后一日显示
[,1]
2007-01-31 0.9533217
2007-02-28 0.9158947
2007-03-31 1.2821450
2007-04-30 1.2805976
2007-05-31 0.9725438
2007-06-30 1.5228904
2007-07-31 0.8737030
2007-08-19 0.8490521
> apply.quarterly(xts.ts,mean) #
按季计算均值，以每季的最后一日显示
[,1]
2007-03-31 0.12642053
2007-06-30 0.09977926
2007-08-19 0.04589268
> apply.yearly(xts.ts,mean) #
按年计算均值，以年季的最后一日显示
[,1]
2007-08-19 0.09849522
使用to.period()函数按间隔分割数据。
> data(sample_matrix)
> to.period(sample_matrix) #
默认按月分割矩阵数据
sample_matrix.Open sample_matrix.High sample_matrix.Low sample_matrix.Close
2007-01-31 50.03978 50.77336 49.76308 50.22578
2007-02-28 50.22448 51.32342 50.19101 50.77091
2007-03-31 50.81620 50.81620 48.23648 48.97490
2007-04-30 48.94407 50.33781 48.80962 49.33974
2007-05-31 49.34572 49.69097 47.51796 47.73780
2007-06-30 47.74432 47.94127 47.09144 47.76719
> class(to.period(sample_matrix))
[1] "matrix"
> samplexts <- as.xts(sample_matrix) #
默认按月分割xts
类型数据
> to.period(samplexts)
samplexts.Open samplexts.High samplexts.Low samplexts.Close
2007-01-31 50.03978 50.77336 49.76308 50.22578
2007-02-28 50.22448 51.32342 50.19101 50.77091
2007-03-31 50.81620 50.81620 48.23648 48.97490
2007-04-30 48.94407 50.33781 48.80962 49.33974
2007-05-31 49.34572 49.69097 47.51796 47.73780
2007-06-30 47.74432 47.94127 47.09144 47.76719
> class(to.period(samplexts))
[1] "xts" "zoo"
使用endpoints()函数，按间隔分割索引数据。
> data(sample_matrix) > endpoints(sample_matrix) # 默认按月分割 [1] 0 30 58 89 119 150 180 > endpoints(sample_matrix, 'days',k=7) # 按每7 日分割 [1] 0 6 13 20 27 34 41 48 55 62 69 76 83 90 97 104 111 118 125 [20] 132 139 146 153 160 167 174 180 > endpoints(sample_matrix, 'weeks') # 按周分割 [1] 0 7 14 21 28 35 42 49 56 63 70 77 84 9 1 98 105 112 119 126 [20] 133 140 147 154 161 168 175 180 > endpoints(sample_matrix, 'months') # 按月分割 [1] 0 30 58 89 119 150 180
使用merge()函数进行数据合并，按列合并。
> (x <- xts(4:10, Sys.Date()+4:10)) #
创建2
个xts
数据集
[,1]
2013-11-22 4
2013-11-23 5
2013-11-24 6
2013-11-25 7
2013-11-26 8
2013-11-27 9
2013-11-28 10
> (y <- xts(1:6, Sys.Date()+1:6))
[,1]
2013-11-19 1
2013-11-20 2
2013-11-21 3
2013-11-22 4
2013-11-23 5
2013-11-24 6
> merge(x,y) #
按列合并数据，空项以NA
填空
x y
2013-11-19 NA 1
2013-11-20 NA 2
2013-11-21 NA 3
2013-11-22 4 4
2013-11-23 5 5
2013-11-24 6 6
2013-11-25 7 NA
2013-11-26 8 NA
2013-11-27 9 NA
2013-11-28 10 NA
> merge(x,y, join='inner') #
按索引合并数据
x y
2013-11-22 4 4
2013-11-23 5 5
2013-11-24 6 6
> merge(x,y, join='left') #
以左侧为基础合并数据
x y
2013-11-22 4 4
2013-11-23 5 5
2013-11-24 6 6
2013-11-25 7 NA
2013-11-26 8 NA
2013-11-27 9 NA
2013-11-28 10 NA
使用rbind()函数进行数据合并，按行合并。
> x <- xts(1:3, Sys.Date()+1:3)
> rbind(x,x) #
按行合并数据
[,1]
2013-11-19 1
2013-11-19 1
2013-11-20 2
2013-11-20 2
2013-11-21 3
2013-11-21 3
使用split()函数进行数据切片，按行切片。
> data(sample_matrix)
> x <- as.xts(sample_matrix)
> split(x)[[1]] #
默认按月进行切片，打印第一个月的数据
Open High Low Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
2007-01-08 50.03555 50.10363 49.96971 49.98806
> split(x, f="weeks")[[1]] #
按周切片，打印前1
周数据
Open High Low Close
2007-01-02 50.03978 50.11778 49.95041 50.11778
2007-01-03 50.23050 50.42188 50.23050 50.39767
2007-01-04 50.42096 50.42096 50.26414 50.33236
2007-01-05 50.37347 50.37347 50.22103 50.33459
2007-01-06 50.24433 50.24433 50.11121 50.18112
2007-01-07 50.13211 50.21561 49.99185 49.99185
2007-01-08 50.03555 50.10363 49.96971 49.98806
NA值处理。
> x <- xts(1:10, Sys.Date()+1:10)
> x[c(1,2,5,9,10)] <- NA
> x
[,1]
2013-11-19 NA
2013-11-20 NA
2013-11-21 3
2013-11-22 4
2013-11-23 NA
2013-11-24 6
2013-11-25 7
2013-11-26 8
2013-11-27 NA
2013-11-28 NA
> na.locf(x) #
取NA
的前一个，替换NA
值
[,1]
2013-11-19 NA
2013-11-20 NA
2013-11-21 3
2013-11-22 4
2013-11-23 4
2013-11-24 6
2013-11-25 7
2013-11-26 8
2013-11-27 8
2013-11-28 8
> na.locf(x, fromLast=TRUE) #
取NA
后一个，替换NA
值
[,1]
2013-11-19 3
2013-11-20 3
2013-11-21 3
2013-11-22 4
2013-11-23 6
2013-11-24 6
2013-11-25 7
2013-11-26 8
2013-11-27 NA
2013-11-28 NA
5.xts对象的数据统计计算
对xts对象可以进行各种数据统计计算，比如取开始时间和结束时间，计算时间区间，按期间计算统计指标。
（1）取xts对象的开始时间和结束时间，具体代码如下：
> xts.ts <- xts(rnorm(231),as.Date(13514:13744,origin="1970-01-01")) > start(xts.ts) # 取开始时间 [1] "2007-01-01" > end(xts.ts) # 取结束时间 [1] "2007-08-19" > periodicity(xts.ts) # 以日为单位，打印开始和结束时间 Daily periodicity from 2007-01-01 to 2007-08-19
（2）计算时间区间函数，具体代码如下：
> data(sample_matrix) > ndays(sample_matrix) # 计算数据有多少日 [1] 180 > nweeks(sample_matrix) # 计算数据有多少周 [1] 26 > nmonths(sample_matrix) # 计算数据有多少月 [1] 6 > nquarters(sample_matrix) # 计算数据有多少季 [1] 2 > nyears(sample_matrix) # 计算数据有多少年 [1] 1
（3）按期间计算统计指标，具体代码如下：
> zoo.data <- zoo(rnorm(31)+10,as.Date(13514:13744,origin="1970-01-01"))
> ep <- endpoints(zoo.data,'weeks') #
按周获得期间索引
> ep
[1] 0 7 14 21 28 35 42 49 56 63 70 77 84 91 98 105 112 119
[19] 126 133 140 147 154 161 168 175 182 189 196 203 210 217 224 231
> period.apply(zoo.data, INDEX=ep, FUN=function(x) mean(x)) #
计算周的均值
2007-01-07 2007-01-14 2007-01-21 2007-01-28 2007-02-04 2007-02-11 2007-02-18
10.200488 9.649387 10.304151 9.864847 10.382943 9.660175 9.857894
2007-02-25 2007-03-04 2007-03-11 2007-03-18 2007-03-25 2007-04-01 2007-04-08
10.495037 9.569531 10.292899 9.651616 10.089103 9.961048 10.304860
2007-04-15 2007-04-22 2007-04-29 2007-05-06 2007-05-13 2007-05-20 2007-05-27
9.658432 9.887531 10.608082 9.747787 10.052955 9.625730 10.430030
2007-06-03 2007-06-10 2007-06-17 2007-06-24 2007-07-01 2007-07-08 2007-07-15
9.814703 10.224869 9.509881 10.187905 10.229310 10.261725 9.855776
2007-07-22 2007-07-29 2007-08-05 2007-08-12 2007-08-19
9.445072 10.482020 9.844531 10.200488 9.649387
> head(period.max(zoo.data, INDEX=ep)) #
计算周的最大值
[,1]
2007-01-07 12.05912
2007-01-14 10.79286
2007-01-21 11.60658
2007-01-28 11.63455
2007-02-04 12.05912
2007-02-11 10.67887
> head(period.min(zoo.data, INDEX=ep)) #
计算周的最小值
[,1]
2007-01-07 8.874509
2007-01-14 8.534655
2007-01-21 9.069773
2007-01-28 8.461555
2007-02-04 9.421085
2007-02-11 8.534655
> head(period.prod(zoo.data, INDEX=ep)) #
计算周的一个指数值
[,1]
2007-01-07 11140398
2007-01-14 7582350
2007-01-21 11930334
2007-01-28 8658933
2007-02-04 12702505
2007-02-11 7702767
6.xts对象的时间序列操作
检查时间类型。
> class(Sys.time());timeBased(Sys.time()) # Sys.time() 是时间类型POSIXct [1] "POSIXct" "POSIXt" [1] TRUE > class(Sys.Date());timeBased(Sys.Date()) # Sys.Date() 是时间类型Date [1] "Date" [1] TRUE > class(20070101);timeBased(20070101) # 20070101 不是时间类型 [1] "numeric" [1] FALSE
使用timeBasedSeq()函数创建时间序列。
> timeBasedSeq('1999/2008') #
按年
[1] "1999-01-01" "2000-01-01" "2001-01-01" "2002-01-01" "2003-01-01"
[6] "2004-01-01" "2005-01-01" "2006-01-01" "2007-01-01" "2008-01-01"
> head(timeBasedSeq('199901/2008')) #
按月
[1] "
十二月 1998" "
一月 1999" "
二月 1999" "
三月 1999" "
四月 1999"
[6] "
五月 1999"
> head(timeBasedSeq('199901/2008/d'),40) #
按日
[1] "
十二月1998" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[6] "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[11] "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[16] "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[21] "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[26] "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999" "
一月 1999"
[31] "
一月 1999" "
一月 1999" "
二月 1999" "
二月 1999" "
二月 1999"
[36] "
二月 1999" "
二月 1999" "
二月 1999" "
二月 1999" "
二月 1999"
> timeBasedSeq('20080101 0830',length=100) #
按数量创建，100
分钟的数据集
$from
[1] "2008-01-01 08:30:00 CST"
$to
[1] NA
$by
[1] "mins"
$length.out
[1] 100
按索引取数据first()和last()。
> x <- xts(1:100, Sys.Date()+1:100)
> head(x)
[,1]
2013-11-19 1
2013-11-20 2
2013-11-21 3
2013-11-22 4
2013-11-23 5
2013-11-24 6
> first(x, 10) #
取前10
条数据
[,1]
2013-11-19 1
2013-11-20 2
2013-11-21 3
2013-11-22 4
2013-11-23 5
2013-11-24 6
2013-11-25 7
2013-11-26 8
2013-11-27 9
2013-11-28 10
> first(x, '1 day') #
取1
天的数据
[,1]
2013-11-19 1
> last(x, '1 weeks') #
取最后1
周的数据
[,1]
2014-02-24 98
2014-02-25 99
2014-02-26 100
计算步长lag()和差分diff()。
> x <- xts(1:5, Sys.Date()+1:5)
> lag(x) #
以1
为步长
[,1]
2013-11-19 NA
2013-11-20 1
2013-11-21 2
2013-11-22 3
2013-11-23 4
> lag(x, k=-1, na.pad=FALSE) #
以-1
为步长，并去掉NA
值
[,1]
2013-11-19 2
2013-11-20 3
2013-11-21 4
2013-11-22 5
> diff(x) # 1
阶差分
[,1]
2013-11-19 NA
2013-11-20 1
2013-11-21 1
2013-11-22 1
2013-11-23 1
> diff(x,lag=2) # 2
阶差分
[,1]
2013-11-19 NA
2013-11-20 NA
2013-11-21 2
2013-11-22 2
2013-11-23 2
使用isOrdered()函数，检查向量是否排序好的。
> isOrdered(1:10, increasing=TRUE) [1] TRUE > isOrdered(1:10, increasing=FALSE) [1] FALSE > isOrdered(c(1,1:10), increasing=TRUE) [1] FALSE > isOrdered(c(1,1:10), increasing=TRUE, strictly=FALSE) [1] TRUE
使用make.index.unique()函数，强制唯一索引。
> x <- xts(1:5, as.POSIXct("2011-01-21") + c(1,1,1,2,3)/1e3)
> x
[,1]
2011-01-21 00:00:00.000 1
2011-01-21 00:00:00.000 2
2011-01-21 00:00:00.000 3
2011-01-21 00:00:00.002 4
2011-01-21 00:00:00.003 5
> make.index.unique(x) #
增加毫秒级精度，保证索引的唯一性
[,1]
2011-01-21 00:00:00.000999 1
2011-01-21 00:00:00.001000 2
2011-01-21 00:00:00.001001 3
2011-01-21 00:00:00.002000 4
2011-01-21 00:00:00.003000 5
查询xts对象时区。
> x <- xts(1:10, Sys.Date()+1:10) > indexTZ(x) # 时区查询 [1] "UTC" > tzone(x) [1] "UTC" > str(x) An 'xts' object on 2013-11-19/2013-11-28 containing: Data: int [1:10, 1] 1 2 3 4 5 6 7 8 9 10 Indexed by objects of class: [Date] TZ: UTC xts Attributes: NULL
xts给了zoo类型时间序列更多的API支持，这样我们就有了更方便的工具，可以做各种时间序列的转换和变形了。
问题
如何把时间序列可视化？
引言
r-bloggers的一篇博文：plot.xts is wonderful!（http://www.r-bloggers.com/plot-xts-is-wonderful/），让我有动力继续发现xts的强大。xts扩展了zoo的基础数据结构，并提供了更丰富的功能函数。xtsExtra补充库从可视化的角度出发，提供了一个简单而效果非凡的作图函数plot.xts。本节将用plot.xts来演示xts对象的时间序列可视化！
2.3.1 xtsExtra介绍
xtsExtra是xts包的功能补充包，该软件包在Google Summer of Code 2012发布，最终将合并到xts包，不过在笔者写这本书的时候还没有合并。xtsExtra提供的主要功能就是plot.xts()函数。xts::plot.xts()函数与xtsExtra::plot.xts()函数还是有差别的，下面我们就详细介绍其中的差别！
2.3.2 xtsExtra安装
本节使用的系统环境是：
·Win764bit
·R:3.0.1x86_64-w64-mingw32/x64b4bit
注 xtsExtra同时支持Windows 7环境和Linux环境。
由于xtsExtra没有发布到CRAN，我们要从R-Forge下载。
~ R #
启动R
程序
> install.packages("xtsExtra", repos="http://R-Forge.R-project.org")
#
从R-Forge
下载xtsExtra
包
> library(xtsExtra) #
加载xtsExtra
包
xtsExtra::plot.xts()函数覆盖了xts::plot.xts()函数。
2.3.3 xtsExtra包的使用
plot.xts()函数的参数列表如下：
> names(formals(plot.xts)) [1] "x" "y" "screens" "layout.screens" "..." [6] "yax.loc" "auto.grid" "major.ticks" "minor.ticks" "major.format" [11] "bar.col.up" "bar.col.dn" "candle.col" "xy.labels" "xy.lines" [16] "ylim" "panel" "auto.legend" "legend.names" "legend.loc" [21] "legend.pars" "events" "blocks" "nc" "nr"
下面画一个简单的时间序列图，如图2-10所示。
图2-10 时间序列
> data(sample_matrix) > sample_xts <- as.xts(sample_matrix) > plot(sample_xts[,1]) > class(sample_xts[,1]) [1] "xts" "zoo"
从图2-10似乎看不出xtsExtra::plot.xts()函数与xts::plot.xts()函数的不同效果。接下来，我们画点稍微复杂的图形。
1.画K线图
下面就来画K线图，默认红白配色，如图2-11所示。
> plot(sample_xts[1:30, ], type = "candles")
图2-11 K线图
画K线图：自定义颜色，如图2-12所示。
> plot(sample_xts[1:30, ], type = "candles", bar.col.up = "blue", bar.col.dn = "violet", candle.col = "green4")
2.对panel配置
画基本面板，如图2-13所示。
> plot(sample_xts[,1:2])
图2-12 自定义颜色的K线图
图2-13 基本面板
画多行面板，如图2-14所示。
> plot(sample_xts[,rep(1:4, each = 3)])
图2-14 多行面板
画自由组合面板，如图2-15所示。
> plot(sample_xts[,1:4], layout.screens = matrix(c(1,1,1,1,2,3,4,4),ncol = 2, byrow = TRUE))
图2-15 自由组合面板
通过画K线图和面板，就能发现plot.xts()函数提供很多种的画图参数设置，让我们画时间序列图可以有更丰富的可视化表示表现形式。
3.对screens配置
画双屏幕显示，每屏幕2条线，如图2-16所示。
> plot(sample_xts, screens = 1:2)
图2-16 双屏显示
画双屏幕显示，指定曲线出现的屏幕和颜色，如图2-17所示。
> plot(sample_xts, screens = c(1,2,1,2), col = c(1,3,2,2))
图2-17 双屏显示且自定义颜色
画双屏幕显示，指定不同的坐标系，如图2-18所示。
> plot(10^sample_xts, screens = 1:2, log= c("","y"))
图2-18 双屏显示，并设置不同坐标系
画双屏幕显示，指定不同的输出图形，如图2-19所示。
> plot(sample_xts[1:75,1:2] - 50.5, type = c("l","h"), lwd = c(1,2))
图2-19 双屏显示且设置不同的输出图形
画多屏幕，并分组显示，如图2-20所示。
> plot(sample_xts[,c(1:4, 3:4)], layout = matrix(c(1,1,1,1,2,2,3,4,5,6), ncol = 2, byrow = TRUE), yax.loc = "left")
图2-20 多屏显示，并多屏分组
4.对events配置
画基本事件分割线，如图2-21所示。
> plot(sample_xts[,1], events = list(time = c("2007-03-15","2007-05-01"), label
= "bad days"), blocks = list(start.time = c("2007-03-05", "2007-04-15"), end.
time = c("2007-03-20","2007-05-30"), col = c("lightblue1", "lightgreen")))
图2-21 事件分割线
5.双坐标的时间序列
画双坐标视图，如图2-22所示。
> plot(sample_xts[,1],sample_xts[,2])
图2-22 双坐标视图
画双坐标梯度视图，如图2-23所示。
> cr <- colorRampPalette(c("#00FF00","#FF0000"))
> plot(sample_xts[,1],sample_xts[,2], xy.labels = FALSE, xy.lines = TRUE, col =
cr(NROW(sample_xts)), type = "l")
图2-23 双坐标梯度视图
6.对xts类型转换作图
以ts数据类型作图，如图2-24所示。
> tser <- ts(cumsum(rnorm(50, 0.05, 0.15)), start = 2007, frequency = 12) > class(tser) [1] "ts" > plot(tser)
图2-24 ts数据类型作图
以xts类型作图，自动增加了背景坐标线，如图2-25所示。
> plot.xts(tser)
图2-25 xts类型作图
7.用barplot作图
barplot()函数作条形图，如图2-26所示。
> x <- xts(matrix(abs(rnorm(72)), ncol = 6), Sys.Date() + 1:12) > colnames(x) <- LETTERS[1:6] > barplot(x)
图2-26 barplot图
我们看到xtsExtra::plot.xts()函数提供了强大的作图功能，很容易做出包含丰富元素的时间序列！
本章主要介绍了R语言性能相关的3个工具包，帮助读者找到程序中性能的瓶颈。
问题
如何提高R程序性能？
引言
缓存技术在计算机系统中运用得非常广泛，对于高并发访问的应用来说，缓存技术是性价比最高的性能提升方案，特别是对于重复性计算，缓存能为我们节省大量的CPU时间，可能达到99%。R语言以统计计算著名，但其中很多算法包都是在进行大量重复的计算。优化正在进行，改变已经开始。以Hadley Wickham为代表的R语言领军人物，正在让R快起来！
3.1.1 memoise介绍
memoise是一个非常简单的缓存包，以本地缓存为基础，减少单机的重复性计算，从而间接地提升了单机的CPU性能。当用相同的参数对同一个函数执行第二次计算的时候，如果可以直接用第一次计算过的结果作为第二次计算的结果，而不是再重算一遍，那么我们就节省了CPU时间。
memoise包的API非常简单，只有2个函数。
·memoize:定义缓存函数，把一个函数计算的结果加载到本地缓存。
·forget:重置缓存函数，把一个函数的计算结果从本地缓存中删除。
3.1.2 memoise安装
本节使用的系统环境是：
·Win764bit
·R:3.0.1x86_64-w64-mingw32/x64b4bit
注 memoise安装同时支持Windows 7环境和Linux环境。
memoise安装的过程如下：
~ R #
启动R
程序
> install.packages("memoise") #
安装memoise
包
> library(memoise) # memoise
加载
3.1.3 memoise使用
下面我们进行缓存实验，假设一个函数运行时CPU耗时1秒。
> fun <- memoise(function(x) { Sys.sleep(1); runif(1) })#
定义缓存函数，运行时停1
秒
> system.time(print(fun())) #
第一次执行fun
函数，等待1
秒完成
[1] 0.05983416
用户
系统
流逝
0 0 1
> system.time(print(fun())) #
第二次执行fun
函数，从缓存中获得结果，瞬时完成
[1] 0.05983416
用户
系统
流逝
0 0 0
> forget(fun) #
重置缓存函数，清空函数的缓存数据
[1] TRUE
> system.time(print(fun())) #
第三次执行fun
函数，等待1
秒完成
[1] 0.6001663
用户
系统
流逝
0 0 1
> system.time(print(fun())) #
第四次执行fun
函数，从缓存中获得结果，瞬时完成
[1] 0.6001663
用户
系统
流逝
0 0 0
对上面执行过程的具体描述如下：
（1）通过memoise()函数，创建缓存函数fun()。
（2）第一次执行fun()函数，等待1秒。
（3）第二次执行fun()函数，无等待，直接从缓存中获得结果。
（4）通过forget()函数，清空缓存中fun()函数的结果。
（5）第三次执行fun()函数，由于fun()函数的结果被清空，所以重新执行fun()函数，等待1秒。
（6）第四次执行fun()函数，无等待，直接从缓存中返回结果。
3.1.4 memoise()函数源代码分析
memoise项目是一个典型的以面向对象编程方式实现的R语言项目，R语言的面向对象编程将在本书姊妹篇《R的极客理想—高级开发篇》一书中给大家详细介绍。memoise项目的源代码在Gihub上可以找到，项目地址是https://github.com/hadley/memoise。
1.memoise()函数源代码
接下来，我们解读memoise()函数的源代码，memoise()函数主要做了3件事。
（1）new_cache()函数创建新的缓存空间，给f()函数。
（2）生成f()函数的hash值，作为key。
（3）返回缓存后的f函数引用。
memoise <- memoize <- function(f) { #
定义同名函数memoise
和memoize
cache <- new_cache() #
创建新的缓存空间
memo_f <- function(...) { #
生成hash
值
hash <- digest(list(...))
if (cache$has_key(hash)) {
cache$get(hash)
} else {
res <- f(...)
cache$set(hash, res)
res
}
}
attr(memo_f, "memoised") <- TRUE
return(memo_f)
}
2.forget()函数源代码
forget()函数主要做了2件事。
（1）检查环境中是否缓存了f()函数。
（2）如果有f()函数的缓存，则清空f()函数的缓存值。
forget <- function(f) {
if (!is.function(f)) return(FALSE) #
判断函数是否在缓存中
env <- environment(f)
if (!exists("cache", env, inherits = FALSE)) return(FALSE)
cache <- get("cache", env) #
清空函数的缓存值
cache$reset()
TRUE
}
3.私有函数：new_cache()函数源代码
new_cache()函数主要做了3件事。
（1）在new_cache()函数里，定义cache对象，保存在env的环境中。
（2）通过new_cache()函数，构造list类型对象，包括reset,set,get,has_key,keys方法。
（3）通过list对象，对cache对象进行访问操作。
new_cache <- function() {
cache <- NULL #
定义cache
对象
cache_reset <- function() { #
定义清空缓存的函数
cache <<- new.env(TRUE, emptyenv())
}
cache_set <- function(key, value) { #
定义设置缓存的函数
assign(key, value, env = cache)
}
cache_get <- function(key) { #
定义取得缓存的函数
get(key, env = cache, inherits = FALSE)
}
cache_has_key <- function(key) { #
检查缓存是否存在
exists(key, env = cache, inherits = FALSE)
}
cache_reset() #
运行清空缓存
list( #
返回缓存对象
reset = cache_reset,
set = cache_set,
get = cache_get,
has_key = cache_has_key,
keys = function() ls(cache)
)
}
从源代码中，我们不仅能发现memoise包的设计思想，还能感受到作者对R语言深刻的理解。memoise包的代码简洁、高效，值得大家用心学习。
问题
如何找到性能瓶颈？
引言
随着R语言使用越来越深入，R语言的计算性能问题越来越突显。如何能清楚地了解一个算法对CPU的耗时，将成为性能优化的关键因素。值得庆幸的是，R的基础库提供了性能监控的函数Rprof()。
3.2.1 Rprof()函数介绍
Rprof()函数是R语言核心包自带的一个性能数据日志函数，可以打印出函数的调用关系和CPU耗时的数据。然后通过summaryRprof()函数分析Rprof()生成的日志数据，获得性能报告。最后使用profr库的plot()函数，对报告进行可视化。
3.2.2 Rprof()函数的定义
本节使用的系统环境是：
·Win764bit
·R:3.0.1x86_64-w64-mingw32/x64b4bit
注 Rprof和profr同时支持Windows 7环境和Linux环境。
Rprof()函数在基础包utils中定义，不用安装，直接可以使用。查看Rprof()函数的定义。
~ R #
启动R
程序
> Rprof #
查看Rprof()
函数定义
function (filename = "Rprof.out", append = FALSE, interval = 0.02,
memory.profiling = FALSE, gc.profiling = FALSE, line.profiling = FALSE,
numfiles = 100L, bufsize = 10000L)
{
if (is.null(filename))
filename <- ""
invisible(.External(C_Rprof, filename, append, interval,
memory.profiling, gc.profiling, line.profiling, numfiles,
bufsize))
}
<bytecode: 0x000000000d8efda8>
Rprof()函数用来生成记录性能指标的日志文件，通常我们只需要指定输出文件（filename）就可以了。
3.2.3 Rprof()函数使用:股票数据分析案例
以本书6.6节股票数据作为测试数据集，数据文件为000000_0.txt，1.38MB，在随书的代码文件中可以找到。本节只是性能测试，关于数据的业务含义，请参看6.6节内容。
> bidpx1<-read.csv(file="000000_0.txt",header=FALSE)
> names(bidpx1)<-c("tradedate","tradetime","securityid","bidpx1","bidsize1",
"off erpx1", "offersize1")
> bidpx1$securityid<-as.factor(bidpx1$securityid)
> head(bidpx1)
tradedate tradetime securityid bidpx1 bidsize1 offerpx1 offersize1
1 20130724 145004 131810 2.620 6960 2.630 13000
2 20130724 145101 131810 2.860 13880 2.890 6270
3 20130724 145128 131810 2.850 327400 2.851 1500
4 20130724 145143 131810 2.603 44630 2.800 10650
5 20130724 144831 131810 2.890 11400 3.000 77990
6 20130724 145222 131810 2.600 1071370 2.601 35750
> object.size(bidpx1)
1299920 bytes
字段的具体解释：tradedate是交易日期，tradetime是交易时间，securityid是股票ID，bidpx1是买一价，bidsize1是买一量，offerpx1是卖一价，offersize1是卖一量。
计算任务：以securityid分组，计算每小时的买一价的平均值和买一量的总量。
> library(plyr) #
加载plyr
包，用于数据处理
> fun1<-function(){ #
数据处理封装到fun1()
函数
+ datehour<-paste(bidpx1$tradedate,substr(bidpx1$tradetime,1,2),sep="")
+ df<-cbind(datehour,bidpx1[,3:5])
+ ddply(bidpx1,.(securityid,datehour),summarize,price=mean(bidpx1),size=sum(bidsize1))
+ }
> head(fun1())
securityid datehour price size
1 131810 2013072210 3.445549 189670150
2 131810 2013072211 3.437179 131948670
3 131810 2013072212 3.421000 920
4 131810 2013072213 3.509442 299554430
5 131810 2013072214 3.578667 195130420
6 131810 2013072215 1.833000 718940
以system.time()函数查看fun1()函数运行时间，运行2次，系统耗时近似，确定第二次没有缓存。
> system.time(fun1()) 用户 系统 流逝 0.08 0.00 0.07 > system.time(fun1()) 用户 系统 流逝 0.06 0.00 0.06
用Rprof()函数记录性能指标数据，并输出到文件。
> file<-"fun1_rprof.out" # 定义性能日志的输出文件位置 > Rprof(file) # 开始性能监控 > fun1() # 执行计算函数 > Rprof(NULL) # 停止性能监控，并输出到文件
查看生成的性能指标文件：fun1_rprof.out。
~ vi fun1_rprof.out # 用vi 编辑器打开文件 sample.interval=20000 "substr" "paste" "fun1" "paste" "fun1" "structure" "splitter_d" "ddply" "fun1" ".fun" "" ".Call" "loop_apply" "llply" "ldply" "ddply" "fun1" ".fun" "" ".Call" "loop_apply" "llply" "ldply" "ddply" "fun1" ".fun" "" ".Call" "loop_apply" "llply" "ldply" "ddply" "fun1" "[[" "rbind.fill" "list_to_dataframe" "ldply" "ddply" "fun1"
我们其实看不明白这个日志，不知道它到底有什么意义。所以，还要通过summaryRprof()函数来解析这个日志。下面通过summaryRprof()函数查看统计报告。
> summaryRprof(file) #
读入文件
$by.self
self.time self.pct total.time total.pct
".fun" 0.06 42.86 0.06 42.86
"paste" 0.02 14.29 0.04 28.57
"[[" 0.02 14.29 0.02 14.29
"structure" 0.02 14.29 0.02 14.29
"substr" 0.02 14.29 0.02 14.29
$by.total
total.time total.pct self.time self.pct
"fun1" 0.14 100.00 0.00 0.00
"ddply" 0.10 71.43 0.00 0.00
"ldply" 0.08 57.14 0.00 0.00
".fun" 0.06 42.86 0.06 42.86
".Call" 0.06 42.86 0.00 0.00
"" 0.06 42.86 0.00 0.00
"llply" 0.06 42.86 0.00 0.00
"loop_apply" 0.06 42.86 0.00 0.00
"paste" 0.04 28.57 0.02 14.29
"[[" 0.02 14.29 0.02 14.29
"structure" 0.02 14.29 0.02 14.29
"substr" 0.02 14.29 0.02 14.29
"list_to_dataframe" 0.02 14.29 0.00 0.00
"rbind.fill" 0.02 14.29 0.00 0.00
"splitter_d" 0.02 14.29 0.00 0.00
$sample.interval
[1] 0.02
$sampling.time
[1] 0.14
数据解释：
·$by.self是当前函数的耗时情况，其中self.time为实际运行时间，total.time为累计运行时间。
·$by.total是整体函数调用的耗时情况，其中self.time为实际运行时间，total.time为累计运行时间。
从$by.self观察，我们发现时间主要花在.fun函数上。
·.fun:实际运行时间0.06，占当前函数时间的42.86%。
·paste：实际运行时间0.02，占当前函数时间的14.29%。
·"[["：实际运行时间0.02，占当前函数时间的13.29%。
·"structure"：实际运行时间0.02，占当前函数时间的14.29%。
·"substr"：实际运行时间0.02，占当前函数时间的14.29%。
对应在$by.total中，从底到上（1～4），按执行顺序排序。
·4fun1:累计运行时间0.14，占总累计运行时间的100%，实际运行时间0.00。
·3.fun:累计运行时间0.06，占总累计运行时间的42.86%，实际运行时间0.06。
·2paste：累计运行时间0.04，占总累计运行时间的28.57%，实际运行时间0.02。
·1splitter_d：累计运行时间0.02，占总累计运行时间的14.297%，实际运行时间0.00。
这样我们就知道了每个调用函数的CPU时间，进行性能优化就从最耗时的函数开始。
3.2.4 Rprof()函数使用:数据下载案例
我们先安装stockPortfolio包，需要通过stockPortfolio下载股票行情的数据，然后使用Rprof()函数监控下载过程的程序耗时情况。
> install.packages("stockPortfolio") #
安装stockPortfolio
包
> library(stockPortfolio) #
加载stockPortfolio
包
> fileName <- "Rprof2.log"
> Rprof(fileName) #
开始性能监控
> gr <- getReturns(c("GOOG", "MSFT", "IBM"), freq="week")
#
下载Google,
微软,IBM
的股票数据
> gr
Time Scale: week
Average Return
GOOG MSFT IBM
0.004871890 0.001270758 0.001851121
> Rprof(NULL) #
完成性能监控
> summaryRprof(fileName)$by.total[1:8,] #
查看性能报告
total.time total.pct self.time self.pct
"getReturns" 6.76 100.00 0.00 0.00
"read.delim" 6.66 98.52 0.00 0.00
"read.table" 6.66 98.52 0.00 0.00
"scan" 4.64 68.64 4.64 68.64
"file" 2.02 29.88 2.02 29.88
"as.Date" 0.08 1.18 0.02 0.30
"strptime" 0.06 0.89 0.06 0.89
"as.Date.character" 0.06 0.89 0.00 0.00
从打印出的$by.total前8条最耗时的记录来看，实际运行时间(self.time)主要花在file(2.02)和scan(4.64)。
3.2.5 用profr包可视化性能指标
profr包提供了对Rprof()函数输出的性能报告可视化展示功能，比summaryRprof()的纯文字的报告，要更友好，有利于阅读。首先，我们要安装profr包。
> install.packages("profr") #
安装profr
> library(profr) #
加载profr
包
> library(ggplot2) #
用ggplot2
包来画图，
加载ggplot2
包
数据可视化的第一个例子，即股票数据分析案例。下面的代码段产生图3-1和图3-2。
> file<-"fun1_rprof.out"> plot(parse_rprof(file)) # plot 画图 > ggplot(parse_rprof(file)) # ggplot2 画图
图3-1 第一个例子中用plot画的图
数据可视化的第二个例子，即数据下载案例。下面的代码段生成图3-3和图3-4。
> fileName <- "Rprof2.log" > plot(parse_rprof(fileName)) > ggplot(parse_rprof(fileName))
图3-2 第一个例子中用ggplot2画的图
图3-3 第二个例子中用plot画的图
图3-4 第二个例子中用ggplot2画的图
3.2.6 Rprof的命令行使用
Rprof的命令行方法，用来方便地查看日志文件。
1.查看Rprof命令行帮助
~ D:\workspace\R\preforemence\Rprof>R CMD Rprof --help Usage: R CMD Rprof [options] [file] Post-process profiling information in file generated by Rprof(). Options: -h, --help print short help message and exit # 打印帮助信息 -v, --version print version info and exit # 打印版本信息 --lines print line information # 打印显示多行 --total print only by total # 只显示总数 --self print only by self # 只显示私有的 --linesonly print only by line (implies --lines) # 只显示单行( 配合–lines 使用) --min%total= minimum % to print for 'by total' # 显示total 的不低于X 的百分比 --min%self= minimum % to print for 'by self' # 显示self 的不低于X 的百分比
2.Rprof命令行使用
显示完整的报告：
~ D:\workspace\R\preforemence\Rprof>R CMD Rprof fun1_rprof.out Each sample represents 0.02 seconds. Total run time: 0.14 seconds. Total seconds: time spent in function and callees. Self seconds: time spent in function alone. % total % self total seconds self seconds name 100.0 0.14 0.0 0.00 "fun1" 71.4 0.10 0.0 0.00 "ddply" 57.1 0.08 0.0 0.00 "ldply" 42.9 0.06 42.9 0.06 ".fun" 42.9 0.06 0.0 0.00 ".Call" 42.9 0.06 0.0 0.00 "" 42.9 0.06 0.0 0.00 "llply" 42.9 0.06 0.0 0.00 "loop_apply" 28.6 0.04 14.3 0.02 "paste" 14.3 0.02 14.3 0.02 "[[" 14.3 0.02 14.3 0.02 "structure" 14.3 0.02 14.3 0.02 "substr" 14.3 0.02 0.0 0.00 "list_to_dataframe" 14.3 0.02 0.0 0.00 "rbind.fill" 14.3 0.02 0.0 0.00 "splitter_d" % self % total self seconds total seconds name 42.9 0.06 42.9 0.06 ".fun" 14.3 0.02 28.6 0.04 "paste" 14.3 0.02 14.3 0.02 "[[" 14.3 0.02 14.3 0.02 "structure" 14.3 0.02 14.3 0.02 "substr"
只显示total指标，占用时间不低于50%的部分。
~ D:\workspace\R\preforemence\Rprof>R CMD Rprof --total --min%total=50 fun1_rprof.out Each sample represents 0.02 seconds. Total run time: 0.14 seconds. Total seconds: time spent in function and callees. Self seconds: time spent in function alone. % total % self total seconds self seconds name 100.0 0.14 0.0 0.00 "fun1" 71.4 0.10 0.0 0.00 "ddply" 57.1 0.08 0.0 0.00 "ldply"
我们可以通过Rprof进行性能调优，提高代码的性能，让R语言的计算不再是瓶颈。
问题
有没有一个可视化的性能监控工具？
引言
数据可视化越来越受到人们的追捧，图形比文字更有表达力，基于HTML的可交互的图形比静态的PNG图片更让人惊喜。R语言已经为数据可视化做了充分的准备，比如，图形可视化包ggplot2，世界地图可视化包ggmap，股票可视化包quantmod，基于HTML的可交互的可视化包googleVis等，简简单单的几行代码，就可以让数据变成图片，再让图片变成会动的图片。本节以“性能报告”为切入点，讲述R语言可视化包lineprof。
3.3.1 lineprof介绍
lineprof是一个数据可视化的项目，目标是更友好地可视化性能监控效果。3.2节通过profr库，可以把性能数据以图片的形式输出，但仅仅是一张静态图片。而lineprof可以做得更好，生成基于shiny的可交互的网页，让用户自己动手发现问题。
lineprof项目也是Hadley Wickham的作品。目前项目只在Github上面发布，项目地址是https://github.com/hadley/lineprof。lineprof库的API主要有以下两类函数。
（1）功能函数：
·focus:设置显示高度zoom。
·auto_focus:自动设置显示高度zoom。
·lineprof:记录CPU和内存的占用。
·shine:用shiny输出。
（2）内部函数：辅助功能函数
·align:源代码对齐。
·find_ex:用于加载demo。
·line_profile:格式化性能监控数据的输出（Rprof）。
·parse_prof:格式化输出。
·reduce_depth:设置输出的深度。
3.3.2 lineprof安装
本节使用的系统环境是：
·Linux:Ubuntu Server 12.04.2LTS 64bit
·R:3.0.1x86_64-pc-linux-gnu
·IP:192.168.1.201
注 lineprof只支持Linux环境。
由于项目lineprof还没有发布到CRAN，在lineprof包安装时，仅支持从Github安装。所以，需要通过devtools包来安装Github上面发布的R语言项目。
~ R #
启动R
程序
> library(devtools) #
加载devtools
> install_github("lineprof") #
通过devtools
工具，安装linepro
> install_github("shiny-slickgrid", "wch") #
安装制定的shiny
控制台
> library(lineprof) #
加载lineprof
3.3.3 lineprof使用
我们使用官方提供的例子，介绍lineprof包使用。
> source(find_ex("read-delim.r")) #
加载脚本文件read-delim.r
> wine <- find_ex("wine.csv") #
加载测试数据集
> x <- lineprof(read_delim(wine, sep = ","), torture = TRUE) #
执行read_delim
算法，
并通过lineprof
记录性能指标Zooming to read-delim.r (97% of total time)
在例子的资源中，read-delim.r是目标函数的脚本文件，wine.csv是测试数据集，x:lineprof是生成的数据报告。
查看文件:read-delim.r。
function(file, header = TRUE, sep = ",", stringsAsFactors = TRUE) {
# Determine number of fields by reading first line
first <- scan(file, what = character(1), nlines = 1, sep = sep, quiet = TRUE)
p <- length(first)
# Load all fields
all <- scan(file, what = as.list(rep("character", p)), sep = sep,
skip = if (header) 1 else 0, quiet = TRUE)
# Convert from strings to appropriate types
all[] <- lapply(all, type.convert, as.is = !stringsAsFactors)
# Set column names
if (header) {
names(all) <- first
rm(first)
} else {
names(all) <- paste0("V", seq_along(all))
}
# Convert list into data frame
class(all) <- "data.frame"
attr(all, "row.names") <- c(NA_integer_, -length(all[[1]]))
all
}
加载文件wine.csv。
> df<-read.csv(file=wine) > object.size(df) # 数据集占内存大小 20440 bytes > head(df ，3) # 显示前3 行数据 type alcohol malic ash alcalinity magnesium phenols flavanoids nonflavanoids proanthocyanins color hue dilution proline 1 A 14.23 1.71 2.43 15.6 127 2.80 3.06 0.28 2.29 5.64 1.04 3.92 1065 2 A 13.20 1.78 2.14 11.2 100 2.65 2.76 0.26 1.28 4.38 1.05 3.40 1050 3 A 13.16 2.36 2.67 18.6 101 2.80 3.24 0.30 2.81 5.68 1.03 3.17 1185
x对象:lineprof生成的数据报告。
> x #
记录性能指标的对象
Reducing depth to 2 (from 8)
Common path:
time alloc release dups ref src
1 0.002 0.001 0.000 0 #3 read_delim
2 0.049 0.009 0.003 11 #3 read_delim/scan
3 0.026 0.001 0.008 0 #4 read_delim
4 0.379 0.072 0.006 14 #7 read_delim/scan
5 0.003 0.000 0.000 0 #11 read_delim
6 0.106 0.015 0.030 3 #11 read_delim/lapply
7 0.008 0.004 0.000 3 #11 read_delim
8 0.210 0.028 0.077 36 #16 read_delim/rm
9 0.004 0.001 0.000 1 #22 read_delim
10 0.035 0.005 0.004 8 #23 read_delim/[[
11 0.002 0.000 0.000 1 #23 read_delim/length
12 0.001 0.000 0.000 1 #23 read_delim/c
13 0.006 0.004 0.000 1 #23 read_delim
14 0.001 0.000 0.000 0 #23 read_delim/attr<-
用shinySlickgrid库，把性能指标数据可视化，并以网页形式输出。
> library(shinySlickgrid) # 加载shinySlickgrid 库 > shine(x) # 启动shiny 程序 Loading required package: shiny Shiny URLs starting with /lineprof will mapped to /home/conan/R/x86_64-pc-linuxgnu- library/3.0/lineprof/www Shiny URLs starting with /slickgrid will mapped to /home/conan/R/x86_64-pc-linuxgnu- library/3.0/shinySlickgrid/slickgrid Listening on port 6742
shiny会自动在后台打开一个Web服务器，默认的Web访问端口是6742，之后就可以通过浏览器远程访问了。打开浏览器地址栏输入http://192.168.1.201:6742，如图3-5所示。
图3-5 可视化性能监控
在图3-5中，网页的表格一共有6列，其中#是行号，source code是监控的目标函数源代码，t是当前行执行的总（total）时间（秒），r是释放（released）内存量，a是分配（allocated）内存量，d是重复（duplicates）次数。下面给出函数性能数据的解释。
·#6（第6行），用于加载数据，总时间占用0.309s，分配内存0.064mb，方法重复次数为14次
·#15（第15行），用于清除数据，总时间占用0.179s，释放内存0.065mb，方法重复次数为37次
通过lineprof可视化工具，我们生成的报告更灵活，更直观，更好看，支持在网页上进行用户交互，甚至是网页版本的。如果再和3.2节R语言性能监控工具Rprof的效果比较一下，你会惊叹R语言是如此强大，它的进步如此神速。当然，R语言具有无限的潜力，需要更多人来推动它进步，希望你也能给它助力。