• Archive by category "R语言实践"

Blog Archives

超高性能数据处理包data.table

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-data-table/

datatable-title

前言

在R语言中,我们最常用的数据类型是data.frame,绝大多数的数据处理的操作都是围绕着data.frame结构来做的。用data.frame可以很方便的进行数据存储和数据查询,配合apply族函数对数据循环计算,也可也用plyr, reshape2, melt等包对数据实现切分、分组、聚合等的操作。在数据量不太大的时候,使用起来很方便。但是,用data.frame结构处理数据时并不是很高效,特别是在稍大一点数据规模的时候,就会明显变慢。

data.table其实提供了一套和data.frame类似的功能,特别增加了索引的设置,让数据操作非常高效,可能会提升1-2数量级。本章就将data.table包的使用方法。

目录

  1. data.table包介绍
  2. data.table包的使用
  3. data.table包性能对比

1. data.table包介绍

data.table包是一个data.frame的扩展工具集,可以通过自定义keys来设置索引,实现高效的数据索引查询、快速分组、快速连接、快速赋值等数据操作。data.table主要通过二元检索法大大提高数据操作的效率,它也兼容适用于data.frame的向量检索法。同时,data.table对于大数据的快速聚合也有很好的效果,官方介绍说对于 100GB规模内存数据处理,运行效率还是很好的。那么,就让我们试验一下吧。

data.table项目地址:https://cran.r-project.org/web/packages/data.table/

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

data.table包是在CRAN发布的标准库,安装起来非常简单,2条命令就可以了。


~ R
> install.packages("data.table")
> library(data.table)

2. data.table包的使用

接下来,开始用data.table包,并熟悉一下data.table包的基本操作。

2.1 用data.table创建数据集

通常情况,我们用data.frame创建一个数据集时,可以使用下面的语法。


# 创建一个data.frame数据框
> df<-data.frame(a=c('A','B','C','A','A','B'),b=rnorm(6))
> df
  a          b
1 A  1.3847248
2 B  0.6387315
3 C -1.8126626
4 A -0.0265709
5 A -0.3292935
6 B -1.0891958

对于data.table来说,创建一个数据集是和data.frame同样语法。


# 创建一个data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a           b
1: A  0.09174236
2: B -0.84029180
3: C -0.08157873
4: A -0.39992084
5: A -1.66034154
6: B -0.33526447

检查df, dt两个对象的类型,可以看到data.table是对data.frame的扩展类型。


# data.frame类型
> class(df)
[1] "data.frame"

# data.table类型
> class(dt)
[1] "data.table" "data.frame"

如果data.table仅仅是对data.frame的做了S3的扩展类型,那么data.table是不可能做到对data.frame从效率有极大的改进的。为了验证,我们需要检查一下data.table代码的结构定义。


# 打印data.table函数定义
> data.table
function (..., keep.rownames = FALSE, check.names = FALSE, key = NULL) 
{
    x <- list(...)
    if (!.R.listCopiesNamed) 
        .Call(CcopyNamedInList, x)
    if (identical(x, list(NULL)) || identical(x, list(list())) || 
        identical(x, list(data.frame(NULL))) || identical(x, 
        list(data.table(NULL)))) 
        return(null.data.table())
    tt <- as.list(substitute(list(...)))[-1L]
    vnames = names(tt)
    if (is.null(vnames)) 
        vnames = rep.int("", length(x))
    vnames[is.na(vnames)] = ""
    novname = vnames == ""
    if (any(!novname)) {
        if (any(vnames[!novname] == ".SD")) 
            stop("A column may not be called .SD. That has special meaning.")
    }
    for (i in which(novname)) {
        if (is.null(ncol(x[[i]]))) {
            if ((tmp <- deparse(tt[[i]])[1]) == make.names(tmp)) 
                vnames[i] <- tmp
        }
    }
    tt = vnames == ""
    if (any(tt)) 
        vnames[tt] = paste("V", which(tt), sep = "")
    n <- length(x)
    if (n < 1L) 
        return(null.data.table())
    if (length(vnames) != n) 
        stop("logical error in vnames")
    vnames <- as.list.default(vnames)
    nrows = integer(n)
    numcols = integer(n)
    for (i in seq_len(n)) {
        xi = x[[i]]
        if (is.null(xi)) 
            stop("column or argument ", i, " is NULL")
        if ("POSIXlt" %chin% class(xi)) {
            warning("POSIXlt column type detected and converted to POSIXct. We do not recommend use of POSIXlt at all because it uses 40 bytes to store one date.")
            x[[i]] = as.POSIXct(xi)
        }
        else if (is.matrix(xi) || is.data.frame(xi)) {
            xi = as.data.table(xi, keep.rownames = keep.rownames)
            x[[i]] = xi
            numcols[i] = length(xi)
        }
        else if (is.table(xi)) {
            x[[i]] = xi = as.data.table.table(xi, keep.rownames = keep.rownames)
            numcols[i] = length(xi)
        }
        nrows[i] <- NROW(xi)
        if (numcols[i] > 0L) {
            namesi <- names(xi)
            if (length(namesi) == 0L) 
                namesi = rep.int("", ncol(xi))
            namesi[is.na(namesi)] = ""
            tt = namesi == ""
            if (any(tt)) 
                namesi[tt] = paste("V", which(tt), sep = "")
            if (novname[i]) 
                vnames[[i]] = namesi
            else vnames[[i]] = paste(vnames[[i]], namesi, sep = ".")
        }
    }
    nr <- max(nrows)
    ckey = NULL
    recycledkey = FALSE
    for (i in seq_len(n)) {
        xi = x[[i]]
        if (is.data.table(xi) && haskey(xi)) {
            if (nrows[i] < nr) 
                recycledkey = TRUE
            else ckey = c(ckey, key(xi))
        }
    }
    for (i in which(nrows < nr)) {
        xi <- x[[i]]
        if (identical(xi, list())) {
            x[[i]] = vector("list", nr)
            next
        }
        if (nrows[i] == 0L) 
            stop("Item ", i, " has no length. Provide at least one item (such as NA, NA_integer_ etc) to be repeated to match the ", 
                nr, " rows in the longest column. Or, all columns can be 0 length, for insert()ing rows into.")
        if (nr%%nrows[i] != 0L) 
            warning("Item ", i, " is of size ", nrows[i], " but maximum size is ", 
                nr, " (recycled leaving remainder of ", nr%%nrows[i], 
                " items)")
        if (is.data.frame(xi)) {
            ..i = rep(seq_len(nrow(xi)), length.out = nr)
            x[[i]] = xi[..i, , drop = FALSE]
            next
        }
        if (is.atomic(xi) || is.list(xi)) {
            x[[i]] = rep(xi, length.out = nr)
            next
        }
        stop("problem recycling column ", i, ", try a simpler type")
        stop("argument ", i, " (nrow ", nrows[i], ") cannot be recycled without remainder to match longest nrow (", 
            nr, ")")
    }
    if (any(numcols > 0L)) {
        value = vector("list", sum(pmax(numcols, 1L)))
        k = 1L
        for (i in seq_len(n)) {
            if (is.list(x[[i]]) && !is.ff(x[[i]])) {
                for (j in seq_len(length(x[[i]]))) {
                  value[[k]] = x[[i]][[j]]
                  k = k + 1L
                }
            }
            else {
                value[[k]] = x[[i]]
                k = k + 1L
            }
        }
    }
    else {
        value = x
    }
    vnames <- unlist(vnames)
    if (check.names) 
        vnames <- make.names(vnames, unique = TRUE)
    setattr(value, "names", vnames)
    setattr(value, "row.names", .set_row_names(nr))
    setattr(value, "class", c("data.table", "data.frame"))
    if (!is.null(key)) {
        if (!is.character(key)) 
            stop("key argument of data.table() must be character")
        if (length(key) == 1L) {
            key = strsplit(key, split = ",")[[1L]]
        }
        setkeyv(value, key)
    }
    else {
        if (length(ckey) && !recycledkey && !any(duplicated(ckey)) && 
            all(ckey %in% names(value)) && !any(duplicated(names(value)[names(value) %in% 
            ckey]))) 
            setattr(value, "sorted", ckey)
    }
    alloc.col(value)
}
<bytecode: 0x0000000017bfb990>
<environment: namespace:data.table>

从上面的整个大段代码来看,data.table的代码定义中并没有使用data.frame结构的依赖的代码,data.table都在自己函数定义中做的数据处理,所以我们可以确认data.table和data.frame的底层结果是不一样的。

那么为什么从刚刚用class函数检查data.table对象时,会看到data.table和data.frame的扩展关系呢?这里就要了解R语言中对于S3面向对象系统的结构设计了,关于S3的面向对象设计,请参考文章R语言基于S3的面向对象编程

从上面代码中,倒数第17行找到 setattr(value, "class", c("data.table", "data.frame")) 这行,发现这个扩展的定义是作者主动设计的,那么其实就可以理解为,data.table包的作者希望data.table使用起来更像data.frame,所以通过一些包装让使用者无切换成本的。

2.2 data.table和data.frame相互转换

如果想把data.frame对象和data.table对象进行转换,转换的代码是非常容易的,直接转换就可以了。

从一个data.frame对象转型到data.table对象。


# 创建一个data.frame对象
> df<-data.frame(a=c('A','B','C','A','A','B'),b=rnorm(6))

# 检查类型
> class(df)
[1] "data.frame"

# 转型为data.table对象
> df2<-data.table(df)

# 检查类型
> class(df2)
[1] "data.table" "data.frame"

从一个data.table对象转型到data.frame对象。


# 创建一个data.table对象
> dt <- data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))

# 检查类型
> class(dt)
[1] "data.table" "data.frame"

# 转型为data.frame对象
> dt2<-data.frame(dt)

# 检查类型
> class(dt2)
[1] "data.frame"

2.3 用data.table进行查询

由于data.table对用户使用上是希望和data.frame的操作尽量相似,所以适用于data.frame的查询方法基本都适用于data.table,同时data.table自己具有的一些特性,提供了自定义keys来进行高效的查询。

下面先看一下,data.table基本的数据查义方法。


# 创建一个data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a          b
1: A  0.7792728
2: B  1.4870693
3: C  0.9890549
4: A -0.2769280
5: A -1.3009561
6: B  1.1076424

按行或按列查询


# 取第二行的数据
> dt[2,]
   a        b
1: B 1.487069

# 不加,也可以
> dt[2]
   a        b
1: B 1.487069


# 取a列的值
> dt$a
[1] "A" "B" "C" "A" "A" "B"

# 取a列中值为B的行
> dt[a=="B",]
   a        b
1: B 1.487069
2: B 1.107642

# 取a列中值为B的行的判断
> dt[,a=='B']
[1] FALSE  TRUE FALSE FALSE FALSE  TRUE

# 取a列中值为B的行的索引
> which(dt[,a=='B'])
[1] 2 6

上面的操作,不管是用索引值,== 和 $ 都是data.frame操作一样的。下面我们取data.table特殊设计的keys来查询。


# 设置a列为索引列
> setkey(dt,a)

# 打印dt对象,发现数据已经按照a列字母对应ASCII码值进行了排序。
> dt
   a          b
1: A  0.7792728
2: A -0.2769280
3: A -1.3009561
4: B  1.4870693
5: B  1.1076424
6: C  0.9890549

按照自定义的索引进行查询。


# 取a列中值为B的行
> dt["B",]
   a        b
1: B 1.487069
2: B 1.107642

# 取a列中值为B的行,并保留第一行
> dt["B",mult="first"]
   a        b
1: B 1.487069

# 取a列中值为B的行,并保留最后一行
> dt["B",mult="last"]
   a        b
1: B 1.107642

# 取a列中值为b的行,没有数据则为NA
> dt["b"]
   a  b
1: b NA

从上面的代码测试中我们可以看出,在定义了keys后,我们要查询的时候就不用再指定列了,默认会把方括号中的第一位置留给keys,作为索引匹配的查询条件。从代码的角度,又节省了一个变量定义的代码。同时,可以用mult参数,对数据集增加过滤条件,让代码本身也变得更高效。如果查询的值,不是索引列包括的值,则返回NA。

2.4 对data.table对象进行增、删、改操作

给data.table对象增加一列,可以使用这样的格式 data.table[, colname := var1]。


# 创建data.table对象
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a           b
1: A  1.51765578
2: B  0.01182553
3: C  0.71768667
4: A  0.64578235
5: A -0.04210508
6: B  0.29767383

# 增加1列,列名为c
> dt[,c:=b+2]
> dt
   a           b        c
1: A  1.51765578 3.517656
2: B  0.01182553 2.011826
3: C  0.71768667 2.717687
4: A  0.64578235 2.645782
5: A -0.04210508 1.957895
6: B  0.29767383 2.297674

# 增加2列,列名为c1,c2
> dt[,`:=`(c1 = 1:6, c2 = 2:7)]
> dt
   a          b        c c1 c2
1: A  0.7545555 2.754555  1  2
2: B  0.5556030 2.555603  2  3
3: C -0.1080962 1.891904  3  4
4: A  0.3983576 2.398358  4  5
5: A -0.9141015 1.085899  5  6
6: B -0.8577402 1.142260  6  7

# 增加2列,第2种写法
> dt[,c('d1','d2'):=list(1:6,2:7)]
> dt
   a          b        c c1 c2 d1 d2
1: A  0.7545555 2.754555  1  2  1  2
2: B  0.5556030 2.555603  2  3  2  3
3: C -0.1080962 1.891904  3  4  3  4
4: A  0.3983576 2.398358  4  5  4  5
5: A -0.9141015 1.085899  5  6  5  6
6: B -0.8577402 1.142260  6  7  6  7

给data.table对象删除一列时,就是给这列赋值为空,使用这样的格式 data.table[, colname := NULL]。我们继续使用刚才创建的dt对象。


# 删除c1列
> dt[,c1:=NULL]
> dt
   a          b        c c2 d1 d2
1: A  0.7545555 2.754555  2  1  2
2: B  0.5556030 2.555603  3  2  3
3: C -0.1080962 1.891904  4  3  4
4: A  0.3983576 2.398358  5  4  5
5: A -0.9141015 1.085899  6  5  6
6: B -0.8577402 1.142260  7  6  7

# 同时删除d1,d2列
> dt[,c('d1','d2'):=NULL]
> dt
   a          b        c c2
