• 粉丝日志首页

R语言跨界调用C++

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

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

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

关于作者:

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

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

rcpp

前言

使用R语言已经很多年了,对很多的R包都已经了解,唯独没有碰和C++相关的部分,这可能很大的原因和我长期使用Java的背景有关。但随着多语言的发展,跨语言应用的流行,打通各语言界限的方法也已经是成熟。让R和C++实现通信,已经变得很简单。

跟上跨语言的步伐,打开R和C++的通道,让C++来解决R性能的诟病吧。

目录

  1. Rcpp的简单介绍
  2. 5分钟上手
  3. 数据类型转换

1. Rcpp的简单介绍

Rcpp包是一个打通R语言和C++语言的通信组件包,提供了R语言和C++函数的相互调用。R语言和C++语言的数据类型通过Rcpp包进行完整的映射。

Rcpp的官方网站:https://cran.r-project.org/web/packages/Rcpp/index.html

本文做为入门教程,只是简单介绍,如何能打通R语言和C++的通信通道,并不做深入地探讨。R语言和其他语言也有类似的通信实现,R语言和JAVA的调用,请参考文章解惑rJava R与Java的高速通道;R语言和Nodejs的调用,请参考文章Nodejs与R跨平台通信

2. 5分钟上手

做为5分钟上手的教程,我们只讲例子不讲API。

本文的系统环境

  • Win10 64bit
  • R version 3.2.3 (2015-12-10)

由于Windows系统的环境下需要Rtools支持,所以要手动下载对应版本的Rtoosl包,下载地址。我的R语言版本是3.2.3,所以我需要安装Rtools33.exe。安装EXE程序就不多说了,双击完成即可。

rtools

下载Rcpp的程序包,进行安装,一行代码搞定。


> install.packages("Rcpp")
trying URL 'https://mirrors.tuna.tsinghua.edu.cn/CRAN/bin/windows/contrib/3.2/Rcpp_0.12.6.zip'
Content type 'application/zip' length 3221864 bytes (3.1 MB)
downloaded 3.1 MB

package ‘Rcpp’ successfully unpacked and MD5 sums checked
Warning in install.packages :
  cannot remove prior installation of package ‘Rcpp’

The downloaded binary packages are in
	C:\Users\tinkpad\AppData\Local\Temp\RtmpKkg8zo\downloaded_packages

2.1 从hello world开始

从一个简单程序hello world开始吧,让R语言程序调用C++中的hello()函数。我用需要新建2个文件,放在同一个目录中。

  • demo.cpp,C++程序的源文件
  • demo.r,R程序源文件

首先,编辑demo.cpp,定义hello()函数。


~ notepad demo.cpp

#include <Rcpp.h>
#include <string>  

using namespace std;
using namespace Rcpp;

//[[Rcpp::export]]
string hello(string name) {
  cout << "hello " << name << endl;  
  return name;
}

/*** R
hello('world')
hello('Conan')
*/

上面Rcpp的代码,我们可以从3部分来看。

  • #include和using部分: 为包引用和命名空间的声明。<Rcpp.h>和namespace Rcpp是必要要加载的,另外由于使用了string的类型作为参数和返回值,所以需要<string>和namespace std。
  • 功能函数部分: 们定义了一个 hello(string name) 函数,有一个参数是string类型,返回值也为string类型。需要强调的是,对R开放的函数必须增加 //[[Rcpp::export]] 的注释声明
  • 代码执行: 用/*** R 和 */ 包含的部分,为R语言的代码,会默认被执行。

编辑demo.r,用来调用demo.cpp的hello()函数。


~ notepad demo.r

library(Rcpp)

sourceCpp(file='demo.cpp')
hello('R')

执行R语言的代码


# 加载Rcpp包
> library(Rcpp)

# 编译和加载demo.cpp文件
> sourceCpp(file='demo.cpp')

# 执行封装在demo.cpp中的R代码
> hello('world')
hello world
[1] "world"

> hello('Conan')
hello Conan
[1] "Conan"

# 执行hello函数
> hello('R')
hello [1]R
 "R"

一个非常简单的helloworld程序,就这样子完成了。

2.2 R和Rcpp的混写代码

上面2行代码,就完成了R对C++程序的调用,sourceCpp()函数真是强大。其实,sourceCpp()函数还提供了一种代码混写的方法,就是在R的代码中,直接嵌入C++代码。