1: A  0.7545555 2.754555  2
2: B  0.5556030 2.555603  3
3: C -0.1080962 1.891904  4
4: A  0.3983576 2.398358  5
5: A -0.9141015 1.085899  6
6: B -0.8577402 1.142260  7

修改data.table对象的值,就是通过索引定位后进行值的替换,通过这样的格式 data.table[condition, colname := 0]。我们继续使用刚才创建的dt对象。


# 给b赋值为30
> dt[,b:=30]
> dt
   a  b        c c2
1: A 30 2.754555  2
2: B 30 2.555603  3
3: C 30 1.891904  4
4: A 30 2.398358  5
5: A 30 1.085899  6
6: B 30 1.142260  7

# 对a列值为B的行,c2列值值大于3的行,的b列赋值为100
> dt[a=='B' & c2>3, b:=100]
> dt
   a   b        c c2
1: A  30 2.754555  2
2: B  30 2.555603  3
3: C  30 1.891904  4
4: A  30 2.398358  5
5: A  30 1.085899  6
6: B 100 1.142260  7

# 还有另一种写法
> dt[,b:=ifelse(a=='B' & c2>3,50,b)]
> dt
   a  b        c c2
1: A 30 2.754555  2
2: B 30 2.555603  3
3: C 30 1.891904  4
4: A 30 2.398358  5
5: A 30 1.085899  6
6: B 50 1.142260  7

2.5 data.table的分组计算

基于data.frame对象做分组计算时,要么使用apply函数自己处理,要么用plyr包的分组计算功能。对于data.table包本身就支持了分组计算,很像SQL的group by这样的功能,这是data.table包主打的优势。

比如,按a列分组,并对b列按分组求和。


# 创建数据
> dt = data.table(a=c('A','B','C','A','A','B'),b=rnorm(6))
> dt
   a          b
1: A  1.4781041
2: B  1.4135736
3: C -0.6593834
4: A -0.1231766
5: A -1.7351749
6: B -0.2528973

# 对整个b列数据求和
> dt[,sum(b)]
[1] 0.1210455

# 按a列分组,并对b列按分组求和
> dt[,sum(b),by=a]
   a         V1
1: A -0.3802474
2: B  1.1606763
3: C -0.6593834

2.6 多个data.table的连接操作

在操作数据的时候,经常会出现2个或多个数据集通过一个索引键进行关联,而我们的算法要把多种数据合并到一起再进行处理,那么这个时候就会用的数据的连接操作,类似关系型数据库的左连接(LEFT JOIN)。

举个例子,学生考试的场景。按照ER设计方法,我们通常会按照实体进行数据划分。这里存在2个实体,一个是学生,一个是成绩。学生实体会包括,学生姓名等的基本资料,而成绩实体会包括,考试的科目,考试的成绩。

假设有6个学生,分别参加A和B两门考试,每门考试得分是不一样的。


# 6个学生
> student <- data.table(id=1:6,name=c('Dan','Mike','Ann','Yang','Li','Kate'));student
   id name
1:  1  Dan
2:  2 Mike
3:  3  Ann
4:  4 Yang
5:  5   Li
6:  6 Kate

# 分别参加A和B两门考试
> score <- data.table(id=1:12,stuId=rep(1:6,2),score=runif(12,60,99),class=c(rep('A',6),rep('B',6)));score
    id stuId    score class
 1:  1     1 89.18497     A
 2:  2     2 61.76987     A
 3:  3     3 74.67598     A
 4:  4     4 64.08165     A
 5:  5     5 85.00035     A
 6:  6     6 95.25072     A
 7:  7     1 81.42813     B
 8:  8     2 82.16083     B
 9:  9     3 69.53405     B
10: 10     4 89.01985     B
11: 11     5 96.77196     B
12: 12     6 97.02833     B

通过学生ID,把学生和考试成绩2个数据集进行连接。


# 设置score数据集,key为stuId
> setkey(score,"stuId")

# 设置student数据集,key为id
> setkey(student,"id")

# 合并两个数据集的数据
> student[score,nomatch=NA,mult="all"]
    id name i.id    score class
 1:  1  Dan    1 89.18497     A
 2:  1  Dan    7 81.42813     B
 3:  2 Mike    2 61.76987     A
 4:  2 Mike    8 82.16083     B
 5:  3  Ann    3 74.67598     A
 6:  3  Ann    9 69.53405     B
 7:  4 Yang    4 64.08165     A
 8:  4 Yang   10 89.01985     B
 9:  5   Li    5 85.00035     A
10:  5   Li   11 96.77196     B
11:  6 Kate    6 95.25072     A
12:  6 Kate   12 97.02833     B

最后我们会看到,两个数据集的结果合并在了一个结果数据集中。这样就完成了,数据连接的操作。从代码的角度来看,1行代码要比用data.frame去拼接方便的多。

3. data.table包性能对比

现在很多时候我们需要处理的数据量是很大的,动辄上百万行甚至上千万行。如果我们要使用R对其进行分析或处理,在不增加硬件的条件下,就需要用一些高性能的数据包进行数据的操作。这里就会发现data.table是非常不错的一个选择。

3.1 data.table和data.frame索引查询性能对比

我们先生成一个稍大数据集,包括2列x和y分别用英文字母进行赋值,100,000,004行,占内存大小1.6G。分别比较data.frame操作和data.table操作的索引查询性能耗时。

使用data.frame创建数据集。


# 清空环境变量
> rm(list=ls())

# 设置大小
> size = ceiling(1e8/26^2)
[1] 147929

# 计算data.frame对象生成的时间 
> t0=system.time(
+   df <- data.frame(x=rep(LETTERS,each=26*size),y=rep(letters,each=size))
+ )

# 打印时间
> t0
用户 系统 流逝 
3.63 0.18 3.80 

# df对象的行数
> nrow(df)
[1] 100000004

# 占用内存
> object.size(df)
1600003336 bytes

# 进行条件查询
> t1=system.time(
+   val1 <- dt[dt$x=="R" & dt$y=="h",]
+ )

# 查询时间
> t1
用户 系统 流逝 
8.53 0.84 9.42 

再使用data.table创建数据集。


# 清空环境变量
> rm(list=ls())

# 设置大小
> size = ceiling(1e8/26^2)
[1] 147929

# 计算data.table对象生成的时间 
> t3=system.time(
+   dt <- data.table(x=rep(LETTERS,each=26*size),y=rep(letters,each=size))
+ )

# 生成对象的时间
> t3
用户 系统 流逝 
3.22 0.39 3.63 

# 对象行数
> nrow(dt)
[1] 100000004

# 占用内存
> object.size(dt)
2000004040 bytes

# 进行条件查询
> t3=system.time(
+ val2 <- dt[x=="R" & y=="h",]
+ )

# 查询时间
> t3
用户 系统 流逝 
6.52 0.26 6.80 

从上面的测试来看,创建对象时,data.table比data.frame显著的高效,而查询效果则并不明显。我们对data.table数据集设置索引,试试有索引查询的效果。


# 设置key索引列为x,y
> setkey(dt,x,y)

# 条件查询
> t4=system.time(
+   val3  <- dt[list("R","h")]
+ )

# 查看时间
> t4
用户 系统 流逝 
0.00 0.00 0.06 

设置索引列后,按索引进行查询,无CPU耗时。震惊了!!

3.2 data.table和data.frame的赋值性能对比

对于赋值操作来说,通常会分为2个动作,先查询再值替换,对于data.frame和data.table都是会按照这个过程来实现的。从上一小节中,可以看到通过索引查询时data.table比data.frame明显的速度要快,对于赋值的操作测试,我们就要最好避免复杂的查询。

对x列值为R的行,对应的y的值进行赋值。首先测试data.frame的计算时间。


> size = 1000000
> df <- data.frame(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   df$y[which(df$x=='R')]<-10
+ )
用户 系统 流逝 
0.75 0.01 0.77 

计算data.table的赋值时间。


> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   dt[x=='R', y:=10]
+ )
用户 系统 流逝 
0.11 0.00 0.11 
> setkey(dt,x)
> system.time(
+   dt['R', y:=10]
+ )
用户 系统 流逝 
0.01 0.00 0.02 

通过对比data.table和data.frame的赋值测试,有索引的data.table性能优势是非常明显的。我们增大数据量,再做一次赋值测试。


> size = 1000000*5
> df <- data.frame(x=rep(LETTERS,each=size),y=rnorm(26*size))
> system.time(
+   df$y[which(df$x=='R')]<-10
+ )
用户 系统 流逝 
3.22 0.25 3.47 

> rm(list=ls())
> size = 1000000*5
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
> system.time(
+   dt['R', y:=10]
+ )
用户 系统 流逝 
0.08 0.01 0.08 

对于增加数据量后data.table,要比data.frame的赋值快更多倍。

3.3 data.table和tapply分组计算性能对比

再对比一下data.table处理数据和tapply的分组计算的性能。测试同样地只做一个简单的计算设定,比如,对一个数据集按x列分组对y列求和。


# 设置数据集大小
> size = 100000
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))

# 设置key为x列
> setkey(dt,x)

# 计算按x列分组,对y列的求和时间
> system.time(
+ r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
0.03 0.00 0.03 

# 用tapply实现,计算求和时间
> system.time(
+ r2<-tapply(dt$y,dt$x,sum)
+ )
用户 系统 流逝 
0.25 0.05 0.30 

# 查看数据集大小, 40mb
> object.size(dt)
41602688 bytes

对于40mb左右的数据来说,tapply比data.table要快,那么我增加数据集的大小,给size*10再测试一下。


> size = 100000*10
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
> val3<-dt[list("R")]
 
> system.time(
+   r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
0.25 0.03 0.28 

> system.time(
+   r2<-tapply(dt$y,dt$x,sum)
+ )
用户 系统 流逝 
2.56 0.36 2.92 

# 400mb数据 
> object.size(dt)
416002688 bytes

对于400mb的数据来说,data.table的计算性能已经明显优于tapply了,再把数据时增加让size*5。


> size = 100000*10*5
> dt <- data.table(x=rep(LETTERS,each=size),y=rnorm(26*size))
> setkey(dt,x)
 
> system.time(
+     r1<-dt[,sum(y),by=x]
+ )
用户 系统 流逝 
1.50 0.11 1.61 

> system.time(
+     r2<-tapply(dt$y,dt$x,sum)
+ )
 用户  系统  流逝 
13.30  3.58 16.90 
 
# 2G数据
> object.size(dt)
2080002688 bytes

对于2G左右的数据来说,tapply总耗时到了16秒,而data.table为1.6秒,从2个的测试来说,大于400mb数据时CPU耗时是线性的。

把上几组测试数据放到一起,下图所示。

data-table

通过上面的对比,我们发现data.table包比tapply快10倍,比data.frame赋值操作快30倍,比data.frame的索引查询快100倍,绝对是值得花精力去学习的一个包。

赶紧用data.table包去优化你的程序吧!

转载请注明出处:
http://blog.fens.me/r-data-table/

打赏作者

掌握R语言中的apply函数族

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-apply/

apply-title

前言

刚开始接触R语言时,会听到各种的R语言使用技巧,其中最重要的一条就是不要用循环,效率特别低,要用向量计算代替循环计算。

那么,这是为什么呢?原因在于R的循环操作for和while,都是基于R语言本身来实现的,而向量操作是基于底层的C语言函数实现的,从性能上来看,就会有比较明显的差距了。那么如何使用C的函数来实现向量计算呢,就是要用到apply的家族函数,包括apply, sapply, tapply, mapply, lapply, rapply, vapply, eapply等。

目录

  1. apply的家族函数
  2. apply函数
  3. lapply函数
  4. sapply函数
  5. vapply函数
  6. mapply函数
  7. tapply函数
  8. rapply函数
  9. eapply函数

1. apply的家族函数

apply函数族是R语言中数据处理的一组核心函数,通过使用apply函数,我们可以实现对数据的循环、分组、过滤、类型控制等操作。但是,由于在R语言中apply函数与其他语言循环体的处理思路是完全不一样的,所以apply函数族一直是使用者玩不转一类核心函数。

很多R语言新手,写了很多的for循环代码,也不愿意多花点时间把apply函数的使用方法了解清楚,最后把R代码写的跟C似得,我严重鄙视只会写for的R程序员。

apply函数本身就是解决数据循环处理的问题,为了面向不同的数据类型,不同的返回值,apply函数组成了一个函数族,包括了8个功能类似的函数。这其中有些函数很相似,有些也不是太一样的。

apply

我一般最常用的函数为apply和sapply,下面将分别介绍这8个函数的定义和使用方法。

2. apply函数

apply函数是最常用的代替for循环的函数。apply函数可以对矩阵、数据框、数组(二维、多维),按行或列进行循环计算,对子元素进行迭代,并把子元素以参数传递的形式给自定义的FUN函数中,并以返回计算结果。

函数定义:

apply(X, MARGIN, FUN, ...)

参数列表:

  • X:数组、矩阵、数据框
  • MARGIN: 按行计算或按按列计算,1表示按行,2表示按列
  • FUN: 自定义的调用函数
  • …: 更多参数,可选

比如,对一个矩阵的每一行求和,下面就要用到apply做循环了。


> x<-matrix(1:12,ncol=3)
> apply(x,1,sum)
[1] 15 18 21 24

下面计算一个稍微复杂点的例子,按行循环,让数据框的x1列加1,并计算出x1,x2列的均值。


# 生成data.frame
> x <- cbind(x1 = 3, x2 = c(4:1, 2:5)); x
     x1 x2
[1,]  3  4
[2,]  3  3
[3,]  3  2
[4,]  3  1
[5,]  3  2
[6,]  3  3
[7,]  3  4
[8,]  3  5

# 自定义函数myFUN,第一个参数x为数据
# 第二、三个参数为自定义参数,可以通过apply的'...'进行传入。
> myFUN<- function(x, c1, c2) {
+   c(sum(x[c1],1), mean(x[c2])) 
+ }

# 把数据框按行做循环,每行分别传递给myFUN函数,设置c1,c2对应myFUN的第二、三个参数
> apply(x,1,myFUN,c1='x1',c2=c('x1','x2'))
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]  4.0    4  4.0    4  4.0    4  4.0    4
[2,]  3.5    3  2.5    2  2.5    3  3.5    4

通过这个上面的自定义函数myFUN就实现了,一个常用的循环计算。

如果直接用for循环来实现,那么代码如下:


# 定义一个结果的数据框
> df<-data.frame()

# 定义for循环
> for(i in 1:nrow(x)){
+   row<-x[i,]                                         # 每行的值
+   df<-rbind(df,rbind(c(sum(row[1],1), mean(row))))   # 计算,并赋值到结果数据框
+ }

# 打印结果数据框
> df
  V1  V2
1  4 3.5
2  4 3.0
3  4 2.5
4  4 2.0
5  4 2.5
6  4 3.0
7  4 3.5
8  4 4.0

通过for循环的方式,也可以很容易的实现上面计算过程,但是这里还有一些额外的操作需要自己处理,比如构建循环体、定义结果数据集、并合每次循环的结果到结果数据集。

对于上面的需求,还有第三种实现方法,那就是完成利用了R的特性,通过向量化计算来完成的。


> data.frame(x1=x[,1]+1,x2=rowMeans(x))
  x1  x2
1  4 3.5
2  4 3.0
3  4 2.5
4  4 2.0
5  4 2.5
6  4 3.0
7  4 3.5
8  4 4.0

那么,一行就可以完成整个计算过程了。

接下来,我们需要再比较一下3种操作上面性能上的消耗。


# 清空环境变量
> rm(list=ls())

# 封装fun1
> fun1<-function(x){
+   myFUN<- function(x, c1, c2) {
+     c(sum(x[c1],1), mean(x[c2])) 
+   }
+   apply(x,1,myFUN,c1='x1',c2=c('x1','x2'))
+ }

# 封装fun2
> fun2<-function(x){
+   df<-data.frame()
+   for(i in 1:nrow(x)){
+     row<-x[i,]
+     df<-rbind(df,rbind(c(sum(row[1],1), mean(row))))
+   }
+ }

# 封装fun3
> fun3<-function(x){
+   data.frame(x1=x[,1]+1,x2=rowMeans(x))
+ }

# 生成数据集
> x <- cbind(x1=3, x2 = c(400:1, 2:500))

# 分别统计3种方法的CPU耗时。
> system.time(fun1(x))
用户 系统 流逝 
0.01 0.00 0.02 

> system.time(fun2(x))
用户 系统 流逝 
0.19 0.00 0.18 

> system.time(fun3(x))
用户 系统 流逝 
   0    0    0 

从CPU的耗时来看,用for循环实现的计算是耗时最长的,apply实现的循环耗时很短,而直接使用R语言内置的向量计算的操作几乎不耗时。通过上面的测试,对同一个计算来说,优先考虑R语言内置的向量计算,必须要用到循环时则使用apply函数,应该尽量避免显示的使用for,while等操作方法。

3. lapply函数

lapply函数是一个最基础循环操作函数之一,用来对list、data.frame数据集进行循环,并返回和X长度同样的list结构作为结果集,通过lapply的开头的第一个字母’l’就可以判断返回结果集的类型。

函数定义:

lapply(X, FUN, ...)

参数列表:

  • X:list、data.frame数据
  • FUN: 自定义的调用函数
  • …: 更多参数,可选

比如,计算list中的每个KEY对应该的数据的分位数。


# 构建一个list数据集x,分别包括a,b,c 三个KEY值。
> x <- list(a = 1:10, b = rnorm(6,10,5), c = c(TRUE,FALSE,FALSE,TRUE));x
$a
 [1]  1  2  3  4  5  6  7  8  9 10
$b
[1]  0.7585424 14.3662366 13.3772979 11.6658990  9.7011387 21.5321427
$c
[1]  TRUE FALSE FALSE  TRUE

# 分别计算每个KEY对应该的数据的分位数。
> lapply(x,fivenum)
$a
[1]  1.0  3.0  5.5  8.0 10.0

$b
[1]  0.7585424  9.7011387 12.5215985 14.3662366 21.5321427

$c
[1] 0.0 0.0 0.5 1.0 1.0

lapply就可以很方便地把list数据集进行循环操作了,还可以用data.frame数据集按列进行循环,但如果传入的数据集是一个向量或矩阵对象,那么直接使用lapply就不能达到想要的效果了。

比如,对矩阵的列求和。


# 生成一个矩阵
> x <- cbind(x1=3, x2=c(2:1,4:5))
> x; class(x)
     x1 x2
[1,]  3  2
[2,]  3  1
[3,]  3  4
[4,]  3  5
[1] "matrix"

# 求和
> lapply(x, sum)
[[1]]
[1] 3

[[2]]
[1] 3

[[3]]
[1] 3

[[4]]
[1] 3

[[5]]
[1] 2

[[6]]
[1] 1

[[7]]
[1] 4

[[8]]
[1] 5

lapply会分别循环矩阵中的每个值,而不是按行或按列进行分组计算。

如果对数据框的列求和。


> lapply(data.frame(x), sum)
$x1
[1] 12

$x2
[1] 12

lapply会自动把数据框按列进行分组,再进行计算。

4. sapply函数

sapply函数是一个简化版的lapply,sapply增加了2个参数simplify和USE.NAMES,主要就是让输出看起来更友好,返回值为向量,而不是list对象。

函数定义:

sapply(X, FUN, ..., simplify=TRUE, USE.NAMES = TRUE)

参数列表:

  • X:数组、矩阵、数据框
  • FUN: 自定义的调用函数
  • …: 更多参数,可选
  • simplify: 是否数组化,当值array时,输出结果按数组进行分组
  • USE.NAMES: 如果X为字符串,TRUE设置字符串为数据名,FALSE不设置

我们还用上面lapply的计算需求进行说明。


> x <- cbind(x1=3, x2=c(2:1,4:5))

# 对矩阵计算,计算过程同lapply函数
> sapply(x, sum)
[1] 3 3 3 3 2 1 4 5

# 对数据框计算
> sapply(data.frame(x), sum)
x1 x2 
12 12 

# 检查结果类型,sapply返回类型为向量,而lapply的返回类型为list
> class(lapply(x, sum))
[1] "list"
> class(sapply(x, sum))
[1] "numeric"

如果simplify=FALSE和USE.NAMES=FALSE,那么完全sapply函数就等于lapply函数了。


> lapply(data.frame(x), sum)
$x1
[1] 12

$x2
[1] 12

> sapply(data.frame(x), sum, simplify=FALSE, USE.NAMES=FALSE)
$x1
[1] 12

$x2
[1] 12

对于simplify为array时,我们可以参考下面的例子,构建一个三维数组,其中二个维度为方阵。


> a<-1:2

# 按数组分组
> sapply(a,function(x) matrix(x,2,2), simplify='array')
, , 1

     [,1] [,2]
[1,]    1    1
[2,]    1    1

, , 2

     [,1] [,2]
[1,]    2    2
[2,]    2    2

# 默认情况,则自动合并分组
> sapply(a,function(x) matrix(x,2,2))
     [,1] [,2]
[1,]    1    2
[2,]    1    2
[3,]    1    2
[4,]    1    2

对于字符串的向量,还可以自动生成数据名。


> val<-head(letters)

# 默认设置数据名
> sapply(val,paste,USE.NAMES=TRUE)
  a   b   c   d   e   f 
"a" "b" "c" "d" "e" "f" 

# USE.NAMES=FALSE,则不设置数据名
> sapply(val,paste,USE.NAMES=FALSE)
[1] "a" "b" "c" "d" "e" "f"

5. vapply函数

vapply类似于sapply,提供了FUN.VALUE参数,用来控制返回值的行名,这样可以让程序更健壮。

函数定义:

vapply(X, FUN, FUN.VALUE, ..., USE.NAMES = TRUE)

参数列表:

  • X:数组、矩阵、数据框
  • FUN: 自定义的调用函数
  • FUN.VALUE: 定义返回值的行名row.names
  • …: 更多参数,可选
  • USE.NAMES: 如果X为字符串,TRUE设置字符串为数据名,FALSE不设置

比如,对数据框的数据进行累计求和,并对每一行设置行名row.names


# 生成数据集
> x <- data.frame(cbind(x1=3, x2=c(2:1,4:5)))

# 设置行名,4行分别为a,b,c,d
> vapply(x,cumsum,FUN.VALUE=c('a'=0,'b'=0,'c'=0,'d'=0))
  x1 x2
a  3  2
b  6  3
c  9  7
d 12 12

# 当不设置时,为默认的索引值
> a<-sapply(x,cumsum);a
     x1 x2
[1,]  3  2
[2,]  6  3
[3,]  9  7
[4,] 12 12

# 手动的方式设置行名
> row.names(a)<-c('a','b','c','d')
> a
  x1 x2
a  3  2
b  6  3
c  9  7
d 12 12

通过使用vapply可以直接设置返回值的行名,这样子做其实可以节省一行的代码,让代码看起来更顺畅,当然如果不愿意多记一个函数,那么也可以直接忽略它,只用sapply就够了。

6. mapply函数

mapply也是sapply的变形函数,类似多变量的sapply,但是参数定义有些变化。第一参数为自定义的FUN函数,第二个参数’…’可以接收多个数据,作为FUN函数的参数调用。

函数定义:

mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE,USE.NAMES = TRUE)

参数列表:

  • FUN: 自定义的调用函数
  • …: 接收多个数据
  • MoreArgs: 参数列表
  • SIMPLIFY: 是否数组化,当值array时,输出结果按数组进行分组
  • USE.NAMES: 如果X为字符串,TRUE设置字符串为数据名,FALSE不设置

比如,比较3个向量大小,按索引顺序取较大的值。


> set.seed(1)

# 定义3个向量
> x<-1:10
> y<-5:-4
> z<-round(runif(10,-5,5))

# 按索引顺序取较大的值。
> mapply(max,x,y,z)
 [1]  5  4  3  4  5  6  7  8  9 10

再看一个例子,生成4个符合正态分布的数据集,分别对应的均值和方差为c(1,10,100,1000)。


> set.seed(1)

# 长度为4
> n<-rep(4,4)

# m为均值,v为方差
> m<-v<-c(1,10,100,1000)

# 生成4组数据,按列分组
> mapply(rnorm,n,m,v)
          [,1]      [,2]      [,3]       [,4]
[1,] 0.3735462 13.295078 157.57814   378.7594
[2,] 1.1836433  1.795316  69.46116 -1214.6999
[3,] 0.1643714 14.874291 251.17812  2124.9309
[4,] 2.5952808 17.383247 138.98432   955.0664

由于mapply是可以接收多个参数的,所以我们在做数据操作的时候,就不需要把数据先合并为data.frame了,直接一次操作就能计算出结果了。

7. tapply函数

tapply用于分组的循环计算,通过INDEX参数可以把数据集X进行分组,相当于group by的操作。

函数定义:

tapply(X, INDEX, FUN = NULL, ..., simplify = TRUE)

参数列表:

  • X: 向量
  • INDEX: 用于分组的索引
  • FUN: 自定义的调用函数
  • …: 接收多个数据
  • simplify : 是否数组化,当值array时,输出结果按数组进行分组

比如,计算不同品种的鸢尾花的花瓣(iris)长度的均值。


# 通过iris$Species品种进行分组
> tapply(iris$Petal.Length,iris$Species,mean)
    setosa versicolor  virginica 
     1.462      4.260      5.552 

对向量x和y进行计算,并以向量t为索引进行分组,求和。


> set.seed(1)

# 定义x,y向量
> x<-y<-1:10;x;y
 [1]  1  2  3  4  5  6  7  8  9 10
 [1]  1  2  3  4  5  6  7  8  9 10

# 设置分组索引t
> t<-round(runif(10,1,100)%%2);t
 [1] 1 2 2 1 1 2 1 0 1 1

# 对x进行分组求和
> tapply(x,t,sum)
 0  1  2 
 8 36 11 

由于tapply只接收一个向量参考,通过’…’可以把再传给你FUN其他的参数,那么我们想去y向量也进行求和,把y作为tapply的第4个参数进行计算。


> tapply(x,t,sum,y)
 0  1  2 
63 91 66 

得到的结果并不符合我们的预期,结果不是把x和y对应的t分组后求和,而是得到了其他的结果。第4个参数y传入sum时,并不是按照循环一个一个传进去的,而是每次传了完整的向量数据,那么再执行sum时sum(y)=55,所以对于t=0时,x=8 再加上y=55,最后计算结果为63。那么,我们在使用’…’去传入其他的参数的时候,一定要看清楚传递过程的描述,才不会出现的算法上的错误。

8. rapply函数

rapply是一个递归版本的lapply,它只处理list类型数据,对list的每个元素进行递归遍历,如果list包括子元素则继续遍历。

函数定义:

rapply(object, f, classes = "ANY", deflt = NULL, how = c("unlist", "replace", "list"), ...)

参数列表:

  • object:list数据
  • f: 自定义的调用函数
  • classes : 匹配类型, ANY为所有类型
  • deflt: 非匹配类型的默认值
  • how: 3种操作方式,当为replace时,则用调用f后的结果替换原list中原来的元素;当为list时,新建一个list,类型匹配调用f函数,不匹配赋值为deflt;当为unlist时,会执行一次unlist(recursive = TRUE)的操作
  • …: 更多参数,可选

比如,对一个list的数据进行过滤,把所有数字型numeric的数据进行从小到大的排序。


> x=list(a=12,b=1:4,c=c('b','a'))
> y=pi
> z=data.frame(a=rnorm(10),b=1:10)
> a <- list(x=x,y=y,z=z)

# 进行排序,并替换原list的值
> rapply(a,sort, classes='numeric',how='replace')
$x
$x$a
[1] 12
$x$b
[1] 4 3 2 1
$x$c
[1] "b" "a"

$y
[1] 3.141593

$z
$z$a
 [1] -0.8356286 -0.8204684 -0.6264538 -0.3053884  0.1836433  0.3295078
 [7]  0.4874291  0.5757814  0.7383247  1.5952808
$z$b
 [1] 10  9  8  7  6  5  4  3  2  1

> class(a$z$b)
[1] "integer"

从结果发现,只有$z$a的数据进行了排序,检查$z$b的类型,发现是integer,是不等于numeric的,所以没有进行排序。

接下来,对字符串类型的数据进行操作,把所有的字符串型加一个字符串’++++’,非字符串类型数据设置为NA。


> rapply(a,function(x) paste(x,'++++'),classes="character",deflt=NA, how = "list")
$x
$x$a
[1] NA
$x$b
[1] NA
$x$c
[1] "b ++++" "a ++++"

$y
[1] NA

$z
$z$a
[1] NA
$z$b
[1] NA