sourceCpp(code='
  #include >Rcpp.h<
  #include >string<
  
  using namespace std;
  using namespace Rcpp;
  
  //[[Rcpp::export]]
  string hello(string name) {
    cout << "hello " << name << endl;  
    return name;
  }
')
hello('R2')

运行代码


> sourceCpp(code='
+   #include >Rcpp.h<
+   #include >string<
+   
+   using namespace std;
+   using namespace Rcpp;
+   
+   //[[Rcpp::export]]
+   string hello(string name) {
+     cout << "hello " << name << endl;  
+     return name;
+   }
+ ')

> hello('R2')
hello R2
[1] "R2"

这种多语言混写的语法虽然不太推荐,但对于这只有几行代码来说,还是很方便的。

2.2 用RStudioIDE生成cpp文件

如果你使用的RStudio IDE,开发起来将会非常方便,可以直接新建C++程序,生成一段标准的代码模板。

rstudio-cpp

生成的代码模板如下


#include <Rcpp.h>
using namespace Rcpp;

// This is a simple example of exporting a C++ function to R. You can
// source this function into an R session using the Rcpp::sourceCpp 
// function (or via the Source button on the editor toolbar). Learn
// more about Rcpp at:
//
//   http://www.rcpp.org/
//   http://adv-r.had.co.nz/Rcpp.html
//   http://gallery.rcpp.org/
//

// [[Rcpp::export]]
NumericVector timesTwo(NumericVector x) {
  return x * 2;
}

// You can include R code blocks in C++ files processed with sourceCpp
// (useful for testing and development). The R code will be automatically 
// run after the compilation.
//

/*** R
timesTwo(42)
*/

通过RStudio可以快速生成一段标准的代码模板,改改马上就能用了。

3. 数据类型转换

上面的例子中,我们测试了字符串类型的调用。R语言有多种的数据类型,我接下来都测试一下!

3.1 基本类型

基本类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
char char_type(char x){
  return x;
}

// [[Rcpp::export]]
int int_type(int x){
  return x;
}

// [[Rcpp::export]]
double double_type(double x){
  return x;
}

// [[Rcpp::export]]
bool bool_type(bool x){
  return x;
}

// [[Rcpp::export]]
void void_return_type(){
  Rprintf( "return void" );
}

执行R语言调用


# char类型
> a1<-char_type('a')
> a1;class(a1)         # 默认对应R的character类型
[1] "a"
[1] "character"
> char_type('bbii')    # 只处理字符串的第一个字节
[1] "b"

# int类型 
> a2<-int_type(111)
> a2;class(a2)         # 默认对应R的integer类型
[1] 111
[1] "integer" 
> int_type(111.1)      # 直接去掉小数位
[1] 111

# double类型 
> a3<-double_type(111.1)
> a3;class(a3)         # 默认对应R的numeric类型
[1] 111.1
[1] "numeric"
> double_type(111)
[1] 111

# boolean类型 
> a4<-bool_type(TRUE)
> a4;class(a4)        # 默认对应R的logical类型
[1] TRUE
[1] "logical"
> bool_type(0)        # 0为FALSE
[1] FALSE
> bool_type(1)        # 非0为TRUE
[1] TRUE

# 无参数无返回值 的函数
> a5<-void_return_type()
return void
> a5;class(a5)         # 返回值为NULL
NULL
[1] "NULL"

3.2 向量类型

向量类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
CharacterVector CharacterVector_type(CharacterVector x){
  return x;
}

// [[Rcpp::export]]
StringVector StringVector_type(StringVector x){
  return x;
}

// [[Rcpp::export]]
NumericVector NumericVector_type(NumericVector x) {
  return x;
}

// [[Rcpp::export]]
IntegerVector IntegerVector_type(IntegerVector x) {
  return x;
}

// [[Rcpp::export]]
DoubleVector DoubleVector_type(DoubleVector x){
  return x;
}

// [[Rcpp::export]]
LogicalVector LogicalVector_type(LogicalVector x){
  return x;
}

// [[Rcpp::export]]
DateVector DateVector_type(DateVector x){
  return x;
}

// [[Rcpp::export]]
DatetimeVector DatetimeVector_type(DatetimeVector x){
  return x;
}

执行R语言调用


# Character向量
> a6<-CharacterVector_type(c('abc','12345'))     
> a6;class(a6)                                    # 默认对应R的character类型
[1] "abc"   "12345"
[1] "character" 
> CharacterVector_type(c('abc',123.5, NA, TRUE))  # NA不处理
[1] "abc"   "123.5" NA      "TRUE" 

# String向量,完全同Character向量
> a7<-StringVector_type(c('abc','12345'))
> a7;class(a7)                                    # 默认对应R的character类型
[1] "abc"   "12345"
[1] "character"
> StringVector_type(c('abc',123.5, NA, TRUE))
[1] "abc"   "123.5" NA      "TRUE" 

# Numeric向量
> a8<-NumericVector_type(rnorm(5))
> a8;class(a8)                                    # 默认对应R的numeric类型
[1] -0.2813472 -0.2235722 -0.6958443 -1.5322172  0.5004307
[1] "numeric"
> NumericVector_type(c(rnorm(5),NA,TRUE))         # NA不处理,TRUE为1
[1]  0.1700925  0.5169612 -0.3622637  1.0763204 -0.5729958
[6]         NA  1.0000000

# Integer向量
> a9<-IntegerVector_type(c(11,9.9,1.2))           # 直接去掉小数位
> a9;class(a9)                                    # 默认对应R的integer类型
[1] 11  9  1
[1] "integer"
> IntegerVector_type(c(11,9.9,1.2,NA,TRUE))       # NA不处理,TRUE为1
[1] 11  9  1 NA  1

# Double向量,同Numeric向量
> a10<-DoubleVector_type(rnorm(5))
> a10;class(a10)                                  # 默认对应R的numeric类型
[1]  0.9400947 -0.8976913  0.2744319 -1.5278219  1.2010569
[1] "numeric"
> DoubleVector_type(c(rnorm(5),NA,TRUE))          # NA不处理,TRUE为1
[1]  2.0657148  0.2810003  2.1080900 -1.2783693  0.2198551
[6]         NA  1.0000000

# Logical向量 
> a11<-LogicalVector_type(c(TRUE,FALSE))
> a11;class(a11)                                  # 默认对应R的logical类型
[1]  TRUE FALSE
[1] "logical"
> LogicalVector_type(c(TRUE,FALSE,TRUE,0,-1,NA))  # NA不处理,0为FALSE, 非0为TRUE
[1]  TRUE FALSE  TRUE FALSE  TRUE    NA

 # Date向量 
> a12<-DateVector_type(c(Sys.Date(),as.Date('2016-10-10')))
> a12;class(a12)                                  # 默认对应R的Date类型
[1] "2016-08-01" "2016-10-10"
[1] "Date"
> DateVector_type(c(Sys.Date(),as.Date('2016-10-10'),NA,TRUE,FALSE))   # NA不处理,TRUE为1970-01-02, FALSE为1970-01-01
[1] "2016-08-01" "2016-10-10" NA           "1970-01-02"
[5] "1970-01-01"

 # Datetime向量 
> a13<-DatetimeVector_type(c(Sys.time(),as.POSIXct('2016-10-10')))
> a13;class(a13)                                  # 默认对应R的POSIXct类型
[1] "2016-08-01 20:05:25 CST" "2016-10-10 00:00:00 CST"
[1] "POSIXct" "POSIXt" 
> DatetimeVector_type(c(Sys.time(),as.POSIXct('2016-10-10'),NA,TRUE,FALSE))  # NA不处理
[1] "2016-08-01 20:05:25 CST" "2016-10-10 00:00:00 CST"
[3] NA                        "1970-01-01 08:00:01 CST"
[5] "1970-01-01 08:00:00 CST"

3.3 矩阵类型

矩阵类型,C++对应R语言的默认映射关系。C++的代码部分,如下所示:


// [[Rcpp::export]]
CharacterMatrix CharacterMatrix_type(CharacterMatrix x){
  return x;
}

// [[Rcpp::export]]
StringMatrix StringMatrix_type(StringMatrix x){
  return x;
}

// [[Rcpp::export]]
NumericMatrix NumericMatrix_type(NumericMatrix x){
  return x;
}

// [[Rcpp::export]]
IntegerMatrix IntegerMatrix_type(IntegerMatrix x){
  return x;
}

// [[Rcpp::export]]
LogicalMatrix LogicalMatrix_type(LogicalMatrix x){
  return x;
}

// [[Rcpp::export]]
ListMatrix ListMatrix_type(ListMatrix x){
  return x;
}

执行R语言调用


# Character矩阵
> a14<-CharacterMatrix_type(matrix(LETTERS[1:20],ncol=4))
> a14;class(a14)
     [,1] [,2] [,3] [,4]
[1,] "A"  "F"  "K"  "P" 
[2,] "B"  "G"  "L"  "Q" 
[3,] "C"  "H"  "M"  "R" 
[4,] "D"  "I"  "N"  "S" 
[5,] "E"  "J"  "O"  "T" 
[1] "matrix"                        

# String矩阵,同Character矩阵
> a15<-StringMatrix_type(matrix(LETTERS[1:20],ncol=4))
> a15;class(a15)
     [,1] [,2] [,3] [,4]
[1,] "A"  "F"  "K"  "P" 
[2,] "B"  "G"  "L"  "Q" 
[3,] "C"  "H"  "M"  "R" 
[4,] "D"  "I"  "N"  "S" 
[5,] "E"  "J"  "O"  "T" 
[1] "matrix"

# Numeric矩阵
> a16<-NumericMatrix_type(matrix(rnorm(20),ncol=4))
> a16;class(a16)
           [,1]       [,2]       [,3]       [,4]
[1,]  1.2315498  2.3234269  0.5974143  0.9072356
[2,]  0.3484811  0.3814024 -0.2018324  0.8717205
[3,] -0.2025285  2.1076947 -0.3433948  1.1523710
[4,] -1.4948252 -0.7724951 -0.7681800 -0.5406494
[5,]  0.4815904  1.4930873 -1.1444258  0.2537099
[1] "matrix"

# Integer矩阵 
> a17<-IntegerMatrix_type(matrix(seq(1,10,length.out = 20),ncol=4))
> a17;class(a17)
     [,1] [,2] [,3] [,4]
[1,]    1    3    5    8
[2,]    1    3    6    8
[3,]    1    4    6    9
[4,]    2    4    7    9
[5,]    2    5    7   10
[1] "matrix"

# Logical矩阵 
> a18<-LogicalMatrix_type(matrix(c(rep(TRUE,5),rep(FALSE,5),rnorm(10)),ncol=4))
> a18;class(a18)
     [,1]  [,2] [,3] [,4]
[1,] TRUE FALSE TRUE TRUE
[2,] TRUE FALSE TRUE TRUE
[3,] TRUE FALSE TRUE TRUE
[4,] TRUE FALSE TRUE TRUE
[5,] TRUE FALSE TRUE TRUE
[1] "matrix"

# List矩阵,支持多类型的矩阵
> a19<-ListMatrix_type(matrix(rep(list(a=1,b='2',c=NA,d=TRUE),10),ncol=5))
> a19;class(a19)
     [,1] [,2] [,3] [,4] [,5]
[1,] 1    1    1    1    1   
[2,] "2"  "2"  "2"  "2"  "2" 
[3,] NA   NA   NA   NA   NA  
[4,] TRUE TRUE TRUE TRUE TRUE
[5,] 1    1    1    1    1   
[6,] "2"  "2"  "2"  "2"  "2" 
[7,] NA   NA   NA   NA   NA  
[8,] TRUE TRUE TRUE TRUE TRUE
[1] "matrix"

3.4 其他数据类型

其他数据类型包括了,R语言特有的数据类型数据框(data.frame),环境空间(Environment)S3,S4,RC等的对象类型。


// [[Rcpp::export]]
Date Date_type(Date x){
  return x;
}

// [[Rcpp::export]]
Datetime Datetime_type(Datetime x){
  return x;
}

// [[Rcpp::export]]
S4 S4_type(S4 x){
  return x;
}

// [[Rcpp::export]]
RObject RObject_type(RObject x){
  return x;
}

// [[Rcpp::export]]
SEXP SEXP_type(SEXP x){
  return x;
}

// [[Rcpp::export]]
Environment Environment_type(Environment x){
  return x;
}

执行R语言调用


# data.frame类型
> a19<-DataFrame_type(data.frame(a=rnorm(3),b=1:3))
> a19;class(a19)
           a b
1 -1.8844994 1
2  0.6053935 2
3 -0.7693985 3
[1] "data.frame"

# list类型 
> a20<-List_type(list(a=1,b='2',c=NA,d=TRUE))
> a20;class(a20)
$a
[1] 1
$b
[1] "2"
$c
[1] NA
$d
[1] TRUE
[1] "list"

# Date类型
> a21<-Date_type(Sys.Date())
> a21;class(a21)
[1] "2016-08-01"
[1] "Date"
> Date_type(Sys.time())                # 不能正确处理POSIXct类型的数据
[1] "4026842-05-26"

# POSIXct类型 
> a22<-Datetime_type(Sys.time())
> a22;class(a22)
[1] "2016-08-01 20:27:37 CST"
[1] "POSIXct" "POSIXt" 
> Datetime_type(Sys.Date())            # 不能正确处理Date类型的数据
[1] "1970-01-01 12:43:34 CST"

# S3面向对象类型,对应S4的类型定义
> setClass("Person",slots=list(name="character",age="numeric"))
> s4<-new("Person",name="F",age=44)
> a23<-S4_type(s4)
> a23;class(a23)
An object of class "Person"
Slot "name":
[1] "F"
Slot "age":
[1] 44
[1] "Person"
attr(,"package")
[1] ".GlobalEnv"

# S3面向对象类型 ,没有对应的类型,通过RObject来传值
> s3<-structure(2, class = "foo")
> a24<-RObject_type(s3)
> a24;class(a24)
[1] 2
attr(,"class")
[1] "foo"
[1] "foo"

# RObject也可以处理S4对象
> a25<-RObject_type(s4)
> a25;class(a25)
An object of class "Person"
Slot "name":
[1] "F"
Slot "age":
[1] 44
[1] "Person"
attr(,"package")
[1] ".GlobalEnv"

# RObject也可以处理RC对象 
> User<-setRefClass("User",fields=list(name="character"))
> rc<-User$new(name="u1")
> a26<-RObject_type(rc)
> a26;class(a26)
Reference class object of class "User"
Field "name":
[1] "u1"
[1] "User"
attr(,"package")
[1] ".GlobalEnv"

# RObject也可以处理function类型
> a27<-RObject_type(function(x) x+2)
> a27;class(a27)
function(x) x+2
[1] "function"

# environment类型
> a28<-Environment_type(new.env())
> a28;class(a28)
<environment: 0x0000000015350a80>
[1] "environment"

# SEXP为任意类型,通过具体调用时再进行类型判断
> SEXP_type('fdafdaa')
[1] "fdafdaa"

> SEXP_type(rc)
Reference class object of class "User"
Field "name":
[1] "u1"

> SEXP_type(data.frame(a=rnorm(3),b=1:3))
           a b
1 -0.5396140 1
2  0.1694799 2
3 -1.8818596 3

> SEXP_type(function(x) x+2)
function(x) x+2

最后总结一下,R和Rcpp中类型对应的关系。

C++类型 R类型
char character
int integer
double numeric
bool logical
Rcpp::Date Date
Rcpp::Datetime POSIXct
Rcpp::CharacterVector character
Rcpp::StringVector character
Rcpp::NumericVector numeric
Rcpp::IntegerVector integer
Rcpp::DoubleVector numeric
Rcpp::LogicalVector logical
Rcpp::DateVector Date
Rcpp::DatetimeVector POSIXct
Rcpp::CharacterMatrix matrix
Rcpp::StringMatrix matrix
Rcpp::NumericMatrix matrix
Rcpp::IntegerMatrix matrix
Rcpp::LogicalMatrix matrix
Rcpp::ListMatrix matrix
Rcpp::DataFrame data.frame
Rcpp::List list
Rcpp::S4 S4
Rcpp::Environment environment
Rcpp::RObject 任意类型
Rcpp::SEXP 任意类型

本文简单地介绍了通过R语言Rcpp包调用C++程序的一种方法,调用的关键点就在于数据类型的匹配,而从保证R语言和C++之间的数据传输。从上面测试来看,R语言中的所有数据类型,都可以通过Rcpp包进行映射到C++的程序中。接下来,我们就可以根据自己的需求,把一些更关注的性能的程序放到C++中来实现,从而提高计算效率。

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

打赏作者

R语言解读多元线性回归模型

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

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

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

关于作者:

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

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

reg-multi-liner

前言

本文接上一篇R语言解读一元线性回归模型。在许多生活和工作的实际问题中,影响因变量的因素可能不止一个,比如对于知识水平越高的人,收入水平也越高,这样的一个结论。这其中可能包括了因为更好的家庭条件,所以有了更好的教育;因为在一线城市发展,所以有了更好的工作机会;所处的行业赶上了大的经济上行周期等。要想解读这些规律,是复杂的、多维度的,多元回归分析方法更适合解读生活的规律。

由于本文为非统计的专业文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍多元线性回归的知识,同时配合R语言的实现。

目录

  1. 多元线性回归介绍
  2. 元线性回归建模
  3. 模型优化
  4. 案例:黑色系期货日K线数据验证

1. 多元线性回归介绍

对比一元线性回归,多元线性回归是用来确定2个或2个以上变量间关系的统计分析方法。多元线性回归的基本的分析方法与一元线性回归方法是类似的,我们首先需要对选取多元数据集并定义数学模型,然后进行参数估计,对估计出来的参数进行显著性检验,残差分析,异常点检测,最后确定回归方程进行模型预测。

由于多元回归方程有多个自变量,区别于一元回归方程,有一项很重要的操作就是自变量的优化,挑选出相关性最显著的自变量,同时去除不显著的自变量。在R语言中,有很方便地用于优化函数,可以很好的帮助我们来改进回归模型。

下面就开始多元线性回归的建模过程。

2. 多元线性回归建模

做过商品期货研究的人,都知道黑色系品种是具有产业链上下游的关系。铁矿石是炼钢的原材料,焦煤和焦炭是炼钢的能源资源,热卷即热轧卷板是以板坯为原料经加热后制成的钢板,螺纹钢是表面带肋的钢筋。

由于有产业链的关系,假设我们想要预测螺纹钢的价格,那么影响螺纹钢价格的因素可以会涉及到原材料,能源资源和同类材料等。比如,铁矿石价格如果上涨,螺纹钢就应该要涨价了。

2.1 数据集和数学模型

先从数据开始介绍,这次的数据集,我选择的期货黑色系的品种的商品期货,包括了大连期货交易所的 焦煤(JM),焦炭(J),铁矿石(I),上海期货交易所的 螺纹钢(RU) 和 热卷(HC)。

数据集为2016年3月15日,当日白天开盘的交易数据,为黑色系的5个期货合约的分钟线的价格数据。


# 数据集已存在df变量中
> head(df,20)
                       x1    x2    x3   x4    y
2016-03-15 09:01:00 754.5 616.5 426.5 2215 2055
2016-03-15 09:02:00 752.5 614.5 423.5 2206 2048
2016-03-15 09:03:00 753.0 614.0 423.0 2199 2044
2016-03-15 09:04:00 752.5 613.0 422.5 2197 2040
2016-03-15 09:05:00 753.0 615.5 424.0 2198 2043
2016-03-15 09:06:00 752.5 614.5 422.0 2195 2040
2016-03-15 09:07:00 752.0 614.0 421.5 2193 2036
2016-03-15 09:08:00 753.0 615.0 422.5 2197 2043
2016-03-15 09:09:00 754.0 615.5 422.5 2197 2041
2016-03-15 09:10:00 754.5 615.5 423.0 2200 2044
2016-03-15 09:11:00 757.0 616.5 423.0 2201 2045
2016-03-15 09:12:00 756.0 615.5 423.0 2200 2044
2016-03-15 09:13:00 755.5 615.0 423.0 2197 2042
2016-03-15 09:14:00 755.5 615.0 423.0 2196 2042
2016-03-15 09:15:00 756.0 616.0 423.5 2200 2045
2016-03-15 09:16:00 757.5 616.0 424.0 2205 2052
2016-03-15 09:17:00 758.5 618.0 424.0 2204 2051
2016-03-15 09:18:00 759.5 618.5 424.0 2205 2053
2016-03-15 09:19:00 759.5 617.5 424.5 2206 2053
2016-03-15 09:20:00 758.5 617.5 423.5 2201 2050

数据集包括有6列:

  • 索引, 为时间
  • x1, 为焦炭(j1605)合约的1分钟线的报价数据
  • x2, 为焦煤(jm1605)合约的1分钟线的报价数据
  • x3, 为铁矿石(i1605)合约的1分钟线的报价数据
  • x4, 为热卷(hc1605)合约的1分钟线的报价数据
  • y, 为螺纹钢(rb1605)合约的1分钟线的报价数据

假设螺纹钢的价格与其他4个商品的价格有线性关系,那么我们建立以螺纹钢为因变量,以焦煤、焦炭、铁矿石和热卷的为自变量的多元线性回归模型。用公式表示为:

y = a + b * x1 + c * x2 + d * x3 + e * x4 + ε
  • y,为因变量,螺纹钢
  • x1,为自变量,焦煤
  • x2,为自变量,焦炭
  • x3,为自变量,铁矿石
  • x4,为自变量,热卷
  • a,为截距
  • b,c,d,e,为自变量系数
  • ε, 为残差,是其他一切不确定因素影响的总和,其值不可观测。假定ε服从正态分布N(0,σ^2)。

通过对多元线性回归模型的数学定义,接下来让我们利用数据集做多元回归模型的参数估计。

2.2. 回归参数估计

上面公式中,回归参数 a, b, c, d,e都是我们不知道的,参数估计就是通过数据来估计出这些参数,从而确定自变量和因变量之前的关系。我们的目标是要计算出一条直线,使直线上每个点的Y值和实际数据的Y值之差的平方和最小,即(Y1实际-Y1预测)^2+(Y2实际-Y2预测)^2+ …… +(Yn实际-Yn预测)^2 的值最小。参数估计时,我们只考虑Y随X自变量的线性变化的部分,而残差ε是不可观测的,参数估计法并不需要考虑残差。

类似于一元线性回归,我们用R语言来实现对数据的回归模型的参数估计,用lm()函数来实现多元线性回归的建模过程。


# 建立多元线性回归模型
> lm1<-lm(y~x1+x2+x3+x4,data=df)

# 打印参数估计的结果
> lm1

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

这样我们就得到了y和x关系的方程。

y = 212.8780 + 0.8542 * x1 + 0.6672 * x2 - 0.6674 * x3 + 0.4821 * x4

2.3. 回归方程的显著性检验

参考一元线性回归的显著性检验,多元线性回归的显著性检验,同样是需要经过 T检验,F检验,和R^2(R平方)相关系统检验。在R语言中这三种检验的方法都已被实现,我们只需要把结果解读,我们可以summary()函数来提取模型的计算结果。


> summary(lm1)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9648 -1.3241 -0.0319  1.2403  5.4194 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 212.87796   58.26788   3.653 0.000323 ***
x1            0.85423    0.10958   7.795 2.50e-13 ***
x2            0.66724    0.12938   5.157 5.57e-07 ***
x3           -0.66741    0.15421  -4.328 2.28e-05 ***
x4            0.48214    0.01959  24.609  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.028 on 221 degrees of freedom
Multiple R-squared:  0.9725,	Adjusted R-squared:  0.972 
F-statistic:  1956 on 4 and 221 DF,  p-value: < 2.2e-16
  • T检验:所自变量都是非常显著***
  • F检验:同样是非常显著,p-value < 2.2e-16
  • 调整后的R^2:相关性非常强为0.972

最后,我们通过的回归参数的检验与回归方程的检验,得到最后多元线性回归方程为:


y = 212.87796 + 0.85423 * x1 + 0.66724 * x2 - 0.66741 * x3 + 0.48214 * x4

即

螺纹钢 = 212.87796 + 0.85423 * 焦炭 + 0.66724 * 焦煤 - 0.66741 * 铁矿石 + 0.48214 * 热卷

2.4 残差分析和异常点检测

在得到的回归模型进行显著性检验后,还要在做残差分析(预测值和实际值之间的差),检验模型的正确性,残差必须服从正态分布N(0,σ^2)。直接用plot()函数生成4种用于模型诊断的图形,进行直观地分析。


> par(mfrow=c(2,2))
> plot(lm1)

m01

  • 残差和拟合值(左上),残差和拟合值之间数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 残差QQ图(右上),数据点按对角直线排列,趋于一条直线,并被对角直接穿过,直观上符合正态分布。
  • 标准化残差平方根和拟合值(左下),数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征。
  • 标准化残差和杠杆值(右下),没有出现红色的等高线,则说明数据中没有特别影响回归结果的异常点。

结论,没有明显的异常点,残差符合假设条件。

2.5. 模型预测

我们得到了多元线性回归方程的公式,就可以对数据进行预测了。我们可以用R语言的predict()函数来计算预测值y0和相应的预测区间,并把实际值和预测值一起可视化化展示。


> par(mfrow=c(1,1))  #设置画面布局

# 预测计算
> dfp<-predict(lm1,interval="prediction")

# 打印预测时
> head(dfp,10)
                fit      lwr      upr
2014-03-21 3160.526 3046.425 3274.626
2014-03-24 3193.253 3078.868 3307.637
2014-03-25 3240.389 3126.171 3354.607
2014-03-26 3228.565 3114.420 3342.710
2014-03-27 3222.528 3108.342 3336.713
2014-03-28 3262.399 3148.132 3376.666
2014-03-31 3291.996 3177.648 3406.344
2014-04-01 3305.870 3191.447 3420.294
2014-04-02 3275.370 3161.018 3389.723
2014-04-03 3297.358 3182.960 3411.755

# 合并数据
> mdf<-merge(df$y,dfp)	 

# 画图
> draw(mdf)

m02

图例说明

  • y, 实际价格,红色线
  • fit, 预测价格,绿色线
  • lwr,预测最低价,蓝色线
  • upr,预测最高价,紫色线

从图中看出,实际价格y和预测价格fit,在大多数的时候都是很贴近的。我们的一个模型就训练好了!

3. 模型优化

上文中,我们已经很顺利的发现了一个非常不错的模型。如果要进行模型优化,可以用R语言中update()函数进行模型的调整。我们首先检查一下每个自变量x1,x2,x3,x4和因变量y之间的关系。

pairs(as.data.frame(df))

m03

从图中,我们可以发现x2与Y的关系,可能是最偏离线性的。那么,我们尝试对多元线性回归模型进行调整,从原模型中去掉x2变量。



# 模型调整
> lm2<-update(lm1, .~. -x2)

> summary(lm2)

Call:
lm(formula = y ~ x1 + x3 + x4, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0039 -1.3842  0.0177  1.3513  4.8028 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 462.47104   34.26636   13.50  < 2e-16 ***
x1            1.08728    0.10543   10.31  < 2e-16 ***
x3           -0.40788    0.15394   -2.65  0.00864 ** 
x4            0.42582    0.01718   24.79  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.142 on 222 degrees of freedom
Multiple R-squared:  0.9692,	Adjusted R-squared:  0.9688 
F-statistic:  2330 on 3 and 222 DF,  p-value: < 2.2e-16

当把自变量x2去掉后,自变量x3的T检验反而变大了,同时Adjusted R-squared变小了,所以我们这次调整是有问题的。

如果通过生产和原材料的内在逻辑分析,焦煤与焦炭属于上下游关系。焦煤是生产焦炭的一种原材料,焦炭是焦煤与其他炼焦煤经过配煤焦化形成的产品,一般生产 1 吨焦炭需要1.33 吨炼焦煤,其中焦煤至少占 30% 。

我们把焦煤 和 焦炭的关系改变一下,增加x1*x2的关系匹配到模型,看看效果。


# 模型调整
> lm3<-update(lm1, .~. + x1*x2)
> summary(lm3)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.8110 -1.3501 -0.0595  1.2019  5.3884 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7160.32231 7814.50048   0.916    0.361    
x1            -8.45530   10.47167  -0.807    0.420    
x2           -10.58406   12.65579  -0.836    0.404    
x3            -0.64344    0.15662  -4.108 5.63e-05 ***
x4             0.48363    0.01967  24.584  < 2e-16 ***
x1:x2          0.01505    0.01693   0.889    0.375    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.029 on 220 degrees of freedom
Multiple R-squared:  0.9726,	Adjusted R-squared:  0.972 
F-statistic:  1563 on 5 and 220 DF,  p-value: < 2.2e-16

从结果中发现,增加了x1*x2列后,原来的x1,x2和Intercept的T检验都不显著。继续调整模型,从模型中去掉x1,x2两个自变量。


# 模型调整
> lm4<-update(lm3, .~. -x1-x2)
> summary(lm4)

Call:
lm(formula = y ~ x3 + x4 + x1:x2, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.9027 -1.2516 -0.0167  1.2748  5.8683 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.950e+02  1.609e+01  43.183  < 2e-16 ***
x3          -6.284e-01  1.530e-01  -4.108 5.61e-05 ***
x4           4.959e-01  1.785e-02  27.783  < 2e-16 ***
x1:x2        1.133e-03  9.524e-05  11.897  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.035 on 222 degrees of freedom
Multiple R-squared:  0.9722,	Adjusted R-squared:  0.9718 
F-statistic:  2588 on 3 and 222 DF,  p-value: < 2.2e-16

从调整后的结果来看,效果还不错。不过,也并没有比最初的模型有所提高。

对于模型调整的过程,如果我们手动调整测试时,一般都会基于业务知识来操作。如果是按照数据指标来计算,我们可以用R语言中提供的逐步回归的优化方法,通过AIC指标来判断是否需要参数优化。


#对lm1模型做逐步回归
> step(lm1)
Start:  AIC=324.51
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

通过计算AIC指标,lm1的模型AIC最小时为324.51,每次去掉一个自变量都会让AIC的值变大,所以我们还是不调整比较好。

对刚才的lm3模型做逐步回归的模型调整。


#对lm3模型做逐步回归
> step(lm3)
Start:  AIC=325.7               #当前AIC
y ~ x1 + x2 + x3 + x4 + x1:x2

        Df Sum of Sq    RSS    AIC
- x1:x2  1      3.25  908.8 324.51
                905.6 325.70
- x3     1     69.47  975.1 340.41
- x4     1   2487.86 3393.5 622.25

Step:  AIC=324.51               #去掉x1*x2项的AIC
y ~ x1 + x2 + x3 + x4

       Df Sum of Sq    RSS    AIC
               908.8 324.51
- x3    1     77.03  985.9 340.90
- x2    1    109.37 1018.2 348.19
- x1    1    249.90 1158.8 377.41
- x4    1   2490.56 3399.4 620.65

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Coefficients:
(Intercept)           x1           x2           x3           x4  
   212.8780       0.8542       0.6672      -0.6674       0.4821  

通过AIC的判断,去掉X1*X2项后AIC最小,最后的检验结果告诉我们,还是原初的模型是最好的。

4. 案例:黑色系期货日K线数据验证

最后,我们用上面5个期货合约的日K线数据测试一下,找到多元回归关系。


> lm9<-lm(y~x1+x2+x3+x4,data=df)  # 日K线数据
> summary(lm9)

Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-173.338  -37.470    3.465   32.158  178.982 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 386.33482   31.07729  12.431  < 2e-16 ***
x1            0.75871    0.07554  10.045  < 2e-16 ***
x2           -0.62907    0.14715  -4.275 2.24e-05 ***
x3            1.16070    0.05224  22.219  < 2e-16 ***
x4            0.46461    0.02168  21.427  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 57.78 on 565 degrees of freedom
Multiple R-squared:  0.9844,	Adjusted R-squared:  0.9843 
F-statistic:  8906 on 4 and 565 DF,  p-value: < 2.2e-16

数据集的基本统计信息。


> summary(df)
     Index                           x1               x2       
 Min.   :2014-03-21 00:00:00   Min.   : 606.5   Min.   :494.0  
 1st Qu.:2014-10-21 06:00:00   1st Qu.: 803.5   1st Qu.:613.1  
 Median :2015-05-20 12:00:00   Median : 939.0   Median :705.8  
 Mean   :2015-05-21 08:02:31   Mean   : 936.1   Mean   :695.3  
 3rd Qu.:2015-12-16 18:00:00   3rd Qu.:1075.0   3rd Qu.:773.0  
 Max.   :2016-07-25 00:00:00   Max.   :1280.0   Max.   :898.0  

       x3              x4             y       
 Min.   :284.0   Min.   :1691   Min.   :1626  
 1st Qu.:374.1   1st Qu.:2084   1st Qu.:2012  
 Median :434.0   Median :2503   Median :2378  
 Mean   :476.5   Mean   :2545   Mean   :2395  
 3rd Qu.:545.8   3rd Qu.:2916   3rd Qu.:2592  
 Max.   :825.0   Max.   :3480   Max.   :3414  

m04

对于日K线数据,黑色系的5个品种,同样具有非常强的相关关系,那么我们就可以把这个结论应用到实际的交易中了。

本文通过多元回归的统计分析方法,介绍多元回归在金融市场的基本应用。我们通过建立因变量和多个自变量的模型,从而发现生活中更复杂的规律,并建立有效的验证指标。让我们我们的技术优势,去金融市场抢钱吧。

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

打赏作者

R语言解读一元线性回归模型

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

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

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

关于作者:

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

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

reg-liner

前言

在我们的日常生活中,存在大量的具有相关性的事件,比如大气压和海拔高度,海拔越高大气压强越小;人的身高和体重,普遍来看越高的人体重也越重。还有一些可能存在相关性的事件,比如知识水平越高的人,收入水平越高;市场化的国家经济越好,则货币越强势,反而全球经济危机,黄金等避险资产越走强。

如果我们要研究这些事件,找到不同变量之间的关系,我们就会用到回归分析。一元线性回归分析是处理两个变量之间关系的最简单模型,是两个变量之间的线性相关关系。让我们一起发现生活中的规律吧。

由于本文为非统计的专业文章,所以当出现与教课书不符的描述,请以教课书为准。本文力求用简化的语言,来介绍一元线性回归的知识,同时配合R语言的实现。

目录

  1. 一元线性回归介绍
  2. 数据集和数学模型
  3. 回归参数估计
  4. 回归方程的显著性检验
  5. 残差分析和异常点检测
  6. 模型预测

1. 一元线性回归介绍

回归分析(Regression Analysis)是用来确定2个或2个以上变量间关系的一种统计分析方法。如果回归分析中,只包括一个自变量X和一个因变量Y时,且它们的关系是线性的,那么这种回归分析称为一元线性回归分析。

回归分析属于统计学的基本模型,涉及统计学基础,就会有一大堆的名词和知识点需要介绍。

在回归分析中,变量有2类:因变量 和 自变量。因变量通常是指实际问题中所关心的指标,用Y表示。而自变量是影响因变量取值的一个变量,用X表示,如果有多个自变量则表示为X1, X2, …, Xn。

回归分析研究的主要步骤:

  1. 确定因变量Y 与 自变量X1, X2, …, Xn 之间的定量关系表达式,即回归方程。
  2. 对回归方程的置信度检查。
  3. 判断自变量Xn(n=1,2,…,m)对因变量的影响。
  4. 利用回归方程进行预测。

本文会根据回归分析的的主要步骤,进行结构梳理,介绍一元线性回归模型的使用方法。

reg

2. 数据集和数学模型

先让我们通过一个例子开始吧,用一组简单的数据来说明一元线性回归分析的数学模型的原理和公式。找出下面数据集中Y与X的定量关系。

数据集为2016年3月1日,白天开盘的交易数据,为锌的2个期货合约的分钟线的价格数据。数据集包括有3列,索引列为时间,zn1.Close为ZN1604合约的1分钟线的报价数据,zn2.Close为ZN1605合约的1分钟线的报价数据。

数据集如下:


                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180
2016-03-01 09:07:00     14110     14170
2016-03-01 09:08:00     14110     14175
2016-03-01 09:09:00     14105     14170
2016-03-01 09:10:00     14105     14170
2016-03-01 09:11:00     14120     14180
2016-03-01 09:12:00     14105     14170
2016-03-01 09:13:00     14105     14170
2016-03-01 09:14:00     14110     14175
2016-03-01 09:15:00     14105     14175
2016-03-01 09:16:00     14120     14185
2016-03-01 09:17:00     14125     14190
2016-03-01 09:18:00     14115     14185
2016-03-01 09:19:00     14135     14195
2016-03-01 09:20:00     14125     14190
2016-03-01 09:21:00     14135     14205
2016-03-01 09:22:00     14140     14210
2016-03-01 09:23:00     14140     14200
2016-03-01 09:24:00     14135     14205
2016-03-01 09:25:00     14140     14205
2016-03-01 09:26:00     14135     14205
2016-03-01 09:27:00     14130     14205

我们以zn1.Close列的价格为X,zn2.Close列的价格为Y,那么试试找到自变量X和因变量Y的关系的表达式。

为了直观起见,我们可以先画出一张散点图,以X为横坐标,Y为纵坐标,每个点对应一个X和一个Y。


# 数据集已存在df变量中
> head(df)
                    zn1.Close zn2.Close
2016-03-01 09:01:00     14075     14145
2016-03-01 09:02:00     14095     14160
2016-03-01 09:03:00     14095     14160
2016-03-01 09:04:00     14095     14165
2016-03-01 09:05:00     14120     14190
2016-03-01 09:06:00     14115     14180

# 分别给x,y赋值
> x<-as.numeric(df[,1])
> y<-as.numeric(df[,2])

# 画图
> plot(y~x+1)

01

从散点图上发现 X和Y 的排列基本是在一条直线附近,那么我们可以假设X和Y的关系是线性,可以用公式表式为。

Y = a + b * X + c
  • Y,为因变量
  • X,为自变量
  • a,为截距
  • b,为自变量系数
  • a+b*X, 表示Y随X的变化而线性变化的部分
  • c, 为残差或随机误差,是其他一切不确定因素影响的总和,其值不可观测。假定c是符合均值为0方差为σ^2的正态分布 ,记作c~N(0,σ^2)

对于上面的公式,称函数f(X) = a + b * X 为一元线性回归函数,a为回归常数,b为回归系数,统称回归参数。X 为回归自变量或回归因子,Y 为回归因变量或响应变量。如果(X1,Y1),(X2,Y2)…(Xn,Yn)是(X,Y)的一组观测值,则一元线性回归模型可表示为


Yi = a + b * X + ci,     i= 1,2,...n

其中E(ci)=0, var(ci)=σ^2, i=1,2,...n

通过对一元线性回归模型的数学定义,接下来让我们利用数据集做回归模型的参数估计。

3. 回归参数估计

对于上面的公式,回归参数a,b是我们不知道的,我们需要用参数估计的方法来计算出a,b的值,而从得到数据集的X和Y的定量关系。我们的目标是要计算出一条直线,使直接线上每个点的Y值和实际数据的Y值之差的平方和最小,即(Y1实际-Y1预测)^2+(Y2实际-Y2预测)^2+ …… +(Yn实际-Yn预测)^2 的值最小。参数估计时,我们只考虑Y随X的线性变化的部分,而残差c是不可观测的,参数估计法并不需要考虑残差,对于残差的分析在后文中介绍。

令公式变形为a和b的函数Q(a,b), 即 (Y实际-Y测试)的平方和,变成到(Y实际 – (a+b*X))的平方和。

reg2

公式一 回归参数变形公式

通过最小二乘估计推导出a和b的求解公式,详细的推导过程请参考文章一元线性回归的细节

reg3

公式二 回归参数计算公式

其中 x和y的均值,计算方法如下

reg4

公式三 均值计算公式

有了这个公式,我们就可以求出a和b两个的回归参数的解了。

接下来,我们用R语言来实现对上面数据的回归模型的参数估计,R语言中可以用lm()函数来实现一元线性回归的建模过程。


# 建立线性回归模型
> lm.ab<-lm(y ~ 1+x)

# 打印参数估计的结果
> lm.ab

Call:
lm(formula = y ~ 1 + x)

Coefficients:
(Intercept)            x  
   -349.493        1.029  

如果你想动手来计算也可以自己实现公式。


# x均值
> Xm<-mean(x);Xm 
[1] 14034.82

# y均值
> Ym<-mean(y);Ym
[1] 14096.76

# 计算回归系数
> b <- sum((x-Xm)*(y-Ym)) / sum((x-Xm)^2) ;b
[1] 1.029315

# 计算回归常数
> a <- Ym - b * Xm;a
[1] -349.493

回归参数a和b的计算结果,与lm()函数的计算结果是相同的。有了a和b的值,我们就可以画出这条近似的直接线。

计算公式为:

Y= a + b * X = -349.493 + 1.029315 * X 

画出回归线。


> plot(y~x+1)
> abline(lm.ab)

02

这条直线是我们用数据拟合出来的,是一个近似的值。我们看到有些点在线上,有些点不在线上。那么要评价这条回归线拟合的好坏,我们就需要对回归模型进行显著性检验。

4. 回归方程的显著性检验

从回归参数的公式二可知,在计算过程中并不一定要知道Y和X是否有线性相关的关系。如果不存相关关系,那么回归方程就没有任何意义了,如果Y和X是有相关关系的,即Y会随着X的变化而线性变化,这个时候一元线性回归方程才有意义。所以,我们需要用假设检验的方法,来验证相关性的有效性。

通常会采用三种显著性检验的方法。

  • T检验法:T检验是检验模型某个自变量Xi对于Y的显著性,通常用P-value判断显著性,小于0.01更小时说明这个自变量Xi与Y相关关系显著。
  • F检验法:F检验用于对所有的自变量X在整体上看对于Y的线性显著性,也是用P-value判断显著性,小于0.01更小时说明整体上自变量与Y相关关系显著。
  • R^2(R平方)相关系统检验法:用来判断回归方程的拟合程度,R^2的取值在0,1之间,越接近1说明拟合程度越好。

在R语言中,上面列出的三种检验的方法都已被实现,我们只需要把结果解读。上文中,我们已经通过lm()函数构建一元线性回归模型,然后可以summary()函数来提取模型的计算结果。


> summary(lm.ab)      # 计算结果

Call:
lm(formula = y ~ 1 + x)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.9385  -2.2317  -0.1797   3.3546  10.2766 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.495e+02  7.173e+01  -4.872 2.09e-06 ***
x            1.029e+00  5.111e-03 201.390  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.232 on 223 degrees of freedom
Multiple R-squared:  0.9945,	Adjusted R-squared:  0.9945 
F-statistic: 4.056e+04 on 1 and 223 DF,  p-value: < 2.2e-16

模型解读:

  • Call,列出了回归模型的公式。
  • Residuals,列出了残差的最小值点,1/4分位点,中位数点,3/4分位点,最大值点。
  • Coefficients,表示参数估计的计算结果。
  • Estimate,为参数估计列。Intercept行表示常数参数a的估计值 ,x行表示自变量x的参数b的估计值。
  • Std. Error,为参数的标准差,sd(a), sd(b)
  • t value,为t值,为T检验的值
  • Pr(>|t|) ,表示P-value值,用于T检验判定,匹配显著性标记
  • 显著性标记,***为非常显著,**为高度显著, **为显著,·为不太显著,没有记号为不显著。
  • Residual standard error,表示残差的标准差,自由度为n-2。
  • Multiple R-squared,为相关系数R^2的检验,越接近1则越显著。
  • Adjusted R-squared,为相关系数的修正系数,解决多元回归自变量越多,判定系数R^2越大的问题。
  • F-statistic,表示F统计量,自由度为(1,n-2),p-value:用于F检验判定,匹配显著性标记。

通过查看模型的结果数据,我们可以发现通过T检验的截距和自变量x都是非常显著,通过F检验判断出整个模型的自变量是非常显著,同时R^2的相关系数检验可以判断自变量和因变量是高度相关的。

最后,我们通过的回归参数的检验与回归方程的检验,得到最后一元线性回归方程为:

Y = -349.493 + 1.029315 * X 

5. 残差分析和异常点检测

在得到的回归模型进行显著性检验后,还要在做残差分析(预测值和实际值之间的差),检验模型的正确性,残差必须服从正态分布N(0,σ^2)。

我们可以自己计算数据残差,并进行正态分布检验。


# 残差
> y.res<-residuals(lm.ab)

# 打印前6条数据
> head(y.res)
        1         2         3         4         5         6 
6.8888680 1.3025744 1.3025744 6.3025744 5.5697074 0.7162808 

# 正态分布检验
> shapiro.test(y.res)

	Shapiro-Wilk normality test

data:  y.res
W = 0.98987, p-value = 0.1164

# 画出残差散点图
> plot(y.res)

09

对残差进行Shapiro-Wilk正态分布检验,W接近1,p-value>0.05,证明数据集符合正态分布!关于正态分布的介绍,请参考文章常用连续型分布介绍及R语言实现

同时,我们也可以用R语言的工具生成4种用于模型诊断的图形,简化自己写代码计算的操作。


# 画图,回车展示下一张
> plot(lm.ab)    
Hit  to see next plot:   # 残差拟合图
Hit  to see next plot:   # 残差QQ图
Hit  to see next plot:   # 标准化的残差对拟合值 
Hit  to see next plot:   # 标准化残差对杠杆值

04

图1,残差和拟合值对比图

对残差和拟合值作图,横坐标是拟合值,纵坐标是残差。残差和拟合值之间,数据点均匀分布在y=0两侧,呈现出随机的分布,红色线呈现出一条平稳的曲线并没有明显的形状特征,说明残差数据表现非常好。

05

图2,残差QQ图

残差QQ图,用来描述残差是否符合正态分布。图中的数据点按对角直线排列,趋于一条直线,并被对角直接穿过,直观上符合正态分布。对于近似服从正态分布的标准化残差,应该有 95% 的样本点落在 [-2,2] 区间内。

06

图3,标准化残差平方根和拟合值对比图

对标准化残差平方根和拟合值作图,横坐标是拟合值,纵坐标是标准化后的残差平方根。与残差和拟合值对比图(图1)的判断方法类似,数据随机分布,红色线呈现出一条平稳的曲线,无明显的形状特征。

07

图4,标准残差和杠杆值对比图

对标准化残差和杠杆值作图,虚线表示的cooks距离等高线,通常用Cook距离度量的回归影响点。本图中没有出现红色的等高线,则说明数据中没有特别影响回归结果的异常点。

如果想把把4张图画在一起进行展示,可以改变画布布局。


> par(mfrow=c(2,2))
> plot(lm.ab)

08

看到上面4幅中,每幅图上都有一些点被特别的标记出来了,这些点是可能存在的异常值点,如果要对模型进行优化,我们可以从这些来入手。但终于本次残差分析的结果已经很好了,所以对于异常点的优化,可能并不能明显的提升模型的效果。

从图中发现,索引编号为27和192的2个点在多幅图中出现。我们假设这2个点为异常点,从数据中去掉这2个点,再进行显著性检验和残差分析。


# 查看27和192
> df[c(27,192),]
                    zn1.Close zn2.Close
2016-03-01 09:27:00     14130     14205
2016-03-01 14:27:00     14035     14085

# 新建数据集,去掉27和192
> df2<-df[-c(27,192),]

回归建模和显著性检验。


> x2<-as.numeric(df2[,1])
> y2<-as.numeric(df2[,2])
> lm.ab2<-lm(y2 ~ 1+x2)
> summary(lm.ab2)

Call:
lm(formula = y2 ~ 1 + x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0356 -2.1542 -0.2727  3.3336  9.5879 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.293e+02  7.024e+01  -4.688 4.83e-06 ***
x2           1.028e+00  5.004e-03 205.391  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.117 on 221 degrees of freedom
Multiple R-squared:  0.9948,	Adjusted R-squared:  0.9948 
F-statistic: 4.219e+04 on 1 and 221 DF,  p-value: < 2.2e-16

对比这次的显著性检验结果和之前结果,T检验,F检验 和 R^2检验,并没有明显的效果提升,结果和我预想的是一样的。所以,通过残差分析和异常点分析,我认为模型是有效的。

6. 模型预测

最后,我们获得了一元线性回归方程的公式,就可以对数据进行预测了。比如,对给定X=x0时,计算出y0=a+b*x0的值,并计算出置信度为1-α的预测区间。

当X=x0,Y=y0时,置信度为1-α的预测区间为

reg5

reg6

我们可以用R语言的predict()函数来计算预测值y0,和相应的预测区间。程序算法如下。


> new<-data.frame(x=14040)
> lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)

# 预测结果
> lm.pred
       fit      lwr      upr
1 14102.09 14093.73 14110.44

当x0=14040时,在预测区间为0.95的概率时,y0的值为 14102,预测区间为[14093.73,14110.44]。

我们通过图形来表示。


> plot(y~x+1)
> abline(lm.ab,col='red')
> points(rep(newX$x,3),y=lm.pred,pch=19,col=c('red','blue','green'))

03

其中,红色点为y0的值,蓝色点为预测区间最小值,绿色点为预测区间最大值。

对于统计模型中最核心部分就在结果解读,本文介绍了一元回归模型的基本的建模过程和模型的详细解读方法。在我们掌握了这种方法以后,就可以更容易地理解和学习 多元回归,非线性回归 等更多的模型,并把这些模型应用到实际的工作中了。下一篇文章R语言解读多元线性回归模型,请继续阅读。

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

打赏作者

R语言中文分词包jiebaR

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-word-jiebar/

jiebaR

前言

本文挖掘是数据挖掘中一个非常重要的部分,有非常广阔的使用场景,比如我们可以对新闻事件进行分析,了解国家大事;也可以对微博信息进行分析,通过社交舆情看看大家的关注点。通过文本挖掘找到文章中的隐藏信息,对文章的结构进行分析,判断是不是同一个作者写文章;同时可以对邮件分析,结合bayes算法判断哪些是垃圾邮件,哪些是有用的邮件。

本文挖掘的第一步,就是要进行分词,分词将直接影响文本挖掘的效果。R语言在分词方面有很好的支持,接下来就给大家介绍一个不错的R语言中文分词包“结巴分词”(jiebaR)。

目录

  1. jiebaR包介绍
  2. 5分钟上手
  3. 分词引擎
  4. 配置词典
  5. 停止词过滤
  6. 关键词提取

1. jiebaR包介绍

结巴分词(jiebaR),是一款高效的R语言中文分词包,底层使用的是C++,通过Rcpp进行调用很高效。结巴分词基于MIT协议,就是免费和开源的,感谢国人作者的给力支持,让R的可以方便的处理中文文本。

官方Github的地址:https://github.com/qinwf/jiebaR

本文所使用的系统环境

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

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


~ R
> install.packages("jiebaR")
> library("jiebaR")

如果想要安装开发版本,可以使用devtools来进行安装,devtools的介绍请参考文章:在巨人的肩膀前行 催化R包开发


> library(devtools)
> install_github("qinwf/jiebaRD")
> install_github("qinwf/jiebaR")
> library("jiebaR")

开发版本安装,官方建议使用Linux系统 gcc >= 4.6 编译,Windows需要安装 Rtools。

2. 5分钟上手

5分钟上手,直接看第一个例子吧,对一段文字进行分词。


> wk = worker()

> wk["我是《R的极客理想》图书作者"]
[1] "我是" "R"    "的"   "极客" "理想" "图书" "作者"

> wk["我是R语言的深度用户"]
[1] "我"   "是"   "R"    "语言" "的"   "深度" "用户"

很简单地,2行代码,就完成了中文分词。

jiebaR提供了3种分词语句的写法,例子上面的用[]符号的语法,还可以使用<=符合语法,或者使用segment()函数。虽然形式不同,但是分词效果是一样的。 使用<=符号的语法,如下


> wk<='另一种符合的语法'
[1] "另"   "一种" "符合" "的"   "语法"

使用segment()函数的语法,如下


> segment( "segment()函数语句的写法" , wk )
[1] "segment" "函数"    "语句"    "的"      "写法" 

如果你觉得很神奇,想了解如何自定义操作符的,可以检查项目的源代码quick.R文件


# <= 符号定义
`<=.qseg`<-function(qseg, code){
  if(!exists("quick_worker",envir = .GlobalEnv ,inherits = F) || 
       .GlobalEnv$quick_worker$PrivateVarible$timestamp != TIMESTAMP){
    
    if(exists("qseg",envir = .GlobalEnv,inherits = FALSE ) ) 
      rm("qseg",envir = .GlobalEnv)
    
    modelpath  = file.path(find.package("jiebaR"),"model","model.rda")
    quickparam = readRDS(modelpath)
    
    if(quickparam$dict == "AUTO") quickparam$dict = DICTPATH
    if(quickparam$hmm == "AUTO") quickparam$hmm = HMMPATH
    if(quickparam$user == "AUTO") quickparam$user = USERPATH
    if(quickparam$stop_word == "AUTO") quickparam$stop_word = STOPPATH
    if(quickparam$idf == "AUTO") quickparam$idf = IDFPATH
    
    createquickworker(quickparam)
    setactive()
  } 

  //..代码省略
}

# [ 符号定义
`[.qseg`<- `<=.qseg`

我们也可以直接对文本文件进行分词,在当前目录新建一个文本文件idea.txt。


~ notepad idea.txt

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

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

当然,我们运行分词程序,会在当前目录生成一个新的分词结果的文件。


> wk['./idea.txt']
[1] "./idea.segment.2016-07-20_23_25_34.txt"

打开文件idea.segment.2016-07-20_23_25_34.txt,整个本文以空格进行分词。


~ notepad idea.segment.2016-07-20_23_25_34.txt

R 的 极客 理想 系列 文章 涵盖 了 R 的 思想 使用 工具 创新 等 的 一系列 要点 以 我 个人 的 学习 和 体验 去 诠释 R 的 强大 R 语言 作为 统计学 一门 语言 一直 在 小众 领域 闪耀着 光芒 直到 大 数据 的 爆发 R 语言 变成 了 一门 炙手可热 的 数据分析 的 利器 随着 越来越 多 的 工程 背景 的 人 的 加入 R 语言 的 社区 在 迅速 扩大 成长 现在 已 不仅仅 是 统计 领域 教育 银行 电商 互联网 都 在 使用 R 语言

是不是很简单,5分钟实践就能完成分词的任务。

3. 分词引擎

在调用worker()函数时,我们实际是在加载jiebaR库的分词引擎。jiebaR库提供了7种分词引擎。

  • 混合模型(MixSegment):是四个分词引擎里面分词效果较好的类,结它合使用最大概率法和隐式马尔科夫模型。
  • 最大概率法(MPSegment) :负责根据Trie树构建有向无环图和进行动态规划算法,是分词算法的核心。
  • 隐式马尔科夫模型(HMMSegment):是根据基于人民日报等语料库构建的HMM模型来进行分词,主要算法思路是根据(B,E,M,S)四个状态来代表每个字的隐藏状态。 HMM模型由dict/hmm_model.utf8提供。分词算法即viterbi算法。
  • 索引模型(QuerySegment):先使用混合模型进行切词,再对于切出来的较长的词,枚举句子中所有可能成词的情况,找出词库里存在。
  • 标记模型(tag)
  • Simhash模型(simhash)
  • 关键词模型(keywods)

如果你不太关心引擎的事,那么直接用官方推荐的混合模型(默认选择)就行了。查看worker()函数的定义。


worker(type = "mix", dict = DICTPATH, hmm = HMMPATH, user = USERPATH,
  idf = IDFPATH, stop_word = STOPPATH, write = T, qmax = 20, topn = 5,
  encoding = "UTF-8", detect = T, symbol = F, lines = 1e+05,
  output = NULL, bylines = F, user_weight = "max")

参数列表:

  • type, 引擎类型
  • dict, 系统词典
  • hmm, HMM模型路径
  • user, 用户词典
  • idf, IDF词典
  • stop_word, 关键词用停止词库
  • write, 是否将文件分词结果写入文件,默认FALSE
  • qmax, 最大成词的字符数,默认20个字符
  • topn, 关键词数,默认5个
  • encoding, 输入文件的编码,默认UTF-8
  • detect, 是否编码检查,默认TRUE
  • symbol, 是否保留符号,默认FALSE
  • lines, 每次读取文件的最大行数,用于控制读取文件的长度。大文件则会分次读取。
  • output, 输出路径
  • bylines, 按行输出
  • user_weight, 用户权重

我们在调用worker()时,就加载了分词引擎,可以打印出来,查看分词的引擎的配置。


> wk = worker()
> wk
Worker Type:  Jieba Segment

Default Method  :  mix     # 混合模型
Detect Encoding :  TRUE    # 检查编码
Default Encoding:  UTF-8   # UTF-8
Keep Symbols    :  FALSE   # 不保留符号
Output Path     :          # 输出文件目录
Write File      :  TRUE    # 写文件
By Lines        :  FALSE   # 不行输出
Max Word Length :  20      # 最大单单词长度
Max Read Lines  :  1e+05   # 最大读入文件行数

Fixed Model Components:  

$dict                      # 系统词典
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/jieba.dict.utf8"

$user                      # 用户词典
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/user.dict.utf8"

$hmm                       # 隐式马尔科夫模型模型
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict/hmm_model.utf8"

$stop_word                 # 停止词,无
NULL

$user_weight               # 用户词典权重
[1] "max"

$timestamp                 # 时间戳
[1] 1469027302

$default $detect $encoding $symbol $output $write $lines $bylines can be reset.

如果我们想改变分词引擎的配置项,可以在调用worker()创建分词引擎时,也可以通过wk$XX来进行设置。如果想了解wk是什么类型的对象,我们通过pryr包的otype的函数来检查wk对象的类型。关于pryr包的详细使用,请参考文章撬动R内核的高级工具包pryr


# 加载 pryr包
> library(pryr)
> otype(wk)  # 面向对象的类型检查
[1] "S3"

> class(wk)  # 查看class是属性
[1] "jiebar"  "segment" "jieba" 

4. 配置词典

对于分词的结果好坏的关键因素是词典,jiebaR默认有配置标准的词典。对于我们的使用来说,不同行业或不同的文字类型,最好用专门的分词词典。在jiebaR中通过show_dictpath()函数可以查看默认的标准词典,可以通过上一小节介绍的配置项,来指定我们自己的词典。日常对话的常用词典,比如搜狗输入法的词库。


# 查看默认的词库位置
> show_dictpath()
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict"

# 查看目录
> dir(show_dictpath())
[1] "D:/tool/R-3.2.3/library/jiebaRD/dict"
 [1] "backup.rda"      "hmm_model.utf8"  "hmm_model.zip"  
 [4] "idf.utf8"        "idf.zip"         "jieba.dict.utf8"
 [7] "jieba.dict.zip"  "model.rda"       "README.md"      
[10] "stop_words.utf8" "user.dict.utf8" 

看到词典目录中,包括了多个文件。

  • jieba.dict.utf8, 系统词典文件,最大概率法,utf8编码的
  • hmm_model.utf8, 系统词典文件,隐式马尔科夫模型,utf8编码的
  • user.dict.utf8, 用户词典文件,utf8编码的
  • stop_words.utf8,停止词文件,utf8编码的
  • idf.utf8,IDF语料库,utf8编码的
  • jieba.dict.zip,jieba.dict.utf8的压缩包
  • hmm_model.zip,hmm_model.utf8的压缩包
  • idf.zip,idf.utf8的压缩包
  • backup.rda,无注释
  • model.rda,无注释
  • README.md,说明文件

打开系统词典文件jieba.dict.utf8,并打印前50行。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/jieba.dict.utf8",
+           what=character(),nlines=50,sep='\n',
+           encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "1号店 3 n"  "1號店 3 n"  "4S店 3 n"   "4s店 3 n"  
 [5] "AA制 3 n"   "AB型 3 n"   "AT&T 3 nz"  "A型 3 n"   
 [9] "A座 3 n"    "A股 3 n"    "A輪 3 n"    "A轮 3 n"   
[13] "BB机 3 n"   "BB機 3 n"   "BP机 3 n"   "BP機 3 n"  
[17] "B型 3 n"    "B座 3 n"    "B股 3 n"    "B超 3 n"   
[21] "B輪 3 n"    "B轮 3 n"    "C# 3 nz"    "C++ 3 nz"  
[25] "CALL机 3 n" "CALL機 3 n" "CD机 3 n"   "CD機 3 n"  
[29] "CD盒 3 n"   "C座 3 n"    "C盘 3 n"    "C盤 3 n"   
[33] "C語言 3 n"  "C语言 3 n"  "D座 3 n"    "D版 3 n"   
[37] "D盘 3 n"    "D盤 3 n"    "E化 3 n"    "E座 3 n"   
[41] "E盘 3 n"    "E盤 3 n"    "E通 3 n"    "F座 3 n"   
[45] "F盘 3 n"    "F盤 3 n"    "G盘 3 n"    "G盤 3 n"   
[49] "H盘 3 n"    "H盤 3 n"

我们发现系统词典每一行都有三列,并以空格分割,第一列为词项,第二列为词频,第三列为词性标记。

打开用户词典文件user.dict.utf8,并打印前50行。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/user.dict.utf8",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 5 items
[1] "云计算"   "韩玉鉴赏" "蓝翔 nz"  "CEO"      "江大桥"  

用户词典第一行有二列,,第一列为词项,第二列为词性标记,没有词频的列。用户词典默认词频为系统词库中的最大词频。

jiebaR包关于词典词性标记,采用ictclas的标记方法。ICTCLAS 汉语词性标注集。

代码 名称 帮助记忆的诠释
Ag 形语素 形容词性语素。形容词代码为a,语素代码g前面置以A。
a 形容词 取英语形容词adjective的第1个字母。
ad 副形词 直接作状语的形容词。形容词代码a和副词代码d并在一起。
an 名形词 具有名词功能的形容词。形容词代码a和名词代码n并在一起。
b 区别词 取汉字"别"的声母。
c 连词 取英语连词conjunction的第1个字母。
Dg 副语素 副词性语素。副词代码为d,语素代码g前面置以D。
d 副词 取adverb的第2个字母,因其第1个字母已用于形容词。
e 叹词 取英语叹词exclamation的第1个字母。
f 方位词 取汉字"方"的声母。
g 语素 绝大多数语素都能作为合成词的"词根",取汉字"根"的声母。
h 前接成分 取英语head的第1个字母。
i 成语 取英语成语idiom的第1个字母。
j 简称略语 取汉字"简"的声母。
k 后接成分
l 习用语 习用语尚未成为成语,有点"临时性",取"临"的声母。
m 数词 取英语numeral的第3个字母,n,u已有他用。
Ng 名语素 名词性语素。名词代码为n,语素代码g前面置以N。
n 名词 取英语名词noun的第1个字母。
nr 人名 名词代码n和"人(ren)"的声母并在一起。
ns 地名 名词代码n和处所词代码s并在一起。
nt 机构团体 "团"的声母为t,名词代码n和t并在一起。
nz 其他专名 "专"的声母的第1个字母为z,名词代码n和z并在一起。
o 拟声词 取英语拟声词onomatopoeia的第1个字母。
p 介词 取英语介词prepositional的第1个字母。
q 量词 取英语quantity的第1个字母。
r 代词 取英语代词pronoun的第2个字母,因p已用于介词。
s 处所词 取英语space的第1个字母。
Tg 时语素 时间词性语素。时间词代码为t,在语素的代码g前面置以T。
t 时间词 取英语time的第1个字母。
u 助词 取英语助词auxiliary 的第2个字母,因a已用于形容词。
Vg 动语素 动词性语素。动词代码为v。在语素的代码g前面置以V。
v 动词 取英语动词verb的第一个字母。
vd 副动词 直接作状语的动词。动词和副词的代码并在一起。
vn 名动词 指具有名词功能的动词。动词和名词的代码并在一起。
w 标点符号
x 非语素字 非语素字只是一个符号,字母x通常用于代表未知数、符号。
y 语气词 取汉字"语"的声母。
z 状态词 取汉字"状"的声母的前一个字母。

下面我们自定义一个用户词典,来试试效果。编写词典文件,user.utf8。


~ notepad user.utf8

R语言
R的极客理想
大数据
数据

使用我们的自定义的用户词典,对刚才的文本再进行分词。


> wk = worker(user='user.utf8')
> wk['./idea.txt']
[1] "./idea.segment.2016-07-21_11_14_24.txt"

对比2次产生的分词结果,idea.segment.2016-07-20_23_25_34.txt 和 idea.segment.2016-07-21_11_14_24.txt。

jiebaR-cut

在实际使用中,jiebaR默认提供的用户词典只有5个单词,太简单了,肯定是不够用的。我们可以用搜狗词典,来丰富用户自己的词库。接下来,让我们配置搜狗词典。你需要安装一个搜狗输入法,具体的安装过程不再解释。

我安装的是搜狗五笔输入法,找到搜狗的安装目录,并找到词典文件。我的搜狗词典,在下面的安装位置。


C:\Program Files (x86)\SogouWBInput\2.1.0.1288\scd\17960.scel

把17960.scel文件复制到自己的项目目录里,用文本编辑器打开文件,发现是二进制的。那么我需要用工具进行转换,把二进制的词典转成我们可以使用的文本文件。jiebaR包的作者,同时开发了一个cidian项目,可以转换搜狗的词典,那么我们只需要安装cidian包即可。

安装cidian项目


> install.packages("devtools")
> install.packages("stringi")
> install.packages("pbapply")
> install.packages("Rcpp")
> install.packages("RcppProgress")
> library(devtools)
> install_github("qinwf/cidian")
> library(cidian)

转换二进制词典到文本文件。


# 转换
> decode_scel(scel = "./17960.scel",cpp = TRUE)
output file: ./17960.scel_2016-07-21_00_22_11.dict

# 查看生成的词典文件
> scan(file="./17960.scel_2016-07-21_00_22_11.dict",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "阿坝州 n"         "阿百川 n"         "阿班 n"          
 [4] "阿宾 n"           "阿波菲斯 n"       "阿不都热希提 n"  
 [7] "阿不都西库尔 n"   "阿不力克木 n"     "阿尔姆格伦 n"    
[10] "阿尔沙文 n"       "阿肥星 n"         "阿菲正传 n"      
[13] "阿密特 n"         "阿穆 n"           "阿穆隆 n"        
[16] "阿帕鲁萨镇 n"     "阿披实 n"         "阿衰 n"          
[19] "阿霞 n"           "艾奥瓦 n"         "爱不疚 n"        
[22] "爱的错位 n"       "爱得得体 n"       "爱的火焰 n"      
[25] "爱的流刑地 n"     "爱得起 n"         "埃夫隆 n"        
[28] "爱搞网 n"         "爱国红心 n"       "爱呼 n"          
[31] "爱就宅一起 n"     "埃克希儿 n"       "爱没有错 n"      
[34] "埃蒙斯 n"         "爱奴新传 n"       "爱起点 n"        
[37] "爱情的牙齿 n"     "爱情海滨 n"       "爱情节 n"        
[40] "爱情美的样子 n"   "爱情无限谱 n"     "爱情占线 n"      
[43] "爱情转移 n"       "爱情左灯右行 n"   "爱上你是一个错 n"
[46] "矮哨兵 n"         "爱是妥协 n"       "爱似水仙 n"      
[49] "爱太痛 n"         "爱无界 n"    

接下来,直接把搜狗词典配置到我们的分词库中,就可以直接使用了。把搜狗词典文件改名,从17960.scel_2016-07-21_00_22_11.dict到user.dict.utf8,然后替换D:\tool\R-3.2.3\library\jiebaRD\dict目录下面的user.dict.utf8。这样默认的用户词典,就是搜狗词典了。很酷吧!

5. 停止词过滤

停止词就是分词过程中,我们不需要作为结果的词,像英文的语句中有很多的a,the,or,and等,中文语言中也有很多,比如 的,地,得,我,你,他。这些词因为使用频率过高,会大量出现在一段文本中,对于分词后的结果,在统计词频的时候会增加很多的噪音,所以我们通常都会将这些词进行过滤。

在jiebaR中,过滤停止词有2种方法,一种是通过配置stop_word文件,另一种是使用filter_segment()函数。

首先我们先来看,通过配置stop_word文件的方法。新建一个stop_word.txt文件。


~ notepad stop_word.txt

我
我是

加载分词引擎,并配置停止词过滤。


> wk = worker(stop_word='stop_word.txt')
> segment<-wk["我是《R的极客理想》图书作者"]
> segment
[1] "R"    "的"   "极客" "理想" "图书" "作者"

上面的文本,我们把"我是"通过停止词进行了过滤。如果还想过滤“作者”一词,可以动态的调用filter_segment()函数。


> filter<-c("作者")
> filter_segment(segment,filter)
[1] "R"    "的"   "极客" "理想" "图书"

6. 关键词提取

关键词提取是文本处理非常重要的一个环节,一个经典算法是TF-IDF算法。其中,TF(Term Frequency)代表词频,IDF(Inverse Document Frequency)表示逆文档频率。如果某个词在文章中多次出现,而且不是停止词,那么它很可能就反应了这段文章的特性,这就是我们要找的关键词。再通过IDF来算出每个词的权重,不常见的词出现的频率越高,则权重越大。计算TF-IDF的公式为:

TF-IDF = TF(词频) * 逆文档频率(IDF)

对文档中每个词计算TF-IDF的值,把结果从大到小排序,就得到了这篇文档的关键性排序列表。关于IF-IDF的解释,参考了文章TF-IDF与余弦相似性的应用(一):自动提取关键词

jiebaR包的关键词提取提取的实现,也是使用了TF-IDF的算法。在安装目录中的idf.utf8文件,为IDF的语料库。查看idf.utf8内容。


> scan(file="D:/tool/R-3.2.3/library/jiebaRD/dict/idf.utf8",
+      what=character(),nlines=50,sep='\n',
+      encoding='utf-8',fileEncoding='utf-8')
Read 50 items
 [1] "劳动防护 13.900677652"      "生化学 13.900677652"       
 [3] "奥萨贝尔 13.900677652"      "考察队员 13.900677652"     
 [5] "岗上 11.5027823792"         "倒车档 12.2912397395"      
 [7] "编译 9.21854642485"         "蝶泳 11.1926274509"        
 [9] "外委 11.8212361103"         "故作高深 11.9547675029"    
[11] "尉遂成 13.2075304714"       "心源性 11.1926274509"      
[13] "现役军人 10.642581114"      "杜勃留 13.2075304714"      
[15] "包天笑 13.900677652"        "贾政陪 13.2075304714"      
[17] "托尔湾 13.900677652"        "多瓦 12.5143832909"        
[19] "多瓣 13.900677652"          "巴斯特尔 11.598092559"     
[21] "刘皇帝 12.8020653633"       "亚历山德罗夫 13.2075304714"
[23] "社会公众 8.90346537821"     "五百份 12.8020653633"      
[25] "两点阈 12.5143832909"       "多瓶 13.900677652"         
[27] "冰天 12.2912397395"         "库布齐 11.598092559"       
[29] "龙川县 12.8020653633"       "银燕 11.9547675029"        
[31] "历史风貌 11.8212361103"     "信仰主义 13.2075304714"    
[33] "好色 10.0088573539"         "款款而行 12.5143832909"    
[35] "凳子 8.36728816325"         "二部 9.93038573842"        
[37] "卢巴 12.1089181827"         "五百五 13.2075304714"      
[39] "畅叙 11.598092559"          "吴栅子 13.2075304714"      
[41] "智力竞赛 13.900677652"      "库邦 13.2075304714"        
[43] "非正义 11.3357282945"       "编订 10.2897597393"        
[45] "悲号 12.8020653633"         "陈庄搭 13.2075304714"      
[47] "二郎 9.62401153296"         "电光石火 11.8212361103"    
[49] "抢球 11.9547675029"         "南澳大利亚 10.9562386728"  

idf.utf8文件每一行有2列,第一列是词项,第二列为权重。然后,我通过计算文档的词频(TF),与语料库的IDF值相乘,就可以得到TF-IDF值,从而提取文档的关键词。

比如,我们对下面的文本内容进行关键词的提取。


> wk = worker()
> segment<-wk["R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。"]

# 计算词频
> freq(segment)
     char freq
1    创新    1
2      了    1
3    文章    1
4    强大    1
5       R    3
6    个人    1
7      的    5
8    诠释    1
9      和    1
10 一系列    1
11   使用    1
12     以    1
13     等    1
14   极客    1
15   理想    1
16   思想    1
17   涵盖    1
18   系列    1
19     去    1
20     我    1
21   工具    1
22   学习    1
23   体验    1
24   要点    1

# 取TF-IDF的前5的关键词
> keys = worker("keywords",topn=5)

# 计算关键词
> vector_keywords(segment,keys)
11.7392 8.97342 8.23425  8.2137 7.43298 
 "极客"  "诠释"  "要点"  "涵盖"  "体验" 

使用jiebaR包处理分词确实简单,几行的代码就能实现分词的各种算法操作。有了这个工具,我们就可以文档中,发现各种语言规则进行文本挖掘了。下篇文章让我们挖掘一下上市公司的公告吧,说不定能发现什么市场规则。

本文只是抛砖引玉地介绍了jiebaR包的使用方法,详细使用操作,请参考包作者的官方介绍。再次感谢jiebaR作者@qinwenfeng,为R语言在中文分词中提供了一套非常不错的工具包!

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

打赏作者

2016天善智能交流会第22场: R语言为量化而生

跨界知识聚会系列文章,“知识是用来分享和传承的”,各种会议、论坛、沙龙都是分享知识的绝佳场所。我也有幸作为演讲嘉宾参加了一些国内的大型会议,向大家展示我所做的一些成果。从听众到演讲感觉是不一样的,把知识分享出来,你才能收获更多。

关于作者

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

转载请注明出处:
http://blog.fens.me/meeting-hellobi-20160701/

meeting-hellobi

前言

感谢天善智能社区的邀请,有幸参加每周一期的跟数据有关的行业、工具、技术的交流盛宴,活动的口号是“Friday BI Fly 周五BI飞起来”。

目录

  1. 我的分享主题:R语言为量化而生
  2. 会议体验
  3. 自由讨论

1. 我的分享主题:R语言为量化而生

本次分享的主题 R语言为量化而生,主要内容来自我的一篇博客文章:R语言为量化而生。希望能够解释清楚,在量化投资中为什么要用R语言。从程序员的角度看,C++,Java,Python, C#都是可行方案;从数据人员的角度看,Excel, SAS, Matlab更是不错的。那么为什么是R语言呢,R语言的优势在哪里体现?

这类的问题,总是会被问到。那么答案,就在于你对量化这件事情的了解,和对各种编程语言的理解。最近3年,互联网在量化领域的大发展,以Quantopian为代表的在线策略研发平台,用Python做为核心语言,国内同样支持Python的平台也有 优矿聚宽米筐。这些平台主是面向程序员群体的平台,希望通过挖掘草根明星,来推动量化的发展。传统的量化交易软件,像文华MC, TB, TS 都有自己一套的脚本化的编程语言。有实力的专业团队,通常会自成体系的独立开发一套自己的系统。如果面向更广泛的人群,最常用的方法就是Wind导数据,Excel中拉个表出来。

所以,其实用什么语言不重要,关键是怎么理解做量化这件事情。那么R语言的天生优势就是数学计算,数据处理,免费开源,大量支持库。试试吧,你一定会喜欢的。

2. 会议体验

本次分享受天善智能社区的邀请,我真的非常高兴。天善智能是新一代的商业智能和大数据的垂直社区,聚集了大量的数据分析从业人员。活动介绍,https://ask.hellobi.com/blog/tianshansoft/4229。 本次活动同时有30个微信群进行直播,参加的人员,至少有2000人以上。可以天善智能社区,在行业的影响力是非常大的。

发个截图,体会一下微信同步直播的震撼吧!

wx

本此的分享基于微信的直播,我也第一次体验,要用纯文字的方式来进行介绍。想把一个事情说清楚,又增加了不少的难度。由于不能分享屏幕,代码部分会通过图片截屏。

本次活动的总结,https://ask.hellobi.com/blog/tianshansoft/4271,感谢天善社区的工作人员进行整理。

远程分享,就是没能与大家合照,有点遗憾!!贴张自己的照片吧。

01

3. 自由讨论

分享后,很多朋友都对于R语言都是非常的好奇,提了很多的问题,用户的参与性非常强。下列直接贴出用户的问题和我的回复。

1、替新手问一个,请教一下,R语言的数据分析应该从哪方面入手练习啊?因为目前工作上不是用R的,看完书之后想具体去试一下。

张丹: R其实上手很快,找一本书,认真操作练习一遍就上手了。

2、玉琴:不建议用for loop的原因是考虑到性能问题吗

张丹:for loop是调用的R的循环库,apply是调用C的循环库,性能差距还是很大的

3、来自20群的提问:提个问题,微软对R的收购会对R语言的发展产生什么影响?

张丹:我觉得这是正向发展的,是好事情。大公司看到了R的潜力!

4、尚林栋:R语言金融建模的具体步骤能说一下吗

金融建模的具体步骤,你可以参考这篇文章,http://blog.fens.me/finance-stock-ma/

5、刘嘉丰Alan:丹哥,现在有很多量化平台,提供打包好的函数,在线回测,和自己造轮子拿R语言相比,您觉得各有什么优势呢?

张丹:R的优势就是在数学计算,数据处理上。行业标准还没有统一,所以不一定在线平台的轮子就一定好用。但另外,我们从开发或使用的角度,更多的用到的R包,都是RStudio公司的产品,我觉得是RStudio在推动R的整个的进化过程。

6、我也觉得r语言不错,但经常想不到商业场景,到现在,我只是用它统计考勤,各种绩效kpi,每月算一次奖金,已经这样过去2年了,r语言路在何方哪?

张丹: 你所说的统计,只能说简单计数。比如,你要预测下个月的考勤情况,从而设计预算方案。你可能就需要做个回归分析,这时R就能给你很大的帮助了。生活和工作中,随处都是数据分析的场景。

7、Allen:r在拟合上感觉比python用起来更爽一些,其返回的结果较多

张丹:那么R和python比,R更面向数据,特别是对于没有编程基础的人。PYTHON,还是程序语言,还要了解程序结构,程序架构,代码量不会少。

有IT背景程序员,可能更倾向于PYTHON;如果没有IT背景,R更容易上手。

8、越中女儿:请教一个问题:quantmod对美股的实时接口很好用,对A股不支持,且A股基本面数据才更新到2013.09,请问有好用的ETL包么,类似于python的tushare那样对A股友好的,各种etl啊清洗的脏活累活感觉python更好啊,R就是安安静静做做统计,玩玩图形。

张丹: quantmod使用的是yahoo等国外的数据源,这些数据源本身没有A股数据,如果需要A股数据,用tushare还是不错的。 R特有的data.frame,matrix 等类型和操作方法,在python也需要单独去实现。

9、柠檬味的香草:最近想研究一些互联网文本数据与指数或各股走势的关系,但是在使用R语言处理文本数据不是很方便,丹哥可有一些强大的library推荐,对于非结构,文本数据的处理。

张丹:“尽量使用向量计算或矩阵计算的计算方法”,可以这样理解,对于一个二维结构,for需要2次,0(N^2)的时间复杂度。如果我们把数据,直接就按矩阵存储, 你让矩阵里的每个点都加1, 只需要算一次。Hadley提供的包,源代码我都看过,写很棒,也很实用。

r在拟合上感觉比python用起来更爽一些,其返回的结果较多

其实R有很多的第三方的包,已经有了大量的算法包,而其他语言相对较少。只是我们平时接触的不多,所以觉得用不到。R有大量的统计包,你可以从官方网站找到,输出的结果,大部分也都是统计的结果。

R所支持的行业领域,非常广泛。而工程的语言,不会做细粒度的区分,只是通用的解决方法。

10、郑州—金融数据:python有pandas.DataFrame,pandas应该是第三方的数据库结构吧?R的data.frame是内置的。

张丹:pandas.DataFrame,在底层处理,还需要对原PYTHON的数据结构做映射。当然他可以解决的很好,但你看到的内存结构,可能并不是真正的内存结构。

R内置数据类型,就可以理解是内存结构。不需要再考虑转换了。找一个自己熟悉的语言,大多数的功能,每种语言都是能实现。只有很细的领域,才会进一步区分。

11、RHaoop采用分布式并行计算,那请问如何解决需要嵌套循环的算法。

张丹:对于基于hadoop大数据的MR计算,建议做数学变成,通过数学的角度处理。我写过2个例子,一个是pagerank, 一个是itemcf。

12、@柠檬味的香草:想听听丹哥对传统数据挖掘转量化投资的建议。比如前景?竞争力?

张丹:量化投资,其实是IT人都想转的行业。你写的代码,不是通过工资来赚钱,而直接通过交易赚钱,代码的效用是最大化的。但这个行业竞争很大,聪明人都在这里,要么你的技术牛,要么你了解市场,要么的算法是独特的,不然也很难。

JhT: 做量化交易和策略的都是高智商的

越中女儿:我觉得量化对金融市场的理解比对技术本身更重要,R的需求应该会很快凸显出来。因为数据基础都有了,后面就是差会分析的人了。通常懂数据分析的程序员,比纯程序员待遇高。

13、老师,有好的spark或者hadoop入门的书吗,计算机能力弱和java不懂啊

张丹:hadoop有很多书了,我当初看的是 权威指南。spark的书不了解,我的是网上文档。

14、@Mia.W 学RHadoop需要对Hadoop或Mapreduce了解到什么程度,需要从头学hadoop或java吗

张丹:hadoop的MR的原理要了解,找到懂JAVA的同事,帮你把环境搭好。

15、@JhT 我是刚进来的,R的优势是什么?

张丹:R是免费开源的,CRAN上有8000多个包,遍布各行各业。R语言的3个特性,数学计算,数据建模,可视化。

16、@郑州—金融数据个人感觉商业上matlab比R和python支持度都要好,不管是分析,统计,挖掘还是量化方便,收费的毕竟是收费的

张丹:有商业推动,当然要比免费的好了。不过,像SAS和Matlab也在打通和R的接口,毕竟由全球第三方贡献包,要比一家公司提供的包要多很多的。

17、@越中女儿 有用R做过实盘风控么

张丹:有做,其实不太复杂。你把需要的实时数据,都同步存到redis中,用R在秒级调reids取数据,计算完成再写回去。

18、@Jason.k计算机8g内存,数据虽然行数不多,但是很多列,所以数据csv格式大小会高达几个G,这个规模数据量,内存应该是不够的。

张丹:R的机制,会把数据一次性加载到内存中。就算能读到内存,每次计算时,也会有中间变量,所以你的基础内存是不够的。而且对于win性能会更差。

最后,再次感谢 天善社区的小伙伴们的努力,谢谢大家!

转载请注明出处:
http://blog.fens.me/meeting-hellobi-20160701/

打赏作者