只有$x$c为字符串向量,都合并了一个新字符串。那么,有了rapply就可以对list类型的数据进行方便的数据过滤了。

9. eapply函数

对一个环境空间中的所有变量进行遍历。如果我们有好的习惯,把自定义的变量都按一定的规则存储到自定义的环境空间中,那么这个函数将会让你的操作变得非常方便。当然,可能很多人都不熟悉空间的操作,那么请参考文章 揭开R语言中环境空间的神秘面纱解密R语言函数的环境空间

函数定义:

eapply(env, FUN, ..., all.names = FALSE, USE.NAMES = TRUE)

参数列表:

  • env: 环境空间
  • FUN: 自定义的调用函数
  • …: 更多参数,可选
  • all.names: 匹配类型, ANY为所有类型
  • USE.NAMES: 如果X为字符串,TRUE设置字符串为数据名,FALSE不设置

下面我们定义一个环境空间,然后对环境空间的变量进行循环处理。


# 定义一个环境空间
> env


# 向这个环境空间中存入3个变量
> env$a <- 1:10
> env$beta <- exp(-3:3)
> env$logic <- c(TRUE, FALSE, FALSE, TRUE)
> env


# 查看env空间中的变量
> ls(env)
[1] "a"     "beta"  "logic"

# 查看env空间中的变量字符串结构
> ls.str(env)
a :  int [1:10] 1 2 3 4 5 6 7 8 9 10
beta :  num [1:7] 0.0498 0.1353 0.3679 1 2.7183 ...
logic :  logi [1:4] TRUE FALSE FALSE TRUE

计算env环境空间中所有变量的均值。


> eapply(env, mean)
$logic
[1] 0.5
$beta
[1] 4.535125
$a
[1] 5.5

再计算中当前环境空间中的所有变量的占用内存大小。


# 查看当前环境空间中的变量
> ls()
 [1] "a"     "df"     "env"    "x"     "y"    "z"    "X"  

# 查看所有变量的占用内存大小
> eapply(environment(), object.size)
$a
2056 bytes

$df
1576 bytes

$x
656 bytes

$y
48 bytes

$z
952 bytes

$X
1088 bytes

$env
56 bytes

eapply函数平时很难被用到,但对于R包开发来说,环境空间的使用是必须要掌握的。特别是当R要做为工业化的工具时,对变量的精确控制和管理是非常必要的。

本文全面地介绍了,R语言中的数据循环处理的apply函数族,基本已经可以应对所有的循环处理的情况了。同时,在apply一节中也比较了,3种数据处理方面的性能,R的内置向量计算,要优于apply循环,大幅优于for循环。那么我们在以后的R的开发和使用过程中,应该更多地把apply函数使用好。

忘掉程序员的思维,换成数据的思维,也许你就一下子开朗了。

转载请注明出处:
http://blog.fens.me/r-apply/

打赏作者

R语言高效的管道操作magrittr

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-magrittr/

magrittr

前言

使用R语言进行数据处理是非常方便的,几行代码就可以完成很复杂的操作。但是,对于数据的连续处理,还是有人觉得代码不好看,要么是长长的函数嵌套调用,有点像Lisp感觉,括号包一切;要么就是每次操作赋值一个临时变量,啰嗦。为什么就不能像Linux的管道一样优雅呢?

magrittr包在这样场景中被开发出来,通过管道的方式让连续复杂数据的处理操作,代码更短,更容易读,甚至一行代码可以搞定原来10行代码的事情。

目录

  1. magrittr介绍
  2. magrittr安装
  3. magrittr包的基本使用
  4. magrittr包的扩展功能

1. magrittr介绍

magrittr包被定义为一个高效的管道操作工具包,通过管道的连接方式,让数据或表达式的传递更高效,使用操作符%>%,可以直接把数据传递给下一个函数调用或表达式。magrittr包的主要目标有2个,第一是减少代码开发时间,提高代码的可读性和维护性;第二是让你的代码更短,再短,短短短…

magrittr包,主要定义了4个管道操作符,分另是%>%, %T>%, %$% 和 %<>%。其中,操作符%>%是最常用的,其他3个操作符,与%>%类似,在特殊的使用场景会起到更好的作用。当正确掌握这几个操作符后,你一定会爱不释手的,快去把所有的代码都重构吧,砍掉原来大段冗长的代码是一件多么令人激动的事情啊。

magrittr的项目主页:https://github.com/smbache/magrittr

2. magrittr安装

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

magrittr是在CRAN发布的标准库,安装起来非常简单,2条命令就可以了。


~ R
> install.packages('magrittr')
> library(magrittr)

3. magrittr包的使用

对于magrittr包的使用,其实就是掌握这4个操作符的用法,向右操作符%>%, 向左操作符%T>%, 解释操作符%$% 和 复合赋值操作符%<>%。

3.1 %>% 向右操作符(forward-pipe operator)

%>%是最常用的一个操作符,就是把左侧准备的数据或表达式,传递给右侧的函数调用或表达式进行运行,可以连续操作就像一个链条一样。

现实原理如下图所示,使用%>%把左侧的程序的数据集A传递右侧程序的B函数,B函数的结果数据集再向右侧传递给C函数,最后完成数据计算。

mag1

比如,我们要做下面的事情。(这是一个YY的需求。)

  1. 取10000个随机数符合,符合正态分布。
  2. 求这个10000个数的绝对值,同时乘以50。
  3. 把结果组成一个100*100列的方阵。
  4. 计算方阵中每行的均值,并四舍五入保留到整数。
  5. 把结果除以7求余数,并话出余数的直方图。

我们发现上面的5个过程是连续的,正常的代码我要怎么实现呢。


# 设置随机种子
> set.seed(1)

# 开始 
> n1<-rnorm(10000)            # 第1步
> n2<-abs(n1)*50              # 第2步
> n3<-matrix(n2,ncol = 100)   # 第3步
> n4<-round(rowMeans(n3))     # 第4步
> hist(n4%%7)                 # 第5步

输出的直方图:

01

上面的代码写法是,每一行实现一个条件,但中间多了不少的临时变量。再看另外一种的写法,括号包一切。


# 设置随机种子
> set.seed(1)
> hist(round(rowMeans(matrix(abs(rnorm(10000))*50,ncol=100)))%%7)

输出的直方图:

02

我分别用两种常见的代码风格,实现了我们的需求。再看看%>%的方式,有多么的不一样。


# 设置随机种子
> set.seed(1)

# 开始
> rnorm(10000) %>%
+   abs %>% `*` (50)  %>%
+   matrix(ncol=100)  %>%
+   rowMeans %>% round %>% 
+   `%%`(7) %>% hist

输出的直方图:

03

一行代码,不仅搞定所有的事情,而且结构清楚,可读性非常强。这就是管道代码风格,带来的优雅和简约。

3.2 %T>% 向左操作符(tee operator)

%T>%向左操作符,其实功能和 %>% 基本是一样的,只不过它是把左边的值做为传递的值,而不是右边的值。这种情况的使用场景也是很多的,比如,你在数据处理的中间过程,需要打印输出或图片输出,这时整个过程就会被中断,用向左操作符,就可以解决这样的问题。

现实原理如下图所示,使用%T>%把左侧的程序的数据集A传递右侧程序的B函数,,B函数的结果数据集不再向右侧传递,而是把B左侧的A数据集再次向右传递给C函数,最后完成数据计算。

mag2

我们把上面的需求稍微进行调整,在最后增加一个要求,就会用到向左操作符。

  1. 取10000个随机数符合,符合正态分布。
  2. 求这个10000个数的绝对值,同时乘以50。
  3. 把结果组成一个100*100列的方阵。
  4. 计算方阵中每行的均值,并四舍五入保留到整数。
  5. 把结果除以7求余数,并话出余数的直方图。
  6. 对余数求和

由于输出直方图后,返回值为空,那么再继续管道,就会把空值向右进行传递,这样计算最后一步时就会出错。这时我们需求的是,把除以7的余数向右传递给最后一步求和,那么就可以用到 %T>% 了

直接使用%>%向右传值,出现异常。


> set.seed(1)
> rnorm(10000) %>%
+   abs %>% `*` (50)  %>%
+   matrix(ncol=100)  %>%
+   rowMeans %>% round %>% 
+   `%%`(7) %>% hist %>% sum
Error in sum(.) : invalid 'type' (list) of argument

使用 %T>% 把左边的值,再向右传值,则结果正确。


> rnorm(10000) %>%
+   abs %>% `*` (50)  %>%
+   matrix(ncol=100)  %>%
+   rowMeans %>% round %>% 
+   `%%`(7) %T>% hist %>% sum
[1] 328

3.3 %$% 解释操作符(exposition pipe-operator)

%$% 的作用是把左侧数据的属性名传给右侧,让右侧的调用函数直接通过名字,就可以获取左侧的数据。比如,我们获得一个data.frame类型的数据集,通过使用 %$%,在右侧的函数中可以直接使用列名操作数据。

现实原理如下图所示,使用%$%把左侧的程序的数据集A传递右侧程序的B函数,同时传递数据集A的属性名,作为B函数的内部变量方便对A数据集进行处理,最后完成数据计算。

mag3

下面定义一个3列10行的data.frame,列名分别为x,y,z,或缺x列大于5的数据集。使用 %$% 把列名x直接传到右侧进行判断。这里.代表左侧的完整数据对象。一行代码就实现了需求,而且这里不需要显示的定义中间变量。


> set.seed(1)
> data.frame(x=1:10,y=rnorm(10),z=letters[1:10]) %$% .[which(x>5),]
    x          y z
6   6 -0.8204684 f
7   7  0.4874291 g
8   8  0.7383247 h
9   9  0.5757814 i
10 10 -0.3053884 j

如果不使用%$%,我们通常的代码写法为:


> set.seed(1)
> df<-data.frame(x=1:10,y=rnorm(10),z=letters[1:10])
> df[which(df$x>5),]
    x          y z
6   6 -0.8204684 f
7   7  0.4874291 g
8   8  0.7383247 h
9   9  0.5757814 i
10 10 -0.3053884 j

从代码中可以发现,通常的写法是需要定义变量df的,df一共要被显示的使用3次,就是这一点点的改进,会让代码看起来更干净。

3.4 %<>% 复合赋值操作符(compound assignment pipe-operator)

%<>%复合赋值操作符, 功能与 %>% 基本是一样的,对了一项额外的操作,就是把结果写到左侧对象。比如,我们需要对一个数据集进行排序,那么需要获得排序的结果,用%<>%就是非常方便的。

现实原理如下图所示,使用%<>%把左侧的程序的数据集A传递右侧程序的B函数,B函数的结果数据集再向右侧传递给C函数,C函数结果的数据集再重新赋值给A,完成整个过程。

mag4

定义一个符合正态分布的100个随机数,计算绝对值,并按从小到大的顺序排序,获得并取前10个数字赋值给x。


> set.seed(1)
> x<-rnorm(100) %<>% abs %>% sort %>% head(10)
> x
 [1] 0.001105352 0.016190263 0.028002159 0.039240003 0.044933609 0.053805041 0.056128740
 [8] 0.059313397 0.074341324 0.074564983

是不是太方便了,一行就实现了一连串的操作。但是这里同时有一个陷阱,需要注意一下 %<>% 必须要用在第一个管道的对象处,才能完成赋值的操作,如果不是左侧第一个位置,那么赋值将不起作用。


> set.seed(1)
> x<-rnorm(100)

# 左侧第一个位置,赋值成功
> x %<>% abs %>% sort %>% head(10)
> x
 [1] 0.001105352 0.016190263 0.028002159 0.039240003 0.044933609 0.053805041 0.056128740
 [8] 0.059313397 0.074341324 0.074564983

# 左侧第二个位置,结果被直接打印出来,但是x的值没有变
> x %>% abs %<>% sort %>% head(10)
 [1] 0.001105352 0.016190263 0.028002159 0.039240003 0.044933609 0.053805041 0.056128740
 [8] 0.059313397 0.074341324 0.074564983
> length(x)
[1] 10

# 左侧第三个位置,结果被直接打印出来,但是x的值没有变
> x %>% abs %>% sort %<>% head(10)
 [1] 0.001105352 0.016190263 0.028002159 0.039240003 0.044933609 0.053805041 0.056128740
 [8] 0.059313397 0.074341324 0.074564983
> length(x)
[1] 10

4. magrittr包的扩展功能

我们已经了解了magrittr包的4个操作符的使用,除了操作符,我们再看一下magrittr还有哪些功能。

  • 符号操作符定义
  • %>%对代码块的传递
  • %>%对函数的传递
  • 4.1 符号操作符定义

    为了让链条传递看起来更友好,magrittr对于常见的计算符号操作符进行的重新定义,让每个操作都对应用一个函数,这样所有的传递调用代码都是风格统一的。比如,add()函数和`+`是等价的。

    下面列出对应的列表:

    
    extract	                  `[`
    extract2	          `[[`
    inset	                  `[<-`
    inset2	                  `[[<-`
    use_series	          `$`
    add	                  `+`
    subtract	          `-`
    multiply_by	          `*`
    raise_to_power	          `^`
    multiply_by_matrix	  `%*%`
    divide_by	          `/`
    divide_by_int	          `%/%`
    mod	                  `%%`
    is_in	                  `%in%`
    and	                  `&`
    or	                  `|`
    equals	                  `==`
    is_greater_than	          `>`
    is_weakly_greater_than	  `>=`
    is_less_than	          `<`
    is_weakly_less_than	  `<=`
    not (`n'est pas`)	  `!`
    set_colnames	          `colnames<-`
    set_rownames	          `rownames<-`
    set_names	          `names<-`
    

    我们来看一下使用的效果。对一个包括10个随机数的向量的先*5再+5。

    
    # 使用符号的写法
    > set.seed(1)
    > rnorm(10) %>% `*`(5) %>% `+`(5)
     [1]  1.8677309  5.9182166  0.8218569 12.9764040  6.6475389  0.8976581  7.4371453  8.6916235
     [9]  7.8789068  3.4730581
    
    # 使用函数的写法
    > set.seed(1)
    > rnorm(10) %>% multiply_by(5) %>% add(5)
     [1]  1.8677309  5.9182166  0.8218569 12.9764040  6.6475389  0.8976581  7.4371453  8.6916235
     [9]  7.8789068  3.4730581
    

    上面计算结果是完全一样的,用函数替换了符号。其实,这种转换的操作在面向对象的封装时是非常有用的,像hibernate封装了所有的SQL,XFire封装了WebServices协议等。

    4.2 %>%传递到代码块

    有些时候,我们对同一个数据块的要进行次行的处理,一条语句是很难完成的,这些就需要一个代码块也进行处理。把数据集传递到{}代码块中,传入的数据集以.来表示,通过一段代码来完成操作,而不是一句话完成操作。

    比如,对一个包括10个随机数的向量的先*5再+5,求出向量的均值和标准差,并从小到大排序后返回前5条。

    
    > set.seed(1)
    > rnorm(10)    %>%
    +   multiply_by(5) %>%
    +   add(5)         %>%
    +   { 
    +     cat("Mean:", mean(.), 
    +         "Var:", var(.), "\n")
    +     sort(.) %>% head
    +   }
    Mean: 5.661014 Var: 15.23286 
    [1] 0.8218569 0.8976581 1.8677309 3.4730581 5.9182166 6.6475389
    

    通过{}包装的代码块,就可以很方便的完成多少处理的复杂操作。

    4.3 %>%传递到函数

    传递到函数和传递到代码块设计是类似的,是把一个数据集传给一个匿名函数,进行复杂的数据数据的操作。在这里,我们会显示的定义数据集的名字作为匿名函数的参数。

    比如,对鸢尾花数据集进行处理,只保留第一行和最后一行作为结果。

    
    > iris %>%
    +     (function(x) {
    +         if (nrow(x) > 2) 
    +             rbind(head(x, 1), tail(x, 1))
    +         else x
    +     })
        Sepal.Length Sepal.Width Petal.Length Petal.Width   Species
    1            5.1         3.5          1.4         0.2    setosa
    150          5.9         3.0          5.1         1.8 virginica
    

    这里x就是iris数据集,作为了函数的显示参数,被应用于后续的数据处理过程。

    通过对magrittr的学习,我们掌握了一些特殊的R语言代码的编程技巧,用magrittr包写出的R语言程序,与传统的R语言代码是有区别,可以你的程序很简单、很高效。

    天性“懒惰”的程序员总是会想各种办法,来减少自己的代码,让代码变得优雅,同时还能让程序更可靠。什么时候能把代码写得越来越少,那么你就越来越接近高手!

    转载请注明出处:
    http://blog.fens.me/r-magrittr/

    打赏作者

R语言字符串处理包stringr

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹(Conan), 程序员Java,R,PHP,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-stringr/

stringr

前言

用R语言处理字符串,总觉得很麻烦,即不能用向量的方法进行分割,也不能用循环遍历索引。grep()家族函数常常记不住,paste()函数默认以空格分割,各种不顺手啊!随着使用R语言的场景越来越多,字符串处理是必不可少的。给大家推荐一个由 Hadley Wickham 开发的一个灵活的字符串处理包stringr。

目录

  1. stringr介绍
  2. stringr安装
  3. stringr的API介绍

1. stringr介绍

stringr包被定义为一致的、简单易用的字符串工具集。所有的函数和参数定义都具有一致性,比如,用相同的方法进行NA处理和0长度的向量处理。

字符串处理虽然不是R语言中最主要的功能,却也是必不可少的,数据清洗、可视化等的操作都会用到。对于R语言本身的base包提供的字符串基础函数,随着时间的积累,已经变得很多地方不一致,不规范的命名,不标准的参数定义,很难看一眼就上手使用。字符串处理在其他语言中都是非常方便的事情,R语言在这方面确实落后了。stringr包就是为了解决这个问题,让字符串处理变得简单易用,提供友好的字符串操作接口。

stringr的项目主页:https://cran.r-project.org/web/packages/stringr/index.html

2. stringr安装

本文所使用的系统环境

  • Win10 64bit
  • R: 3.2.3 x86_64-w64-mingw32/x64 b4bit

stringr是在CRAN发布的标准库,安装起来非常简单,2条命令就可以了。


~ R
> install.packages('stringr')
> library(stringr)

3. stringr的API介绍

stringr包1.0.0版本,一共提供了30个函数,方便我们对字符串处理。常用的字符串的处理以str_开头来命名,方便更直观理解函数的定义。我们可以根据使用习惯对函数进行分类:

字符串拼接函数

  • str_c: 字符串拼接。
  • str_join: 字符串拼接,同str_c。
  • str_trim: 去掉字符串的空格和TAB(\t)
  • str_pad: 补充字符串的长度
  • str_dup: 复制字符串
  • str_wrap: 控制字符串输出格式
  • str_sub: 截取字符串
  • str_sub<- 截取字符串,并赋值,同str_sub

字符串计算函数

  • str_count: 字符串计数
  • str_length: 字符串长度
  • str_sort: 字符串值排序
  • str_order: 字符串索引排序,规则同str_sort

字符串匹配函数

  • str_split: 字符串分割
  • str_split_fixed: 字符串分割,同str_split
  • str_subset: 返回匹配的字符串
  • word: 从文本中提取单词
  • str_detect: 检查匹配字符串的字符
  • str_match: 从字符串中提取匹配组。
  • str_match_all: 从字符串中提取匹配组,同str_match
  • str_replace: 字符串替换
  • str_replace_all: 字符串替换,同str_replace
  • str_replace_na:把NA替换为NA字符串
  • str_locate: 找到匹配的字符串的位置。
  • str_locate_all: 找到匹配的字符串的位置,同str_locate
  • str_extract: 从字符串中提取匹配字符
  • str_extract_all: 从字符串中提取匹配字符,同str_extract

字符串变换函数

  • str_conv: 字符编码转换
  • str_to_upper: 字符串转成大写
  • str_to_lower: 字符串转成小写,规则同str_to_upper
  • str_to_title: 字符串转成首字母大写,规则同str_to_upper

参数控制函数,仅用于构造功能的参数,不能独立使用。

  • boundary: 定义使用边界
  • coll: 定义字符串标准排序规则。
  • fixed: 定义用于匹配的字符,包括正则表达式中的转义符
  • regex: 定义正则表达式

3.1 字符串拼接函数

3.1.1 str_c,字符串拼接操作,与str_join完全相同,与paste()行为不完全一致。

函数定义:


str_c(..., sep = "", collapse = NULL)
str_join(..., sep = "", collapse = NULL)

参数列表:

  • …: 多参数的输入
  • sep: 把多个字符串拼接为一个大的字符串,用于字符串的分割符。
  • collapse: 把多个向量参数拼接为一个大的字符串,用于字符串的分割符。

把多个字符串拼接为一个大的字符串。


> str_c('a','b')
[1] "ab"
> str_c('a','b',sep='-')
[1] "a-b"
> str_c(c('a','a1'),c('b','b1'),sep='-')
[1] "a-b"   "a1-b1"

把多个向量参数拼接为一个大的字符串。


> str_c(head(letters), collapse = "")
[1] "abcdef"
> str_c(head(letters), collapse = ", ")
[1] "a, b, c, d, e, f"

# collapse参数,对多个字符串无效
> str_c('a','b',collapse = "-")   
[1] "ab"
> str_c(c('a','a1'),c('b','b1'),collapse='-')
[1] "ab-a1b1"

拼接有NA值的字符串向量时,NA还是NA


> str_c(c("a", NA, "b"), "-d")
[1] "a-d" NA    "b-d"

对比str_c()函数和paste()函数之间的不同点。


# 多字符串拼接,默认的sep参数行为不一致
> str_c('a','b')
[1] "ab"
> paste('a','b')
[1] "a b"

# 向量拼接字符串,collapse参数的行为一致
> str_c(head(letters), collapse = "")
[1] "abcdef"
> paste(head(letters), collapse = "")
[1] "abcdef"
 
#拼接有NA值的字符串向量,对NA的处理行为不一致
> str_c(c("a", NA, "b"), "-d")
[1] "a-d" NA    "b-d"
> paste(c("a", NA, "b"), "-d")
[1] "a -d"  "NA -d" "b -d" 

3.1.2 str_trim:去掉字符串的空格和TAB(\t)

函数定义:

str_trim(string, side = c("both", "left", "right"))

参数列表:

  • string: 字符串,字符串向量。
  • side: 过滤方式,both两边都过滤,left左边过滤,right右边过滤

去掉字符串的空格和TAB(\t)


#只过滤左边的空格
> str_trim("  left space\t\n",side='left') 
[1] "left space\t\n"

#只过滤右边的空格
> str_trim("  left space\t\n",side='right')
[1] "  left space"

#过滤两边的空格
> str_trim("  left space\t\n",side='both')
[1] "left space"

#过滤两边的空格
> str_trim("\nno space\n\t")
[1] "no space"

3.1.3 str_pad:补充字符串的长度

函数定义:

str_pad(string, width, side = c("left", "right", "both"), pad = " ")

参数列表:

  • string: 字符串,字符串向量。
  • width: 字符串填充后的长度
  • side: 填充方向,both两边都填充,left左边填充,right右边填充
  • pad: 用于填充的字符

补充字符串的长度。


# 从左边补充空格,直到字符串长度为20
> str_pad("conan", 20, "left")
[1] "               conan"

# 从右边补充空格,直到字符串长度为20
> str_pad("conan", 20, "right")
[1] "conan               "

# 从左右两边各补充空格,直到字符串长度为20
> str_pad("conan", 20, "both")
[1] "       conan        "

# 从左右两边各补充x字符,直到字符串长度为20
> str_pad("conan", 20, "both",'x')
[1] "xxxxxxxconanxxxxxxxx"

3.1.4 str_dup: 复制字符串

函数定义:

str_dup(string, times)

参数列表:

  • string: 字符串,字符串向量。
  • times: 复制数量

复制一个字符串向量。


> val <- c("abca4", 123, "cba2")

# 复制2次
> str_dup(val, 2)
[1] "abca4abca4" "123123"     "cba2cba2"  

# 按位置复制
> str_dup(val, 1:3)
[1] "abca4"        "123123"       "cba2cba2cba2"

3.1.5 str_wrap,控制字符串输出格式

函数定义:

str_wrap(string, width = 80, indent = 0, exdent = 0)

参数列表:

  • string: 字符串,字符串向量。
  • width: 设置一行所占的宽度。
  • indent: 段落首行的缩进值
  • exdent: 段落非首行的缩进值

 txt<-'R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。'

# 设置宽度为40个字符
> cat(str_wrap(txt, width = 40), "\n")
R语言作为统计学一门语言,一直在小众领域
闪耀着光芒。直到大数据的爆发,R语言变成
了一门炙手可热的数据分析的利器。随着越来
越多的工程背景的人的加入,R语言的社区在
迅速扩大成长。现在已不仅仅是统计领域,教
育,银行,电商,互联网….都在使用R语言。 

# 设置宽度为60字符,首行缩进2字符
> cat(str_wrap(txt, width = 60, indent = 2), "\n")
  R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数
据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来
越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已
不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。 

# 设置宽度为10字符,非首行缩进4字符
> cat(str_wrap(txt, width = 10, exdent = 4), "\n")
R语言作为
    统计学一
    门语言,
    一直在小
    众领域闪
    耀着光芒。
    直到大数据
    的爆发,R
    语言变成了
    一门炙手可
    热的数据分
    析的利器。
    随着越来
    越多的工程
    背景的人的
    加入,R语
    言的社区在
    迅速扩大成
    长。现在已
    不仅仅是统
    计领域,教
    育,银行,
    电商,互联
    网….都在使
    用R语言。 

3.1.6 str_sub,截取字符串

函数定义:

str_sub(string, start = 1L, end = -1L)

参数列表:

  • string: 字符串,字符串向量。
  • start : 开始位置
  • end : 结束位置

截取字符串。


> txt <- "I am Conan."

# 截取1-4的索引位置的字符串
> str_sub(txt, 1, 4)
[1] "I am"

# 截取1-6的索引位置的字符串
> str_sub(txt, end=6)
[1] "I am C"

# 截取6到结束的索引位置的字符串
> str_sub(txt, 6)
[1] "Conan."

# 分2段截取字符串
> str_sub(txt, c(1, 4), c(6, 8))
[1] "I am C" "m Con" 

# 通过负坐标截取字符串
> str_sub(txt, -3)
[1] "an."
> str_sub(txt, end = -3)
[1] "I am Cona"

对截取的字符串进行赋值。


> x <- "AAABBBCCC"

# 在字符串的1的位置赋值为1
> str_sub(x, 1, 1) <- 1; x
[1] "1AABBBCCC"

# 在字符串从2到-2的位置赋值为2345
> str_sub(x, 2, -2) <- "2345"; x
[1] "12345C"

3.2 字符串计算函数

3.2.1 str_count, 字符串计数

函数定义:

str_count(string, pattern = "")

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配的字符。

对字符串中匹配的字符计数


> str_count('aaa444sssddd', "a")
[1] 3

对字符串向量中匹配的字符计数


> fruit <- c("apple", "banana", "pear", "pineapple")
> str_count(fruit, "a")
[1] 1 3 1 1
> str_count(fruit, "p")
[1] 2 0 1 3

对字符串中的'.'字符计数,由于.是正则表达式的匹配符,直接判断计数的结果是不对的。


> str_count(c("a.", ".", ".a.",NA), ".")
[1]  2  1  3 NA

# 用fixed匹配字符
> str_count(c("a.", ".", ".a.",NA), fixed("."))
[1]  1  1  2 NA

# 用\\匹配字符
> str_count(c("a.", ".", ".a.",NA), "\\.")
[1]  1  1  2 NA

3.2.2 str_length,字符串长度

函数定义:

str_length(string)

参数列表:

  • string: 字符串,字符串向量。

计算字符串的长度:


> str_length(c("I", "am", "张丹", NA))
[1]  1  2  2 NA

3.2.3 str_sort, 字符串值排序,同str_order索引排序

函数定义:


str_sort(x, decreasing = FALSE, na_last = TRUE, locale = "", ...)
str_order(x, decreasing = FALSE, na_last = TRUE, locale = "", ...)

参数列表:

  • x: 字符串,字符串向量。
  • decreasing: 排序方向。
  • na_last:NA值的存放位置,一共3个值,TRUE放到最后,FALSE放到最前,NA过滤处理
  • locale:按哪种语言习惯排序

对字符串值进行排序。


# 按ASCII字母排序
> str_sort(c('a',1,2,'11'), locale = "en")  
[1] "1"  "11" "2"  "a" 

# 倒序排序
> str_sort(letters,decreasing=TRUE)         
 [1] "z" "y" "x" "w" "v" "u" "t" "s" "r" "q" "p" "o" "n" "m" "l" "k" "j" "i" "h"
[20] "g" "f" "e" "d" "c" "b" "a"

# 按拼音排序
> str_sort(c('你','好','粉','丝','日','志'),locale = "zh")  
[1] "粉" "好" "你" "日" "丝" "志"

对NA值的排序处理


 #把NA放最后面
> str_sort(c(NA,'1',NA),na_last=TRUE) 
[1] "1" NA  NA
 
#把NA放最前面
> str_sort(c(NA,'1',NA),na_last=FALSE) 
[1] NA  NA  "1"

#去掉NA值 
> str_sort(c(NA,'1',NA),na_last=NA)    
[1] "1"

3.3 字符串匹配函数

3.3.1 str_split,字符串分割,同str_split_fixed

函数定义:


str_split(string, pattern, n = Inf)
str_split_fixed(string, pattern, n)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配的字符。
  • n: 分割个数

对字符串进行分割。


> val <- "abc,123,234,iuuu"

# 以,进行分割
> s1<-str_split(val, ",");s1
[[1]]
[1] "abc"  "123"  "234"  "iuuu"

# 以,进行分割,保留2块
> s2<-str_split(val, ",",2);s2
[[1]]
[1] "abc"          "123,234,iuuu"

# 查看str_split()函数操作的结果类型list
> class(s1)
[1] "list"

# 用str_split_fixed()函数分割,结果类型是matrix
> s3<-str_split_fixed(val, ",",2);s3
     [,1]  [,2]          
[1,] "abc" "123,234,iuuu"

> class(s3)
[1] "matrix"

3.3.2 str_subset:返回的匹配字符串

函数定义:

str_subset(string, pattern)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配的字符。

> val <- c("abc", 123, "cba")

# 全文匹配
> str_subset(val, "a")
[1] "abc" "cba"

# 开头匹配
> str_subset(val, "^a")
[1] "abc"

# 结尾匹配
> str_subset(val, "a$")
[1] "cba"

3.3.3 word, 从文本中提取单词

函数定义:

word(string, start = 1L, end = start, sep = fixed(" "))

参数列表:

  • string: 字符串,字符串向量。
  • start: 开始位置。
  • end: 结束位置。
  • sep: 匹配字符。

> val <- c("I am Conan.", "http://fens.me, ok")

# 默认以空格分割,取第一个位置的字符串
> word(val, 1)
[1] "I"               "http://fens.me,"
> word(val, -1)
[1] "Conan." "ok"    
> word(val, 2, -1)
[1] "am Conan." "ok"       

# 以,分割,取第一个位置的字符串 
> val<-'111,222,333,444'
> word(val, 1, sep = fixed(','))
[1] "111"
> word(val, 3, sep = fixed(','))
[1] "333"

3.3.4 str_detect匹配字符串的字符

函数定义:

str_detect(string, pattern)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配字符。

> val <- c("abca4", 123, "cba2")

# 检查字符串向量,是否包括a
> str_detect(val, "a")
[1]  TRUE FALSE  TRUE

# 检查字符串向量,是否以a为开头
> str_detect(val, "^a")
[1]  TRUE FALSE FALSE

# 检查字符串向量,是否以a为结尾
> str_detect(val, "a$")
[1] FALSE FALSE FALSE

3.3.6 str_match,从字符串中提取匹配组

函数定义:


str_match(string, pattern)
str_match_all(string, pattern)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配字符。

从字符串中提取匹配组


> val <- c("abc", 123, "cba")

# 匹配字符a,并返回对应的字符
> str_match(val, "a")
     [,1]
[1,] "a" 
[2,] NA  
[3,] "a" 

# 匹配字符0-9,限1个,并返回对应的字符
> str_match(val, "[0-9]")
     [,1]
[1,] NA  
[2,] "1" 
[3,] NA  

# 匹配字符0-9,不限数量,并返回对应的字符
> str_match(val, "[0-9]*")
     [,1] 
[1,] ""   
[2,] "123"
[3,] ""  

从字符串中提取匹配组,以字符串matrix格式返回


> str_match_all(val, "a")
[[1]]
     [,1]
[1,] "a" 

[[2]]
     [,1]

[[3]]
     [,1]
[1,] "a" 

> str_match_all(val, "[0-9]")
[[1]]
     [,1]

[[2]]
     [,1]
[1,] "1" 
[2,] "2" 
[3,] "3" 

[[3]]
     [,1]

3.3.7 str_replace,字符串替换

函数定义:

str_replace(string, pattern, replacement)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配字符。
  • replacement: 用于替换的字符。

> val <- c("abc", 123, "cba")

# 把目标字符串第一个出现的a或b,替换为-
> str_replace(val, "[ab]", "-")
[1] "-bc" "123" "c-a"

# 把目标字符串所有出现的a或b,替换为-
> str_replace_all(val, "[ab]", "-")
[1] "--c" "123" "c--"

# 把目标字符串所有出现的a,替换为被转义的字符
> str_replace_all(val, "[a]", "\1\1")
[1] "\001\001bc" "123"        "cb\001\001"

3.3.8 str_replace_na把NA替换为NA字符串

函数定义:

str_replace_na(string, replacement = "NA")

参数列表:

  • string: 字符串,字符串向量。
  • replacement : 用于替换的字符。

把NA替换为字符串


> str_replace_na(c(NA,'NA',"abc"),'x')
[1] "x"   "NA"  "abc"

3.3.9 str_locate,找到的模式在字符串中的位置。

函数定义:

str_locate(string, pattern)
str_locate_all(string, pattern)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配字符。

> val <- c("abca", 123, "cba")

# 匹配a在字符串中的位置
> str_locate(val, "a")
     start end
[1,]     1   1
[2,]    NA  NA
[3,]     3   3

# 用向量匹配
> str_locate(val, c("a", 12, "b"))
     start end
[1,]     1   1
[2,]     1   2
[3,]     2   2

# 以字符串matrix格式返回
> str_locate_all(val, "a")
[[1]]
     start end
[1,]     1   1
[2,]     4   4

[[2]]
     start end

[[3]]
     start end
[1,]     3   3

# 匹配a或b字符,以字符串matrix格式返回
> str_locate_all(val, "[ab]")
[[1]]
     start end
[1,]     1   1
[2,]     2   2
[3,]     4   4

[[2]]
     start end

[[3]]
     start end
[1,]     2   2
[2,]     3   3

3.3.10 str_extract从字符串中提取匹配模式

函数定义:

str_extract(string, pattern)
str_extract_all(string, pattern, simplify = FALSE)

参数列表:

  • string: 字符串,字符串向量。
  • pattern: 匹配字符。
  • simplify: 返回值,TRUE返回matrix,FALSE返回字符串向量

> val <- c("abca4", 123, "cba2")

# 返回匹配的数字
> str_extract(val, "\\d")
[1] "4" "1" "2"

# 返回匹配的字符
> str_extract(val, "[a-z]+")
[1] "abca" NA     "cba" 


> val <- c("abca4", 123, "cba2")
> str_extract_all(val, "\\d")
[[1]]
[1] "4"

[[2]]
[1] "1" "2" "3"

[[3]]
[1] "2"

> str_extract_all(val, "[a-z]+")
[[1]]
[1] "abca"

[[2]]
character(0)

[[3]]
[1] "cba"

3.4 字符串变换函数

3.4.1 str_conv:字符编码转换

函数定义:

str_conv(string, encoding)

参数列表:

  • string: 字符串,字符串向量。
  • encoding: 编码名。

对中文进行转码处理。


# 把中文字符字节化
> x <- charToRaw('你好');x
[1] c4 e3 ba c3

# 默认win系统字符集为GBK,GB2312为GBK字集,转码正常
> str_conv(x, "GBK")
[1] "你好"
> str_conv(x, "GB2312")
[1] "你好"

# 转UTF-8失败
> str_conv(x, "UTF-8")
[1] "���"
Warning messages:
1: In stri_conv(string, encoding, "UTF-8") :
  input data \xffffffc4 in current source encoding could not be converted to Unicode
2: In stri_conv(string, encoding, "UTF-8") :
  input data \xffffffe3\xffffffba in current source encoding could not be converted to Unicode
3: In stri_conv(string, encoding, "UTF-8") :
  input data \xffffffc3 in current source encoding could not be converted to Unicode

把unicode转UTF-8


> x1 <- "\u5317\u4eac"
> str_conv(x1, "UTF-8")
[1] "北京"

3.4.2 str_to_upper,字符串大写转换。

函数定义:


str_to_upper(string, locale = "")
str_to_lower(string, locale = "")
str_to_title(string, locale = "")

参数列表:

  • string: 字符串。
  • locale:按哪种语言习惯排序

字符串大写转换:


> val <- "I am conan. Welcome to my blog! http://fens.me"

# 全大写
> str_to_upper(val)
[1] "I AM CONAN. WELCOME TO MY BLOG! HTTP://FENS.ME"

# 全小写
> str_to_lower(val)
[1] "i am conan. welcome to my blog! http://fens.me"

# 首字母大写
> str_to_title(val)
[1] "I Am Conan. Welcome To My Blog! Http://Fens.Me"

字符串在平常的数据处理中经常用过,需要对字符串进行分割、连接、转换等操作,本篇中通过介绍stringr,灵活的字符串处理库,可以有效地提高代码的编写效率。有了好的工具,在用R语言处理字符串就顺手了。

转载请注明出处:
http://blog.fens.me/r-stringr/

打赏作者

R语言构建配对交易量化模型

用IT技术玩金融系列文章,将介绍如何使用IT技术,处理金融大数据。在互联网混迹多年,已经熟练掌握一些IT技术。单纯地在互联网做开发,总觉得使劲的方式不对。要想靠技术养活自己,就要把技术变现。通过“跨界”可以寻找新的机会,创造技术的壁垒。

金融是离钱最近的市场,也是变现的好渠道!今天就开始踏上“用IT技术玩金融”之旅!

关于作者:

  • 张丹(Conan), 程序员Java,R,Javascript
  • weibo:@Conan_Z
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/finance-pairs-trading/

pair-trading

前言

散户每天都在经历中国股市的上蹿下跳,赚到钱是运气,赔钱是常态。那么是否有方法可以让赚钱变成常态呢?

我们可以通过“统计套利”的方法,发现市场的无效性。配对交易,就统计套利策略的一种,通过对冲掉绝大部分的市场风险,抓住套利机会,积累小盈利汇聚大收益。

目录

  1. 什么是配对交易?
  2. 配对交易的模型
  3. 用R语言实现配对交易

1. 什么是配对交易?

配对交易(Pairs Trading)的理念最早来源于上世纪20年代华尔街传奇交易员Jesse Livermore 的姐妹股票对交易策略。配对交易的基本原理是找到两个相关性较高具备均衡关系的股票或其他金融产品,做空近期相对强势的金融产品,同时做多相对弱势金融产品,等待两者价格重返均衡值时进行平仓,赚取两者的价差变动的收益。

假设两个金融产品在未来的时期会保持良好的均衡关系,一旦两者之间的价格走势出现背离,同时这种背离在未来会被进行修复,那么就可能产生套利的机会。对于配对交易来说,就是找到这样的机会,进行统计套利。

配对交易的特点

配对交易与传统股票交易最大的不同之处在于,它的投资标的是两只股票的价差,是一种相对价值而非绝对价值。由于它在股票多头和空头方同时建仓,对冲掉了绝大部分的市场风险,所以它是一种市场的中性策略。无论大盘上涨还是下跌,配对交易策略收益都是相对平稳的,与大盘走势的相关性很低。

在市场无趋势性机会时,可以通过配对交易避免股市系统风险,获取Alpha绝对收益。趋势性的交易策略,可以参考文章 两条均线打天下

配对交易操作方法

  1. 组合筛选:在市场上寻找用于配对的金融产品或者组合,检查历史价格的走势,判断是否可以用来进行配对。主要用下面几个指标来筛选配对组合:相关系数、模型计算的均值回复速度、协整检验、基本面因素等。通过这些因素来寻找出具有稳定相关关系的组合。
  2. 风险衡量和动态组合的构建:计算配对组合各自的预期收益、预期风险、交易成本;判断两个组合之间的价差服从何种分布;判断是具有长期均衡特性还是短期均衡特性;价差发生跳跃的频率等。
  3. 确定交易规则:根据价差的特性,确定交易的频率(高频交易还是低频交易),交易的触发条件和平仓规则等。
  4. 执行交易及风险控制:除了按照交易规则执行外,还必须动态跟踪价差走势,如果发现突变,应该及时调整套利模式和交易频率。

配对交易缺点

  • 统计套利的规则都是基于历史数据计算的,但历史不能代表未来,当市场发生变化模型也会失效
  • 市场对价格进行修复的时间难以准确判断,只能根据历史大致估计。如果回归的时间过长,对套利者的资金使用成本是个考验,也有可能导致套利失败。

2. 构建配对交易的模型

根据配对交易的原理,我们就可以自己设计配对交易的模型了。首先,需要把配对交易涉及的指标都进行量化,比如如何选择不同的两个具备均衡关系金融产品,什么时候做多,什么时候做空,什么时候平仓等。

根据概念,我们生成两个虚拟的金融产品X,Y,包括时间和价格字段。让X和Y的两个产品都价格符合正态分布,生成100个日期的数据。由于是测试程序,日期字段是包括了自然日,暂时理解为连续的日期。

R语言实现的代码如下:


> set.seed(1)                         #设置随机种子
> dates<-as.Date('2010-01-01')+1:100  #100个日期
> x<-round(rnorm(100,50,40),2)        #随机生成X产品,100个正态分析的收盘价 
> y<-round(rnorm(100,50,40),2)        #随机生成Y产品,100个正态分析的收盘价 
> df<-data.frame(dates,x,y)
> df
         dates      x      y
1   2010-01-02  24.94  25.19
2   2010-01-03  57.35  51.68
3   2010-01-04  16.57  13.56
4   2010-01-05 113.81  56.32
5   2010-01-06  63.18  23.82
6   2010-01-07  17.18 120.69
7   2010-01-08  69.50  78.67
8   2010-01-09  79.53  86.41
9   2010-01-10  73.03  65.37
10  2010-01-11  37.78 117.29
11  2010-01-12 110.47  24.57
12  2010-01-13  65.59  31.53
13  2010-01-14  25.15 107.29
14  2010-01-15 -38.59  23.97
15  2010-01-16  95.00  41.70
16  2010-01-17  48.20  34.29
17  2010-01-18  49.35  37.20
18  2010-01-19  87.75  38.84
19  2010-01-20  82.85  69.77
20  2010-01-21  73.76  42.91
21  2010-01-22  86.76  29.76
22  2010-01-23  81.29 103.72
23  2010-01-24  52.98  41.42
24  2010-01-25 -29.57  42.82
25  2010-01-26  74.79  45.99
26  2010-01-27  47.75  78.51
27  2010-01-28  43.77  47.06
28  2010-01-29  -8.83  48.49
29  2010-01-30  30.87  22.73
30  2010-01-31  66.72  37.03
31  2010-02-01 104.35  52.41
32  2010-02-02  45.89  26.44
33  2010-02-03  65.51  71.26
34  2010-02-04  47.85 -10.74
35  2010-02-05  -5.08  62.26
36  2010-02-06  33.40 -11.46
37  2010-02-07  34.23  37.96
38  2010-02-08  47.63  28.87
39  2010-02-09  94.00  23.92
40  2010-02-10  80.53  47.72
41  2010-02-11  43.42 -26.57
42  2010-02-12  39.87  97.06
43  2010-02-13  77.88 -16.60
44  2010-02-14  72.27  31.46
45  2010-02-15  22.45   5.36
46  2010-02-16  21.70  19.97
47  2010-02-17  64.58 133.49
48  2010-02-18  80.74  50.70
49  2010-02-19  45.51  -1.45
50  2010-02-20  85.24 -15.62
51  2010-02-21  65.92  68.01
52  2010-02-22  25.52  49.26
53  2010-02-23  63.64  37.28
54  2010-02-24   4.83  12.83
55  2010-02-25 107.32  -9.50
56  2010-02-26 129.22   6.99
57  2010-02-27  35.31  90.00
58  2010-02-28   8.23  25.15
59  2010-03-01  72.79  -5.38
60  2010-03-02  44.60 124.77
61  2010-03-03 146.06  67.00
62  2010-03-04  48.43  40.45
63  2010-03-05  77.59  92.34
64  2010-03-06  51.12  85.46
65  2010-03-07  20.27  25.23
66  2010-03-08  57.55 138.24
67  2010-03-09 -22.20  39.80
68  2010-03-10 108.62  -6.98
69  2010-03-11  56.13  44.22
70  2010-03-12 136.90  58.30
71  2010-03-13  69.02 142.32
72  2010-03-14  21.60  54.23
73  2010-03-15  74.43  68.28
74  2010-03-16  12.64  46.91
75  2010-03-17  -0.15  36.64
76  2010-03-18  61.66  48.61
77  2010-03-19  32.27  81.51
78  2010-03-20  50.04 133.01
79  2010-03-21  52.97  91.10
80  2010-03-22  26.42  98.32
81  2010-03-23  27.25   0.75
82  2010-03-24  44.59  89.36
83  2010-03-25  97.12  58.80
84  2010-03-26 -10.94  -8.69
85  2010-03-27  73.76  70.84
86  2010-03-28  63.32  43.65
87  2010-03-29  92.52 108.58
88  2010-03-30  37.83  19.36
89  2010-03-31  64.80  32.79
90  2010-04-01  60.68  12.96
91  2010-04-02  28.30  42.92
92  2010-04-03  98.31  66.08
93  2010-04-04  96.42  20.73
94  2010-04-05  78.01  83.21
95  2010-04-06 113.47   1.68
96  2010-04-07  72.34   8.08
97  2010-04-08  -1.06 107.65
98  2010-04-09  27.07   9.37
99  2010-04-10   1.02  66.48
100 2010-04-11  31.06  34.76

把数据进行可视化,可以更直观地理解数据本身。


# 加载R语言类库
> library(ggplot2)
> library(scales)
> library(reshape2)

# 数据转型 
> df2<-melt(df,c('dates'))

# 画图
> g<-ggplot(data=df2,aes(x=dates,y=value,colour=variable))
> g<-g+geom_line()
> g<-g+scale_x_date(date_breaks = "1 week",date_labels='%m-%d')
> g<-g+labs(x='date',y='Price')
> g

01

上图中,X轴为时间,Y轴是价格,红色线为X的产品的价格,蓝色线为Y产品的价格。我们可以直观的看出,X,Y两个产品无任何关系。

根据配对交易的假设条件,如果两个金融产品的价差是收敛的。我们用X的产品价格减去Y产品的价格,当差值为正的时候,我们认为X的价格过高,则做空X,同时Y的价格过低,则做多Y;当差值为负的时候,我们认为X的价格过低,则做多X,同时Y的价格过高,则做空Y;当差值为0时,则价格被市场所修复,则全部平仓。

为了让差异更明显,我们定义的计算公式如下。


价差Z = X价格-Y价格
Z >  10时,做空X,做多Y ;Z<0时,平仓
Z < -10时,做多X,做空Y ;Z>0时,平仓

计算差价,然后计算交易统计。


# 计算差价
> df$diff<-df$x-df$y

# 找到差价大于10时的点
> idx<-which(df$diff>10)
> idx<-idx[-which(diff(idx)==1)-1]

# 打印差价的索引值
> idx
 [1]  4 11 15 23 25 30 34 36 38 43 48 53 55 59 61 68 76 81 83 86 88 92 95 98

接下来,我们进行模拟交易,取第一个索引值的点,在2010-01-04时做空X,做多Y。当差价小于0在2010-01-06时,进行平仓。


# 打印前20个数据
> head(df,20)
        dates      x      y    diff
1  2010-01-02  24.94  25.19   -0.25
2  2010-01-03  57.35  51.68    5.67
3  2010-01-04  16.57  13.56    3.01
4  2010-01-05 113.81  56.32   57.49
5  2010-01-06  63.18  23.82   39.36
6  2010-01-07  17.18 120.69 -103.51
7  2010-01-08  69.50  78.67   -9.17
8  2010-01-09  79.53  86.41   -6.88
9  2010-01-10  73.03  65.37    7.66
10 2010-01-11  37.78 117.29  -79.51
11 2010-01-12 110.47  24.57   85.90
12 2010-01-13  65.59  31.53   34.06
13 2010-01-14  25.15 107.29  -82.14
14 2010-01-15 -38.59  23.97  -62.56
15 2010-01-16  95.00  41.70   53.30
16 2010-01-17  48.20  34.29   13.91
17 2010-01-18  49.35  37.20   12.15
18 2010-01-19  87.75  38.84   48.91
19 2010-01-20  82.85  69.77   13.08
20 2010-01-21  73.76  42.91   30.85

# 当差价大于10时,做空X,当差价小于0时,平仓。
# 第4行做空,第6行平仓
> xprofit<- df$x[4]-df$x[6];xprofit
[1] 96.63

# 当差价大于10时,做多Y;当差价小于0时,平仓。
# 第4行做空,第6行平仓
> yprofit<- df$y[6]-df$y[4];yprofit
[1] 64.37

从交易结果来看,我们第一笔配对交易就是赚钱的。

这是为什么呢?

根据配对交易的假设条件,如果两个金融产品的价差是收敛的,通过协整性检验的方法,我们可验证数据的收敛性。那么如果数据是收敛的,他还会具备均值回归的特性,请参考文章 均值回归,逆市中的投资机会

画出X,Y的价差图,我们可以明显的看出,价差一直围绕着0上下波动,这是明显收敛的,同时符合均值回归的特性。


> plot(df$diff,type='l')

02

这就是市场的规则,通过配对交易的方法,我们找到市场无效性,从而可以赚去套利的收益。

3. 用R语言实现配对交易

看到上面的赚钱方法,也许大家会很兴奋!但是大部分市场的数据,都不会像我们的假设条件一样,轻而易举就能实现赚钱的目标。我们可以用计算机程序进行全市场的扫描发现交易机会,当然你也可以通过肉眼的方式来观察。

市场上有一些天生就具备均衡关系的金融产品,可以作为我们套利的入手对象。

  • 股票类,同行业、市值和基本面相似的个股,比如,中国银行(601988)和农业银行(601288)。
  • 基金类,以相同指数作为标的的不同基金,比如,证券B(150172),券商B(150201)。
  • 期货类,同一期货品种的不同合约,比如,铜(cu1605, cu1606)。
  • 混合类,跨市场为标的的金融产品,比如,沪深300指数,IF的期货合约

接下来,以相同品种不同合约的期货为例,我们把配对交易用在cu1605和cu1606的两个合约上,试试效果如何。由于期货是支持的T+0日内的交易的,而对于套利的操作,通常都不会持仓过夜,所以我们在尽量的短周期上进行操作,而且日内平仓。下面我将以1分钟做为交易周期。

3.1 数据准备

R语言本身提供了丰富的金融函数工具包,时间序列包zoo和xts,指标计算包TTR,可视包ggplot2等,我们会一起使用这些工具包来完成建模、计算和可视化的工作。关于zoo包和xts包的详细使用可以参考文章,R语言时间序列基础库zoo可扩展的时间序列xts

本文用到的数据,是铜的1分钟线的数据,从2016年日2月1日到2016年日2月29日,日盘的交易数据,以CSV格式保存到本地文件cu1605.csv,cu1606.csv。商品期货的日盘交易时间分为3段:09:00:00-10:14:59,10:30:00-11:29:59,13:30:00-14:59:59。当前测试,不考虑夜盘的数据。

数据格式如下:


2016-02-01 09:00:00,35870,35900,35860,35880
2016-02-01 09:01:00,35890,35890,35860,35870
2016-02-01 09:02:00,35870,35870,35860,35870
2016-02-01 09:03:00,35870,35900,35870,35900
2016-02-01 09:04:00,35900,35900,35870,35870
2016-02-01 09:05:00,35870,35880,35860,35870
2016-02-01 09:06:00,35880,35880,35860,35870

一共5列:

  • 第1列,交易时间,date,2016-02-01 09:00:00
  • 第2列,开盘价,Open,35870
  • 第3列,最高价,High,35900
  • 第4列,最低价,Low,35860
  • 第5列,收盘价,Close,35880

通过R语言加载铜的1分钟线数据,因为我们进行日内交易,所以在加载时我就进行了转换,按日期进行分组,生成R语言的list对象,同时把每日的data.frame类型对象转成XTS时间序列类型对象,方便后续的数据处理。


#加载工具包
> library(xts)
> library(TTR)

# 读取CSV数据文件
> read<-function(file){ 
+     df<-read.table(file=file,header=FALSE,sep = ",", na.strings = "NULL")  # 读文件
+     names(df)<-c("date","Open","High","Low","Close")                       # 设置列名
+     dl<-split(df,format(as.POSIXct(df$date),'%Y-%m-%d'))                   # 按日期分组
+     
+     lapply(dl,function(item){                                              # 换成xts类型数据
+         xts(item[-1],order.by = as.POSIXct(item$date))
+     })
+ }

# 加载数据
> cu1605<-read(file='cu1605.csv')
> cu1606<-read(file='cu1606.csv')

# 查看数据类型
> class(cu1605)
[1] "list"

# 查看数据的日期索引值
> names(cu1605)
 [1] "2016-02-01" "2016-02-02" "2016-02-03" "2016-02-04" "2016-02-05"
 [6] "2016-02-15" "2016-02-16" "2016-02-17" "2016-02-18" "2016-02-19"
[11] "2016-02-22" "2016-02-23" "2016-02-24" "2016-02-25" "2016-02-26"
[16] "2016-02-29"

# 查看每日的数据量
> nrow(cu1605[[1]])
[1] 223

# 查看cu1605合约的数据
> head(cu1605[['2016-02-01']])
                     Open  High   Low Close
2016-02-01 09:00:00 35870 35900 35860 35880
2016-02-01 09:01:00 35890 35890 35860 35870
2016-02-01 09:02:00 35870 35870 35860 35870
2016-02-01 09:03:00 35870 35900 35870 35900
2016-02-01 09:04:00 35900 35900 35870 35870
2016-02-01 09:05:00 35870 35880 35860 35870

把数据准备好了,我们就可以来建立模型了。

3.2 配对交易模型

以2016年02月01日为例进行交易,以1分钟线的close价格来计算cu1605和cu1606的两个合约的价差。下面我们对数据进行操作,合并2个合约在2016年02月01日的数据,并对空值进行处理,最后计算出两个合约的价差。


# 合并数据
> xdf<-merge(cu1605[['2016-02-01']]$Close,cu1606[['2016-02-01']]$Close)
> names(xdf)<-c('x1','x2')

# 用前值替换空值
> xdf<-na.locf(xdf)

# 计算价差
> xdf$diff<-xdf$x1-xdf$x2

# 打印前20行数据
> head(xdf,20)
                     x1     x2     diff
2016-02-01 09:00:00  35880  35900  -20
2016-02-01 09:01:00  35870  35920  -50
2016-02-01 09:02:00  35870  35910  -40
2016-02-01 09:03:00  35900  35940  -40
2016-02-01 09:04:00  35870  35910  -40
2016-02-01 09:05:00  35870  35920  -50
2016-02-01 09:06:00  35870  35910  -40
2016-02-01 09:07:00  35860  35910  -50
2016-02-01 09:08:00  35840  35880  -40
2016-02-01 09:09:00  35790  35840  -50
2016-02-01 09:10:00  35800  35840  -40
2016-02-01 09:11:00  35790  35830  -40
2016-02-01 09:12:00  35820  35860  -40
2016-02-01 09:13:00  35810  35850  -40
2016-02-01 09:14:00  35790  35830  -40
2016-02-01 09:15:00  35780  35830  -50
2016-02-01 09:16:00  35770  35810  -40
2016-02-01 09:17:00  35760  35820  -60
2016-02-01 09:18:00  35750  35800  -50
2016-02-01 09:19:00  35760  35810  -50

数据解释:

  • x1列,为第一腿对应cu1605合约
  • x2列,为第二腿对应cu1606合约。
  • diff列,为cu1605-cu1606

从价差的结果看,每1分钟cu1605合约都小于cu1606合约,从-110到-20价差不等,并且以-63为均值上下反复震荡。


# 计算价差范围
> range(xdf$diff)
[1] -110  -20

# 计算价差均值
> mean(xdf$diff)
[1] -63.90135

# 画出价差分布柱状图
> hist(xdf$diff,10)

画出价差分布柱状图
03

我们假设以-63为均值回归点,当差值为大于-45的时候,认为X的价格过高做空X,同时Y的价格过低做多Y;当差值小于-75的时候,我们认为X的价格过低做多X,同时Y的价格过高做空Y;当差值为-63时,价格被市场所修复,则全部平仓。以cu1605和cu1606的两个合约按照1:1持仓进行配比,1手多单对1手空单。

定义模型指标,计算价值列为diff,均值回归列为mid,最大阈值列为top,最小阈值列为bottom。


target.pair<-function(xdf){
  xdf$diff<-xdf$x1-xdf$x2   #差值
  xdf$mid<- -63             #均值回归点
  xdf$top<- -45             #最大阈值
  xdf$bottom<- -75          #最小阈值
  return(xdf)
}

完成指标的定义后,我们创建配对交易模型,并对合同数据进行回测,产生交易信号后,模拟交易输出清单,并可视化交易结果。

回测过程代码省略,产生的交易信号如下所示。


                  date    x1    x2 diff mid top bottom op
21 2016-02-01 09:00:00 35880 35900  -20 -63 -45    -75 ks
1  2016-02-01 09:25:00 35740 35810  -70 -63 -45    -75 pb
22 2016-02-01 09:40:00 35690 35730  -40 -63 -45    -75 ks
2  2016-02-01 09:47:00 35700 35770  -70 -63 -45    -75 pb
13 2016-02-01 10:00:00 35690 35770  -80 -63 -45    -75 kb
5  2016-02-01 10:01:00 35710 35760  -50 -63 -45    -75 ps
23 2016-02-01 10:02:00 35710 35750  -40 -63 -45    -75 ks
3  2016-02-01 10:07:00 35680 35750  -70 -63 -45    -75 pb
14 2016-02-01 10:37:00 35720 35800  -80 -63 -45    -75 kb
6  2016-02-01 10:42:00 35740 35790  -50 -63 -45    -75 ps
15 2016-02-01 11:20:00 35700 35780  -80 -63 -45    -75 kb
7  2016-02-01 11:21:00 35710 35750  -40 -63 -45    -75 ps
24 2016-02-01 11:21:00 35710 35750  -40 -63 -45    -75 ks
4  2016-02-01 11:23:00 35690 35760  -70 -63 -45    -75 pb
16 2016-02-01 11:29:00 35690 35770  -80 -63 -45    -75 kb
8  2016-02-01 13:36:00 35660 35720  -60 -63 -45    -75 ps
17 2016-02-01 13:45:00 35660 35740  -80 -63 -45    -75 kb
9  2016-02-01 13:46:00 35670 35730  -60 -63 -45    -75 ps
18 2016-02-01 13:52:00 35650 35730  -80 -63 -45    -75 kb
10 2016-02-01 13:53:00 35650 35710  -60 -63 -45    -75 ps
19 2016-02-01 13:56:00 35640 35720  -80 -63 -45    -75 kb
11 2016-02-01 14:49:00 35600 35660  -60 -63 -45    -75 ps
20 2016-02-01 14:52:00 35610 35700  -90 -63 -45    -75 kb
12 2016-02-01 14:58:00 35610 35690  -80 -63 -45    -75 ps

数据解释:

  • date列,为交易时间
  • x1列,为第一腿对应cu1605合约
  • x2列,为第二腿对应cu1606合约。
  • diff列,为cu1605-cu1606
  • mid列,为均值回归点
  • top列,为最大阈值
  • bottom列,为最小阈值
  • op列,为交易信号

交易信号一共有4种。

  • ks, 开仓, 做空(卖),对应反向操作为pb。
  • kb, 开仓, 做多(买),对应反向操作为ps。
  • ps, 平仓, 做空(卖),对应反向操作为kb。
  • pb,平仓, 做多(买),对应反向操作为ks。

一共出现了24个交易信号,由于我们进行的是配对交易,所以当出现ks(开仓做空)信号时,实际上会进行2笔操作,开仓做空第一腿,开仓做多第二腿。

接下来,进行模拟交易,计算出交易清单。


$x1
                       code op price pos    fee  value  margin balance     cash
2016-02-01 09:00:00  cu1605 ks 35880   1 8.9700 179400 26910.0      NA 173081.0
2016-02-01 09:25:00  cu1605 pb 35740   0 8.9350      0     0.0     700 173748.1
2016-02-01 09:40:00  cu1605 ks 35690   1 8.9225 178450 26767.5      NA 173437.7
2016-02-01 09:47:00  cu1605 pb 35700   0 8.9250      0     0.0     -50 173339.9
2016-02-01 10:00:00  cu1605 kb 35690   1 8.9225 178450 26767.5      NA 173552.0
2016-02-01 10:01:00  cu1605 ps 35710   0 8.9275      0     0.0     100 173574.2
2016-02-01 10:02:00  cu1605 ks 35710   1 8.9275 178550 26782.5      NA 173651.3
2016-02-01 10:07:00  cu1605 pb 35680   0 8.9200      0     0.0     150 173753.4
2016-02-01 10:37:00  cu1605 kb 35720   1 8.9300 178600 26790.0      NA 173758.1
2016-02-01 10:42:00  cu1605 ps 35740   0 8.9350      0     0.0     100 173780.2
2016-02-01 11:20:00  cu1605 kb 35700   1 8.9250 178500 26775.0      NA 173887.3
2016-02-01 11:21:00  cu1605 ps 35710   0 8.9275      0     0.0      50 173859.4
2016-02-01 11:21:001 cu1605 ks 35710   1 8.9275 178550 26782.5      NA 174044.1
2016-02-01 11:23:00  cu1605 pb 35690   0 8.9225      0     0.0     100 174096.2
2016-02-01 11:29:00  cu1605 kb 35690   1 8.9225 178450 26767.5      NA 174173.3
2016-02-01 13:36:00  cu1605 ps 35660   0 8.9150      0     0.0    -150 173945.5
2016-02-01 13:45:00  cu1605 kb 35660   1 8.9150 178300 26745.0      NA 174260.1
2016-02-01 13:46:00  cu1605 ps 35670   0 8.9175      0     0.0      50 174232.3
2016-02-01 13:52:00  cu1605 kb 35650   1 8.9125 178250 26737.5      NA 174331.9
2016-02-01 13:53:00  cu1605 ps 35650   0 8.9125      0     0.0       0 174254.1
2016-02-01 13:56:00  cu1605 kb 35640   1 8.9100 178200 26730.0      NA 174403.8
2016-02-01 14:49:00  cu1605 ps 35600   0 8.9000      0     0.0    -200 174125.9
2016-02-01 14:52:00  cu1605 kb 35610   1 8.9025 178050 26707.5      NA 174490.6
2016-02-01 14:58:00  cu1605 ps 35610   0 8.9025      0     0.0       0 174405.3

$x2
                       code op price pos    fee  value  margin balance     cash
2016-02-01 09:00:00  cu1606 kb 35900   1 8.9750 179500 26925.0      NA 146147.1
2016-02-01 09:25:00  cu1606 ps 35810   0 8.9525      0     0.0    -450 200214.2
2016-02-01 09:40:00  cu1606 kb 35730   1 8.9325 178650 26797.5      NA 146631.3
2016-02-01 09:47:00  cu1606 ps 35770   0 8.9425      0     0.0     200 200328.4
2016-02-01 10:00:00  cu1606 ks 35770   1 8.9425 178850 26827.5      NA 146715.6
2016-02-01 10:01:00  cu1606 pb 35760   0 8.9400      0     0.0      50 200442.7
2016-02-01 10:02:00  cu1606 kb 35750   1 8.9375 178750 26812.5      NA 146829.8
2016-02-01 10:07:00  cu1606 ps 35750   0 8.9375      0     0.0       0 200557.0
2016-02-01 10:37:00  cu1606 ks 35800   1 8.9500 179000 26850.0      NA 146899.1
2016-02-01 10:42:00  cu1606 pb 35790   0 8.9475      0     0.0      50 200671.2
2016-02-01 11:20:00  cu1606 ks 35780   1 8.9450 178900 26835.0      NA 147043.4
2016-02-01 11:21:00  cu1606 pb 35750   0 8.9375      0     0.0     150 200835.5
2016-02-01 11:21:001 cu1606 kb 35750   1 8.9375 178750 26812.5      NA 147222.6
2016-02-01 11:23:00  cu1606 ps 35760   0 8.9400      0     0.0      50 200949.8
2016-02-01 11:29:00  cu1606 ks 35770   1 8.9425 178850 26827.5      NA 147336.9
2016-02-01 13:36:00  cu1606 pb 35720   0 8.9300      0     0.0     250 201014.1
2016-02-01 13:45:00  cu1606 ks 35740   1 8.9350 178700 26805.0      NA 147446.2
2016-02-01 13:46:00  cu1606 pb 35730   0 8.9325      0     0.0      50 201078.4
2016-02-01 13:52:00  cu1606 ks 35730   1 8.9325 178650 26797.5      NA 147525.5
2016-02-01 13:53:00  cu1606 pb 35710   0 8.9275      0     0.0     100 201142.7
2016-02-01 13:56:00  cu1606 ks 35720   1 8.9300 178600 26790.0      NA 147604.8
2016-02-01 14:49:00  cu1606 pb 35660   0 8.9150      0     0.0     300 201207.0
2016-02-01 14:52:00  cu1606 ks 35700   1 8.9250 178500 26775.0      NA 147706.7
2016-02-01 14:58:00  cu1606 pb 35690   0 8.9225      0     0.0      50 201221.4

数据解释:

  • $x1部分,为第一腿的交易清单。
  • $x2部分,为第二腿的交易清单。
  • code,合约代码
  • op,交易信号
  • price,成交价格
  • pos,成交数量
  • fee,手续费
  • value,对应价值
  • margin,保证金
  • balance,平仓盈亏
  • cash,账号资金

我通过交易清单,统计交易结果。


> page  
$day     # 交易日期
[1] "2016-02-01"

$capital   # 初始资金
[1] 2e+05

$cash      # 账户余额
[1] 201221.4

$num       # 交易信号数
[1] 24

$record    # 配对交易平仓盈亏
                      x1   x2 balance
2016-02-01 09:25:00  700 -450     250
2016-02-01 09:47:00  -50  200     150
2016-02-01 10:01:00  100   50     150
2016-02-01 10:07:00  150    0     150
2016-02-01 10:42:00  100   50     150
2016-02-01 11:21:00   50  150     200
2016-02-01 11:23:00  100   50     150
2016-02-01 13:36:00 -150  250     100
2016-02-01 13:46:00   50   50     100
2016-02-01 13:53:00    0  100     100
2016-02-01 14:49:00 -200  300     100
2016-02-01 14:58:00    0   50      50

$balance   # 汇总平仓盈亏,第一腿盈亏,第二腿盈亏
[1] 1650  850  800

$fee       # 汇总手费费,第一腿手续费,第二腿手续费
[1] 429 214 215

$profit    # 账户净收益,收益率(占保证金)
[1] 1221.000    0.023

$wins      # 胜率,胜数,败数
[1]  1 12  0

最后,通过可视化输出交易信号。

04

图例解释

  • 棕色线,为价差diff
  • 紫色线,为最大阈值top
  • 红色线,为最小阈值bottom
  • 蓝色线,为均值线mid,平行于top和bottom
  • 浅蓝线,为ks开仓做空的交易
  • 绿色线,为kb开仓做多的交易

从图中看就更直观了,我们进行了12次交易,每次4笔,胜率100%。

最后,我们对2月份整个的数据进行回测。回测结果如下。


         date profit    ret balance fee winRate win fail maxProfit maxLoss avgProfit avgLoss
1  2016-02-01   1221  0.023    1650 429    1.00  12    0       250      50       138     NaN
2  2016-02-02   1077  0.020    1650 573    1.00  15    0       150       0       110     NaN
3  2016-02-03     64  0.001     100  36    1.00   1    0       100     100       100     NaN
4  2016-02-04    113  0.002     150  37    1.00   1    0       150     150       150     NaN
5  2016-02-05    926  0.017    1400 474    1.00  13    0       150     100       108     NaN 
6  2016-02-15   1191  0.022    1550 359    1.00  10    0       250     100       155     NaN 
7  2016-02-16     78  0.001     150  72    1.00   1    0       150       0       150     NaN  
8  2016-02-17    179  0.003     250  71    1.00   2    0       200      50       125     NaN 
9  2016-02-18     14  0.000      50  36    1.00   1    0        50      50        50     NaN  
10 2016-02-19    -36 -0.001       0  36     NaN   0    0         0       0       NaN     NaN   
11 2016-02-22     64  0.001     100  36    1.00   1    0       100     100       100     NaN 
12 2016-02-23    632  0.012     850 218    1.00   6    0       200     100       142     NaN 
13 2016-02-24    470  0.009     650 180    1.00   4    0       200       0       162     NaN 
14 2016-02-25    114  0.002     150  36    1.00   1    0       150     150       150     NaN 
15 2016-02-26    178  0.003     250  72    1.00   2    0       150     100       125     NaN 
16 2016-02-29    511  0.009     800 289    0.88   7    1       150     -50       121     -50 

数据解释:

  • date,交易日期
  • profit,净收益
  • ret,每日收益率
  • balance,平仓盈亏
  • fee,手续费
  • winRate,胜率
  • win,胜数
  • fail,败数
  • maxProfit,单笔最大盈利
  • maxLoss,单笔最大亏损
  • avgProfit,平均盈利
  • avgLoss,平均亏损

从结果来看,多么开心啊,几乎每天都是赚钱的!!

cu1605和cu1606两个合同是完美地具备均衡关系的两个金融产品,大家常常所说的跨期套利就是基于这个思路实现的。本文介绍的配对交易模型,是统计套利的一个基本模型,原理很简单,当大家都掌握后拼的就是交易速度了。

利用市场的无效性来获取利润,是每个套利策略都在寻找的目标。通过统计方法,我们可以发现市场的无效性,再以对冲的操作方式,规避绝大部分的市场风险,等待市场的自我修复后来赚钱利润。说起来很简单,但市场的无效性,可能会在极短时间内就被修复。

“天下武功为快不破”,通过量化的手段,让计算机来发现机会,进行交易,实现收益。一切就和谐了!!

转载请注明出处:
http://blog.fens.me/finance-pairs-trading/

打赏作者