• Archive by category "程序算法"
  • (Page 2)

Blog Archives

均值回归,逆市中的投资机会

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

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

关于作者:

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

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

meanReversion

前言

在股票市场中有两种典型的投资策略:趋势追踪(Trend Following) 和 均值回归(Mean Reversion)。 趋势追踪策略的特点在大行情的波动段找到有效的交易信号,不仅简单而且有效,我之前写的一篇文章 两条均线打天下 就属于趋势追踪策略。而均值回归策略则是一种反趋势策略,一波大幅上涨后容易出现下跌,而一波大幅下跌后容易出现上涨。其特点在振荡的在震荡的市场中非常有效,捕捉小的机会,本文就将介绍这种策略。

目录

  1. 均值回归原理
  2. 均值回归模型和实现
  3. 量化选股

1. 均值回归原理

在金融学中,均值回归是价格偏离均衡价格水平一定程度后向均衡价格靠拢的规律。本质上,均值回归就是哲学思想中所说的物极必反,可以简单地概括为“涨多必跌,跌多必涨”的规律。

均值回归是指股票价格无论高于或低于均值(均衡价格水平)都会以很高的概率向均值回归。根据这个理论,股票价格总是围绕其均值上下波动。一种上涨或者下跌的趋势不管其延续的时间多长都不能永远持续下去,最终均值回归的规律一定会出现:涨得太多了,就会向均值移动下跌;跌得太多了,就会向均值移动上升。如果我们认为事物总要回归常态,并且基于这样的预期来做任何决策的时候,我们就是在应用均值回归的理论。

下面以平安银行(000001)股票日K线图为例,可以非常直观的了解均值回归这种现象, 截取2005年到2015年7月的股票数据,股价为向前复权的价格。

01

上图中有3条曲线,黑色线是平安银行向前复权后的每日股价,红色线为20日均线,蓝色线为60日均线。关于均线的介绍,请参考文章 两条均线打天下。图中还有一条红色的水平线虚线,是这10年的股价平均值等于7.14元。这10年间,平安银行的股价经历了几波上涨和下跌,多次穿越7.14平均值。那么这个现象就是我们要讨论的均值回归。

1.1 均值回归的3个特性

均值回归是价值投资理论成立的一个核心理论,具有3个特性:必然性、不对称性、政府调控。

必然性,股票价格不能总是上涨或下跌,一种趋势不管其持续的时间多长都不能永远持续下去。在一个趋势内,股票价格呈持续上升或下降,我们称之为均值回避(Mean Aversion)。当出现相反趋势时就呈均值回归(Mean Reversion),但回归的周期有随机性是我们不能预测。不同的股票市场,回归的周期会不一样的,就算是相同的市场,回归的周期也是不一样的。

我们换支股票,以苏宁云商(002024)股票日K线图为例, 同样截取2005年到2015年7月的向前复权的股价数据,如下图所示。我们看到苏宁云商在2006年到2007年有一波大涨随后下跌;从2009到2010年时,第二波大涨;2013年下半年迎来第三波大涨;2014年下半年到2015年第四波大涨。从图形上可以直观看到,2015年这波涨的最急,波动率也是最大的;从现象中,我们可以判断一种趋势不管其持续的时间多长都不能永远持续下去。

02

不对称性,股价波动的幅度与速度是不一样的,回归时的幅度与速度具有随机性。对称的均值回归才是不正常的、偶然的,这一点也也可以从股票中所验证。

我们合并平安银行(000001)和苏宁云商(002024)股票日K线图为例,所下图所示。两支股票在2007年中,都赶上了大的上涨行情,曲线基本吻合。到2008年2支股票都遇到了大跌,但波动率和速度都是不一样的,随后在2010年到2012年出现了完成不一样的走势,无规律可寻,体现了均值回归时的随机性和不对称性。

03

政府行为,股票收益率不会偏离价值均值时间太久,市场的内在力量会促使其向内在价值回归。市场在没有政府政策的作用下,股票价格会在市场机制下自然地向均值回归。但这并不否定政府行为对促进市场有效性的作用,因为市场偏离内在价值后并不等于立即就会向内在价值回归,很可能会出现持续地均值回避。政府行为会起到抑制市场调节市场的作用,是必不可少的因素之一,市场失灵也是政府参与调控的直接的结果。

对于政府政策行为,比如升准、降准、升息、降息,在股市中都会有比如明显的体现。房地产股、银行股,都会受到国家宏观调控的直接的影响。下如所示,在图中增加万科A(0000002)的股票,图中3条线分别是平安银行,万科A,苏宁云商3支股票。我们发现地产和银行的股价走势是比较相近的,而电商的走势是不太一样的。

另外,增加2种颜色的辅助线,红色为升息的时间点和利率变动值,黄色为降息的时间点和利率变动值。当2007年股市超涨的时候,国家宏观调控通过升息鼓励存款,抑制高股价;当股票超跌的时候,通过降息推动投资和消费。2015年金融改革,政府一直都在降息拉动股市。从图中,我们看到万科A和平安银行对于升息和降息的调控是比较明显的,对于苏宁云商就不是特别的明显了。

04

通过对市场的回顾,我们基本验证了均值回归的理论是和市场的行为是一致的。那么,接下来我们应该如何应用这个理论来找到投资的切入点呢?

1.2 计算原理和公式

从价值投资的角度,我们发现股价会在平均值上下波动,但如果考虑到资金的时间成本,把钱都压在股市中,等待几年的大行情,也是很不划算的。那么我们就需要对价值均值进行重新定义,以20日均值来代替长期均值,找到短周期的一种投资方法。

计算原理:取日K线,以N日均线做为均值回归的短期均衡价格水平(均值),计算股价到均值的差值,求出差值的N日的平均标准差,从而判断差值的对于均值的偏离,当偏离超过2倍标准差时,我们就认为股价超涨或超跌,股价会遵循均值回归的理论,向均值不停地进行修复。

计算公式:


N日平均值     =  [T日股价 + (T-1)日股价 + ... + (T-(N-1))日股价]/N
差值          =  N日平均值 - N日股价
N日差值均值   =  [T日差值 + (T-1)日差值 + ... + (T-(N-1))日差值]/N
N日差值标准差 =  sqrt([(T日差值 - T日差值均值)^2 + ... + ((T-(N-1))日差值 - (T-(N-1))日差值均值)^2 ]/N)

如果N为20日,则


20日平均值     =  [T日股价 + (T-1)日股价 + ... + (T-19)日股价]/20

计算偏离点


T日差值 > T日差值标准差 * 2

我们以偏离点作为买入信号点,以均线和股价的下一个交点做为卖出信号点。这样我们就把均值回归的投资理论,变成了一个数学模型。

2. 均值回归模型和实现

接下来,我们利用R语言对股票数据的进行操作,来实现一个均值回归模型的实例,从而验证我的们投资理论,是否能发现赚钱的机会。

2.1 数据准备

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

我本次用到的数据是从 况客 直接导出的,况客 会提供各种类型的金融数据API,让开发者可以免费下载。当然,你也可以用quantmod包从Yahoo财经下载。

本文用到的数据,包括A股日K线(向前复权)数据,从2014年7月到2015年日7月,以CSV格式保存到本地文件stock.csv。

数据格式如下:


000001.SZ,2014-07-02,8.14,8.18,8.10,8.17,28604171
000002.SZ,2014-07-02,8.09,8.13,8.05,8.12,40633122
000004.SZ,2014-07-02,13.9,13.99,13.82,13.95,1081139
000005.SZ,2014-07-02,2.27,2.29,2.26,2.28,4157537
000006.SZ,2014-07-02,4.57,4.57,4.50,4.55,5137384
000010.SZ,2014-07-02,6.6,6.82,6.5,6.73,9909143

一共7列:

  • 第1列,股票代码,code,000001.SZ
  • 第2列,交易日期,date,2014-07-02
  • 第3列,开盘价,Open,8.14
  • 第4列,最高价,High,8.18
  • 第5列,最低价,Low,8.10
  • 第6列,收盘价,Close,8.17
  • 第7列,交易量,Volume,28604171

通过R语言加载股票数据,由于数据所有股票都是混合在一起的,而进行计算时又需要按每支票股计算,所以在数据加载时我就进行了转换,按股票代码进行分组,生成R语言的list对象,同时把每支股票的data.frame类型对象转成XTS时间序列类型对象,方便后续的数据处理。


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

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

# 加载数据
> data<-read("stock.csv")

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

# 查看数据的索引值
> head(names(data))
[1] "000001.SZ" "000002.SZ" "000004.SZ" "000005.SZ" "000006.SZ" "000007.SZ"

# 查看包括的股票数量
> length(data)
[1] 2782

# 查看股票000001.SZ
> head(data[['000001.SZ']])
               Open     High      Low    Close   Volume
2014-07-02 8.146949 8.180000 8.105636 8.171737 28604171
2014-07-03 8.171737 8.254364 8.122162 8.229576 44690486
2014-07-04 8.237838 8.270889 8.146949 8.188263 34231126
2014-07-07 8.188263 8.204788 8.097374 8.146949 34306164
2014-07-08 8.130424 8.204788 8.072586 8.204788 34608702
2014-07-09 8.196525 8.196525 7.915596 7.973434 58789114

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

2.2 均值回归模型
为了能拉近我们对市场的了解,我们取从2015年1月1日开始的数据,来创建均值回归模型。以平安银行(000001)的为例,画出平安银行的2015年以来的日K线和均线。


# 获得时间范围
> dateArea<-function(sDate=Sys.Date()-365,eDate= Sys.Date(),before=0){  #开始日期,结束日期,提单开始时
+     if(class(sDate)=='character') sDate=as.Date(sDate)
+     if(class(eDate)=='character') eDate=as.Date(eDate)  
+     return(paste(sDate-before,eDate,sep="/"))
+ }
 
# 计算移动平均线
> ma<-function(cdata,mas=c(5,20,60)){
+     if(nrow(cdata)<=max(mas)) return(NULL)
+     ldata<-cdata
+     for(m in mas){
+         ldata<-merge(ldata,SMA(cdata,m))
+     }
+     names(ldata)<-c('Value',paste('ma',mas,sep=''))
+     return(ldata)
+ }

# 日K线和均线
> title<-'000001.SZ'
> SZ000011<-data[[title]]                             # 获得股票数据
> sDate<-as.Date("2015-01-01")                        # 开始日期
> eDate<-as.Date("2015-07-10")                        # 结束日期
> cdata<-SZ000011[dateArea(sDate,eDate,360)]$Close    # 获得收盘价
> ldata<-ma(cdata,c(5,20,60))                         # 选择移动平均指标

# 打印移动平均指标
> tail(ldata)
           Value    ma5    ma20     ma60
2015-07-03 13.07 13.768 15.2545 15.84355
2015-07-06 13.88 13.832 15.1335 15.82700
2015-07-07 14.65 13.854 15.0015 15.79850
2015-07-08 13.19 13.708 14.8120 15.74267
2015-07-09 14.26 13.810 14.6910 15.70867
2015-07-10 14.86 14.168 14.6100 15.67883

我们设置3条移动平均线,分别是5日平均线,20日平均线,60日平均线,当然也可以按照自己的个性要求设置符合自己的周期。画出日K线和均线图。


> drawLine<-function(ldata,titie="Stock_MA",sDate=min(index(ldata)),eDate=max(index(ldata)),breaks="1 year",avg=FALSE,out=FALSE){
+     if(sDate<min(index(ldata))) sDate=min(index(ldata))
+     if(eDate>max(index(ldata))) eDate=max(index(ldata))  
+     ldata<-na.omit(ldata)
+     
+     g<-ggplot(aes(x=Index, y=Value),data=fortify(ldata[,1],melt=TRUE))
+     g<-g+geom_line()
+     g<-g+geom_line(aes(colour=Series),data=fortify(ldata[,-1],melt=TRUE))
+ 
+     if(avg){
+         meanVal<<-round(mean(ldata[dateArea(sDate,eDate)]$Value),2) # 均值
+         g<-g+geom_hline(aes(yintercept=meanVal),color="red",alpha=0.8,size=1,linetype="dashed")
+         g<-g+geom_text(aes(x=sDate, y=meanVal,label=meanVal),color="red",vjust=-0.4)
+     }
+     g<-g+scale_x_date(labels=date_format("%Y-%m"),breaks=date_breaks(breaks),limits = c(sDate,eDate))
+     g<-g+ylim(min(ldata$Value), max(ldata$Value))
+     g<-g+xlab("") + ylab("Price")+ggtitle(title)
+     g
+ }

> drawLine(ldata,title,sDate,eDate,'1 month',TRUE)    # 画图

05

如图所示,60日的移动平均线是最平滑的,5日的移动平均线是波动最大的。5日平均线和股价的交叉,明显多于60日平均线和股价的交叉。那么可以说在相同的时间周期内,短周期的移动平均线,比长周期的移动平均线更具有均值回归的特点。

我们分别计算不同周期的,股价与移动平均线的差值的平均标准差。


> getMaSd<-function(ldata,mas=20,sDate,eDate){}) # ...代码省略

# 5日平均线的差值、平均标准差
> ldata5<-getMaSd(ldata,5,sDate,eDate)
> head(ldata5)
              Value      ma5        dif        sd  rate
2015-01-05 13.23673 12.78724 -0.4494869 0.1613198 -2.79
2015-01-06 13.03842 12.89961 -0.1388121 0.1909328 -0.73
2015-01-07 12.79055 12.99215  0.2016081 0.3169068  0.64
2015-01-08 12.36089 12.90292  0.5420283 0.4472248  1.21
2015-01-09 12.46004 12.77733  0.3172848 0.3910700  0.81
2015-01-12 12.20390 12.57076  0.3668606 0.2533165  1.45


# 20日平均线的差值、平均标准差
> ldata20<-getMaSd(ldata,20,sDate,eDate)
> head(ldata20)
              Value     ma20         dif        sd  rate
2015-01-05 13.23673 12.18613 -1.05059293 0.6556366 -1.60
2015-01-06 13.03842 12.23778 -0.80064848 0.6021093 -1.33
2015-01-07 12.79055 12.24810 -0.54244141 0.4754686 -1.14
2015-01-08 12.36089 12.29975 -0.06114343 0.5130410 -0.12
2015-01-09 12.46004 12.33651 -0.12352626 0.5150453 -0.24
2015-01-12 12.20390 12.37163  0.16773131 0.5531618  0.30


# 60日平均线的差值、平均标准差
> ldata60<-getMaSd(ldata,60,sDate,eDate)
> head(ldata60)
              Value     ma60       dif       sd  rate
2015-01-05 13.23673 10.06939 -3.167340 1.264792 -2.50
2015-01-06 13.03842 10.14678 -2.891644 1.271689 -2.27
2015-01-07 12.79055 10.22087 -2.569677 1.269302 -2.02
2015-01-08 12.36089 10.28752 -2.073368 1.258813 -1.65
2015-01-09 12.46004 10.35527 -2.104766 1.247967 -1.69
2015-01-12 12.20390 10.41821 -1.785691 1.233989 -1.45

5日的平均线的差值和平均标准差是最小的,而60日的平均线的差值和平均标准差是最大的。如果我们以5日移动平均线做为均值时,会频繁进行交易,但每次收益都很小,可能都不够手续费的成本;另一方面,如果我们以60日移动平均线做为均值时,交易次数会较少,但可能会出现股票成形趋势性上涨或下跌,长时间不能回归的情况,可能会造成现金头寸的紧张。综合上面的2种情况,我们可以选择20日均线作为均值的标的。

根据模型的计算公式,当差值超过2倍的平均标准差时,我们认为股价出现了偏离,以偏离点做为模型的买入信号,当均线和股价再次相交时做为卖出信号。

上一步,我们已经计算出了偏离值,并保存在rate列中。下面我们要找到大于2倍标准化差的点,并画图。


# 差值和平均标准差,大于2倍平均标准差的点
> buyPoint<-function(ldata,x=2,dir=2){})     # ...代码省略

# 画交易信号点
> drawPoint<-function(ldata,pdata,titie,sDate,eDate,breaks="1 year"){
+     ldata<-na.omit(ldata)
+     g<-ggplot(aes(x=Index, y=Value),data=fortify(ldata[,1],melt=TRUE))
+     g<-g+geom_line()
+     g<-g+geom_line(aes(colour=Series),data=fortify(ldata[,-1],melt=TRUE))
+     
+     if(is.data.frame(pdata)){
+         g<-g+geom_point(aes(x=Index,y=Value,colour=op),data=pdata,size=4)
+     }else{
+         g<-g+geom_point(aes(x=Index,y=Value,colour=Series),data=na.omit(fortify(pdata,melt=TRUE)),size=4)  
+     }
+     g<-g+scale_x_date(labels=date_format("%Y-%m"),breaks=date_breaks(breaks),limits = c(sDate,eDate))
+     g<-g+xlab("") + ylab("Price")+ggtitle(title)
+     g
+ }
 
> buydata<-buyPoint(ldata20,2,2)                                       # 多空信号点
> drawPoint(ldata20[,c(1,2)],buydata$Value,title,sDate,eDate,'1 month')  # 画图

06

图中蓝色的点就是买入的信号点,由于股票我们只能进行单向交易,即低买高卖,并不能直接做空,所以我们要过滤股价高于移动平均线的点,只留下股价低于移动平均线的点,就是我们的买入信号点。

画出买入信号点,只保留股价低于移动平均线的点。


> buydata<-buyPoint(ldata20,2,1)        # 做多信号点
> drawPoint(ldata20[,c(1,2)],buydata$Value,title,sDate,eDate,'1 month') # 画图

07

计算卖出的信号点,当买入后,下一个股价与移动平均线的交点就是卖出的信号点,我们看一下是否可以赚到钱?!


# 计算卖出的信号点
> sellPoint<-function(ldata,buydata){})     # ...代码省略
> selldata<-sellPoint(ldata20,buydata)

# 买出信号
> selldata
           Value  ma20   dif        sd  rate
2015-07-10 14.86 14.61 -0.25 0.7384824 -0.34

我们把买入信号和卖出信号,合并到一张图上显示,如图所示。


> bsdata<-merge(buydata$Value,selldata$Value)
> names(bsdata)<-c("buy","sell")
> drawPoint(ldata20[,c(1,2)],bsdata,title,sDate,eDate,'1 month') #画图

08

从图上看,我们在绿色点位置进行买入,而在蓝色点位置进行卖出,确实是赚钱的。那么究竟赚了多少钱呢?我们还需要精确的计算出来。


# 合并交易信号
> signal<-function(buy, sell){})    # ...代码省略

# 交易信号数据
> sdata<-signal(buydata,selldata)                                   
> sdata
           Value    ma20     dif        sd  rate op
2015-06-19 14.63 16.0965  1.4665 0.6620157  2.22  B
2015-06-26 13.77 15.7720  2.0020 0.8271793  2.42  B
2015-06-29 13.56 15.6840  2.1240 0.9271735  2.29  B
2015-07-03 13.07 15.2545  2.1845 1.0434926  2.09  B
2015-07-10 14.86 14.6100 -0.2500 0.7384824 -0.34  S

利用交易信号数据,进行模拟交易。我们设定交易参数和规则:

  • 以10万元人民币为本金
  • 买入信号出现时,以收盘价买入,每次买入价值1万元的股票。如果连续出现买入信号,则一直买入。若现金不足1万元时,则跳过买入信号。
  • 卖出信号出现时,以收盘价卖出,一次性平仓信号对应的股票。
  • 手续费为0元

# 模拟交易
> trade<-function(sdata,capital=100000,fixMoney=10000){})    # ...代码省略

# 交易结果
> result<-trade(sdata,100000,10000)  

来看一下,每笔交易的明细。


> result$ticks
           Value    ma20     dif        sd  rate op      cash amount     asset     diff
2015-06-19 14.63 16.0965  1.4665 0.6620157  2.22  B  90007.71    683 100000.00     0.00
2015-06-26 13.77 15.7720  2.0020 0.8271793  2.42  B  80010.69   1409  99412.62  -587.38
2015-06-29 13.56 15.6840  2.1240 0.9271735  2.29  B  70016.97   2146  99116.73  -295.89
2015-07-03 13.07 15.2545  2.1845 1.0434926  2.09  B  60018.42   2911  98065.19 -1051.54
2015-07-10 14.86 14.6100 -0.2500 0.7384824 -0.34  S 103275.88      0 103275.88  5210.69

一共发生了5笔交易,其中4笔买入,1笔卖出。最后,资金剩余103275.88元,赚了3275.88元,收益率3.275%。

在卖出时,赚钱的交易有1笔。


> result$rise
           Value  ma20   dif        sd  rate op     cash amount    asset    diff
2015-07-10 14.86 14.61 -0.25 0.7384824 -0.34  S 103275.9      0 103275.9 5210.69

在卖出时,赔钱的交易,没有发生。


> result$fall
 [1] Value  ma20   dif    sd     rate   op     cash   amount asset  diff  
<0 行> (或0-长度的row.names)

接下来,我们再对比一下,资产净值和股价。


# 资产净值曲线
> drawAsset<-function(ldata,adata,sDate=FALSE,capital=100000){
+     if(!sDate) sDate<-index(ldata)[1]
+     adata<-rbind(adata,as.xts(capital,as.Date(sDate)))
+     
+     g<-ggplot(aes(x=Index, y=Value),data=fortify(ldata[,1],melt=TRUE))
+     g<-g+geom_line()
+     g<-g+geom_line(aes(x=as.Date(Index), y=Value,colour=Series),data=fortify(adata,melt=TRUE))
+     g<-g+facet_grid(Series ~ .,scales = "free_y")
+     g<-g+scale_y_continuous(labels=dollar_format(prefix = "¥"))
+     g<-g+scale_x_date(labels=date_format("%Y-%m"),breaks=date_breaks("2 months"),limits = c(sDate,eDate))
+     g<-g+xlab("") + ylab("Price")+ggtitle(title)
+     g
+ }

> drawAsset(ldata20,as.xts(result$ticks['asset']))  # 资产净值曲线

09

刚才我们是对一支股票进行了测试,发现是有机会的,那么我再换另外一支股票,看一下是否用同样的效果呢?我们把刚才数据操作的过程,封装到统一的quick函数,就可以快速验证均值回归在其他股票的表现情况了。


> quick<-function(title,sDate,eDate){}  # ...代码省略

我们用乐视网(300104)试一下,看看有没有赚钱的机会!!


> title<-"300104.SZ"
> sDate<-as.Date("2015-01-01") #开始日期
> eDate<-as.Date("2015-07-10") #结束日期

> quick(title,sDate,eDate)
$ticks
           Value    ma20     dif       sd  rate op      cash amount     asset     diff
2015-06-19 55.04 69.9095 14.8695 5.347756  2.78  B  90037.76    181 100000.00     0.00
2015-06-23 54.30 68.8075 14.5075 5.477894  2.65  B  80046.56    365  99866.06  -133.94
2015-06-24 56.21 67.8735 11.6635 5.404922  2.16  B  70097.39    542 100563.21   697.15
2015-06-25 51.80 66.8775 15.0775 5.770806  2.61  B  60099.99    735  98172.99 -2390.22
2015-06-26 46.79 65.9830 19.1930 6.580622  2.92  B  50133.72    948  94490.64 -3682.35
2015-06-29 47.05 64.9445 17.8945 7.096230  2.52  B  40159.12   1160  94737.12   246.48
2015-07-07 47.86 58.8150 10.9550 5.401247  2.03  B  30204.24   1368  95676.72   939.60
2015-07-10 57.92 57.3520 -0.5680 5.625309 -0.10  S 109438.80      0 109438.80 13762.08

$rise
           Value   ma20    dif       sd rate op     cash amount    asset     diff
2015-07-10 57.92 57.352 -0.568 5.625309 -0.1  S 109438.8      0 109438.8 13762.08

$fall
 [1] Value  ma20   dif    sd     rate   op     cash   amount asset  diff  
<0 行> (或0-长度的row.names)

从数据结果看,我们又赚到了。一共发生了8笔交易,其中7笔买入,1笔卖出。最后,资金剩余109438.80元,赚了9438.80元,收益率9.43%。

画出交易信号图


> title<-"300104.SZ"
> sDate<-as.Date("2015-01-01") #开始日期
> eDate<-as.Date("2015-07-10") #结束日期

> stock<-data[[title]]
> cdata<-stock[dateArea(sDate,eDate,360)]$Close
> ldata<-ma(cdata,c(20))
> ldata<-getMaSd(ldata,20,sDate,eDate)
> buydata<-buyPoint(ldata,2,1)  
> selldata<-sellPoint(ldata,buydata)
> bsdata<-merge(buydata$Value,selldata$Value)
> drawPoint(ldata[,c(1,2)],bsdata,title,sDate,eDate,'1 month') #画图

10

在恐慌的6月份,当别人都被套牢30%以上的情况下,我们还朿9%正收益,那么应该是多么舒心的一件事情啊!!

3. 量化选股

上文中,我们用2支股票进行了测试,发现均值回归模型是适合于股票交易的。如果我们利用模型对全市场的股票进行扫描,应用会产生更多的交易信号,找到更多的投资机会,这样我们就能如何能获得更大的收益。

那么,接下来我们就根据均值回归的理论进行量化选股。

根据我们之前的经验,当股价与平均标准差的偏离越大,有可能带来的收益就越大。那么通过量化的手段,在整个的市场2700多支股票中,把每天偏离最大股票的找出来进行交易,就可以有效地分配我们的资金,进行更有效的投资。我们要试一下,市场是否是和我们的思路是一致的。

对全市场股票进行扫描,首先计算差值、平均值和平均标准差。


> sDate<-as.Date("2015-01-01")                # 开始日期
> eDate<-as.Date("2015-07-10")                # 结束日期

# 计算差值、平均值和平均标准差
> data0<-lapply(data,function(stock){})       # 代码省略

# 去掉空数据
> data0<-data0[!sapply(data0, is.null)]      

# 全市场股票
> length(data)
[1] 2782

# 有效的股票
> length(data0)
[1] 2697

# 查看第1支股票
> head(data0[[1]])
              Value     ma20         dif        sd  rate
2015-01-05 13.23673 12.18613 -1.05059293 0.6556366 -1.60
2015-01-06 13.03842 12.23778 -0.80064848 0.6021093 -1.33
2015-01-07 12.79055 12.24810 -0.54244141 0.4754686 -1.14
2015-01-08 12.36089 12.29975 -0.06114343 0.5130410 -0.12
2015-01-09 12.46004 12.33651 -0.12352626 0.5150453 -0.24
2015-01-12 12.20390 12.37163  0.16773131 0.5531618  0.30

第一次扫描后,有2697支股票是符合条件的,有85支股票由于数据样本不足被排除。

接下来,继续对2697支股票进行筛选,找到符合要求的买入信号点。


# 计算买入信号
> buys<-lapply(data0,function(stock){})  # ...代码省略 

# 去掉空数据
> buys<-buys[!sapply(buys, is.null)] 

# 查看有买入信号的股票
> length(buys)
[1] 1819

# 查看买入信号
> head(buys)
$`000001.SZ`
           Value    ma20    dif        sd rate
2015-06-19 14.63 16.0965 1.4665 0.6620157 2.22
2015-06-26 13.77 15.7720 2.0020 0.8271793 2.42
2015-06-29 13.56 15.6840 2.1240 0.9271735 2.29
2015-07-03 13.07 15.2545 2.1845 1.0434926 2.09

$`000002.SZ`
           Value   ma20   dif        sd rate
2015-03-05 11.90 12.568 0.668 0.2644101 2.53
2015-03-06 11.94 12.509 0.569 0.2674732 2.13

$`000004.SZ`
           Value    ma20     dif        sd rate
2015-01-05 15.69 17.7210  2.0310 0.7395717 2.75
2015-07-06 26.03 39.1540 13.1240 6.3898795 2.05
2015-07-07 23.43 38.2025 14.7725 6.9421723 2.13
2015-07-08 22.22 37.2635 15.0435 7.4287088 2.03

$`000005.SZ`
           Value    ma20    dif       sd rate
2015-07-06  6.02 10.9600 4.9400 2.381665 2.07
2015-07-07  5.42 10.5655 5.1455 2.333008 2.21

$`000006.SZ`
              Value     ma20       dif      sd rate
2015-01-19 5.829283 6.519462 0.6901792 0.26929 2.56

$`000007.SZ`
           Value    ma20    dif        sd rate
2015-02-06 12.47 14.4200 1.9500 0.6182860 3.15
2015-02-09 12.52 14.3270 1.8070 0.7440473 2.43
2015-02-10 12.10 14.1845 2.0845 0.8484250 2.46

通过计算发现,有1819支股票,在这半年中产生过买入信号。每支股票产生的买入信号的时间和频率都是不同,这样我们就可以把钱分散投资到不同的股票上,同时分散风险。如果交易信号同一天出现在多支的股票上,而我们资金有限,又想让收益最大化,那么我们可以选择偏离值最大的股票进行交易。

接下来,我们用程序找到每日偏离最大的股票。


# 合并数据,从list转型到data.frame
buydf<-ldply(buys,function(e){})    # ...代码省略

# 选出同一日rate最大的股票,做为买入信号
buydatas<-ddply(buydf, .(date), function(row){}) # ...代码省略

# 查看买入信号
> nrow(buydatas)
[1] 81

# 查看买入信号细节
> head(buydatas)
         .id       date      Value       ma20        dif         sd rate
1  002551.SZ 2015-01-05  16.573846  19.565446  2.9916000 0.74591596 4.01
2  002450.SZ 2015-01-06  18.548809  19.766636  1.2178275 0.34008453 3.58
3  300143.SZ 2015-01-07  11.480000  12.603000  1.1230000 0.32028018 3.51
4  300335.SZ 2015-01-08  12.113677  13.139601  1.0259238 0.21760484 4.71
5  300335.SZ 2015-01-09  12.243288  13.043888  0.8005994 0.22940845 3.49
6  300335.SZ 2015-01-12  11.994036  12.941694  0.9476584 0.23168313 4.09

最后,我们选出81个买入信号点,基本上每个交易日都是买入信号。有了买入信号,继续找到卖出信号。


# 卖出信号
> selldatas<-data.frame()     # ...代码省略

# 卖出信号去重
> selldatas<-unique(selldatas)  
> nrow(selldatas)
[1] 33

# 查看买出信号
> head(selldatas)
                Value      ma20         dif        sd  rate       .id       date op
2015-01-12  19.232308 18.848908 -0.38340000 0.9051374 -0.42 002551.SZ 2015-01-12  S
2015-01-08  19.814257 19.729006 -0.08525126 0.3782955 -0.23 002450.SZ 2015-01-08  S
2015-01-28  11.210000 11.019500 -0.19050000 0.7781848 -0.24 300143.SZ 2015-01-28  S
2015-01-21  13.190448 12.899321 -0.29112706 0.3871871 -0.75 300335.SZ 2015-01-21  S
2015-01-213  7.140000  6.989500 -0.15050000 0.2007652 -0.75 002505.SZ 2015-01-21  S
2015-01-22   5.561561  5.490668 -0.07089242 0.2127939 -0.33 600077.SH 2015-01-22  S

通过计算,一共有33个买出信号点。最后,合并买入信号和卖出信号,并计算收益。


> buydatas$op<-'B'                              # 买入标志
> selldatas$op<-'S'                             # 卖出标志
> sdatas<-rbind(buydatas,selldatas)             # 合并数据
> row.names(sdatas)<-1:nrow(sdatas)             # 重设行号
> sdatas<-sdatas[order(sdatas$.id),]            # 按股票代码排序

# 查看合并的信号
> head(sdatas)
          .id       date Value     ma20       dif         sd  rate op
36  000002.SZ 2015-03-05 11.90 12.56800  0.668000 0.26441011  2.53  B
100 000002.SZ 2015-03-16 12.49 12.38050 -0.109500 0.23702768 -0.46  S
58  000553.SZ 2015-05-06 14.35 15.50882  1.158824 0.38429912  3.02  B
110 000553.SZ 2015-05-21 16.57 15.18903 -1.380972 0.55647152 -2.48  S
26  000725.SZ 2015-02-09  2.80  3.11400  0.314000 0.07934585  3.96  B
94  000725.SZ 2015-02-16  3.09  3.06500 -0.025000 0.08182388 -0.31  S

最后,按照股票进行分组,分别计算个股的收益。


# 计算个股的收益
> slist<-split(sdatas[-1],sdatas$.id)      # 按股票代码分组
> results<-lapply(slist,trade)

# 查看信号的股票
> names(results)
 [1] "000002.SZ" "000553.SZ" "000725.SZ" "000786.SZ" "000826.SZ" "002240.SZ" "002450.SZ"
 [8] "002496.SZ" "002505.SZ" "002544.SZ" "002551.SZ" "002646.SZ" "002652.SZ" "300143.SZ"
[15] "300335.SZ" "300359.SZ" "300380.SZ" "300397.SZ" "300439.SZ" "300440.SZ" "300444.SZ"
[22] "600030.SH" "600038.SH" "600077.SH" "600168.SH" "600199.SH" "600213.SH" "600375.SH"
[29] "600490.SH" "600536.SH" "600656.SH" "600733.SH" "600890.SH" "601179.SH" "601186.SH"
[36] "601628.SH" "601633.SH" "601939.SH" "603019.SH"

我们查看万科A(000002)的股票。


> results[['000002.SZ']]$ticks
          date Value    ma20     dif        sd  rate op     cash amount    asset  diff
36  2015-03-05 11.90 12.5680  0.6680 0.2644101  2.53  B  90004.0    840 100000.0   0.0
100 2015-03-16 12.49 12.3805 -0.1095 0.2370277 -0.46  S 100495.6      0 100495.6 495.6

通过优化的规则设计,一共有2笔交易,赚了495元。如要我们没有进行算法优化,一直交易万科A,那么会发生3笔交易,我们可以赚955.95元。


> quick('000002.SZ',sDate,eDate)$ticks
           Value    ma20     dif        sd  rate op      cash amount    asset   diff
2015-03-05 11.90 12.5680  0.6680 0.2644101  2.53  B  90004.00    840 100000.0   0.00
2015-03-06 11.94 12.5090  0.5690 0.2674732  2.13  B  80010.22   1677 100033.6  33.60
2015-03-16 12.49 12.3805 -0.1095 0.2370277 -0.46  S 100955.95      0 100955.9 922.35

本文到此就要结束了!但其实还有很多的事情要做,比如对模型参数的优化,用10日均线代替20日均线,用3倍标准差偏移代替2倍标准差偏移,对样本进行正态分布的检验,结合其他趋势类模型共同产生信号等,这些就不是一篇文章可以解决的事情了。大家可以况客金融平台的网站上,发现更多不一样的策略。

本文从均值回归的理论的介绍开始,到市场特征检验,再到数学公式,R语言建模,历史数据回测,最后找到投资机会,是一套完整的从理论到实践的学习方法。虽然困难重重,但做为有理想的极客,我们是有能力来克服这些困难的。

本文同时用到了计算机、金融、数学、统计等多学科知识的结合,我认为这是技术复合人才未来的发展方向。如果说过去10年是房地产的黄金10年,那么未来的10年将是金融的黄金10年。当我们IT人掌握了足够的金融知识,一定会有能力去金融市场抢钱的。

抓住机会!!程序员,加油!

######################################################
看文字不过瘾,作者视频讲解,请访问网站:http://onbook.me/video
######################################################

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

打赏作者

R语言中的遗传算法

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

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

ga-r

前言

人类总是在生活中摸索规律,把规律总结为经验,再把经验传给后人,让后人发现更多的规规律,每一次知识的传递都是一次进化的过程,最终会形成了人类的智慧。自然界规律,让人类适者生存地活了下来,聪明的科学家又把生物进化的规律,总结成遗传算法,扩展到了更广的领域中。

本文将带你走进遗传算法的世界。

目录

  1. 遗传算法介绍
  2. 遗传算法原理
  3. 遗传算法R语言实现

1. 遗传算法介绍

遗传算法是一种解决最优化的搜索算法,是进化算法的一种。进化算法最初借鉴了达尔文的进化论和孟德尔的遗传学说,从生物进化的一些现象发展起来,这些现象包括遗传、基因突变、自然选择和杂交等。遗传算法通过模仿自然界生物进化机制,发展出了随机全局搜索和优化的方法。遗传算法其本质是一种高效、并行、全局搜索的方法,它能在搜索过程中自动获取和积累有关搜索空间的知识,并自适应的控制搜索过程,计算出全局最优解。

遗传算法的操作使用适者生存的原则,在潜在的种群中逐次产生一个近似最优解的方案,在每一代中,根据个体在问题域中的适应度值和从自然遗传学中借鉴来的再造方法进行个体选择,产生一个新的近似解。这个过程会导致种群中个体的进化,得到的新个体比原来个体更能适应环境,就像自然界中的改造一样。

如果从生物进化的角度,我们可以这样理解。在一个种群中,个体数量已经有一定规模,为了进化发展,通过选择和繁殖产生下一代的个体,其中繁殖过程包括交配和突变。根据适者生存的原则,选择过程会根据新个体的适应度进行保留或淘汰,但也不是完全以适应度高低作为导向,如果单纯选择适应度高的个体可能会产生局部最优的种群,而非全局最优,这个种群将不会再进化,称为早熟。之后,通过繁殖过程,让个体两两交配产生下一代新个体,上一代个体中优秀的基因会保留给下一代,而劣制的基因将被个体另一半的基因所代替。最后,通过小概率事件发生基因突变,通过突变产生新的下一代个体,实现种群的变异进化。

经过这一系列的选择、交配和突变的过程,产生的新一代个体将不同于初始的一代,并一代一代向增加整体适应度的方向发展,因为最好的个体总是更多的被选择去产生下一代,而适应度低的个体逐渐被淘汰掉。这样的过程不断的重复:每个个体被评价,计算出适应度,两个个体交配,然后突变,产生第三代。周而复始,直到终止条件满足为止。

遗传算法需要注意的问题:

  • 遗传算法在适应度函数选择不当的情况下有可能收敛于局部最优,而不能达到全局最优。
  • 初始种群的数量很重要,如果初始种群数量过多,算法会占用大量系统资源;如果初始种群数量过少,算法很可能忽略掉最优解。
  • 对于每个解,一般根据实际情况进行编码,这样有利于编写变异函数和适应度函数。
  • 在编码过的遗传算法中,每次变异的编码长度也影响到遗传算法的效率。如果变异代码长度过长,变异的多样性会受到限制;如果变异代码过短,变异的效率会非常低下,选择适当的变异长度是提高效率的关键。
  • 变异率是一个重要的参数。
  • 对于动态数据,用遗传算法求最优解比较困难,因为染色体种群很可能过早地收敛,而对以后变化了的数据不再产生变化。对于这个问题,研究者提出了一些方法增加基因的多样性,从而防止过早的收敛。其中一种是所谓触發式超级变异,就是当染色体群体的质量下降(彼此的区别减少)时增加变异概率;另一种叫随机外来染色体,是偶尔加入一些全新的随机生成的染色体个体,从而增加染色体多样性。
  • 选择过程很重要,但交叉和变异的重要性存在争议。一种观点认为交叉比变异更重要,因为变异仅仅是保证不丢失某些可能的解;而另一种观点则认为交叉过程的作用只不过是在种群中推广变异过程所造成的更新,对于初期的种群来说,交叉几乎等效于一个非常大的变异率,而这么大的变异很可能影响进化过程。
  • 遗传算法很快就能找到良好的解,即使是在很复杂的解空间中。
  • 遗传算法并不一定总是最好的优化策略,优化问题要具体情况具体分析。所以在使用遗传算法的同时,也可以尝试其他算法,互相补充,甚至根本不用遗传算法。
  • 遗传算法不能解决那些“大海捞针”的问题,所谓“大海捞针”问题就是没有一个确切的适应度函数表征个体好坏的问题,使得算法的进化失去导向。
  • 对于任何一个具体的优化问题,调节遗传算法的参数可能会有利于更好的更快的收敛,这些参数包括个体数目、交叉率和变异率。例如太大的变异率会导致丢失最优解,而过小的变异率会导致算法过早的收敛于局部最优点。对于这些参数的选择,现在还没有实用的上下限。
  • 适应度函数对于算法的速度和效果也很重要。

遗传算法的应用领域包括计算机自动设计、生产调度、电路设计、游戏设计、机器人学习、模糊控制、时间表安排,神经网络训练等。然而,我准备把遗传算法到金融领域,比如回测系统的参数寻优方案,我会在以后的文章中,介绍有关金融解决方案。

2. 遗传算法原理

在遗传算法里,优化问题的解是被称为个体,它表示为一个变量序列,叫做染色体或者基因串。染色体一般被表达为简单的字符串或数字串,也有其他表示法,这一过程称为编码。首先要创建种群,算法随机生成一定数量的个体,有时候也可以人工干预这个过程进行,以提高初始种群的质量。在每一代中,每一个个体都被评价,并通过计算适应度函数得到一个适应度数值。种群中的个体被按照适应度排序,适应度高的在前面。

接下来,是产生下一代个体的种群,通过选择过程和繁殖过程完成。

选择过程,是根据新个体的适应度进行的,但同时并不意味着完全的以适应度高低作为导向,因为单纯选择适应度高的个体将可能导致算法快速收敛到局部最优解而非全局最优解,我们称之为早熟。作为折中,遗传算法依据原则:适应度越高,被选择的机会越高,而适应度低的,被选择的机会就低。初始的数据可以通过这样的选择过程组成一个相对优化的群体。

繁殖过程,表示被选择的个体进入交配过程,包括交配(crossover)和突变(mutation),交配对应算法中的交叉操作。一般的遗传算法都有一个交配概率,范围一般是0.6~1,这个交配概率反映两个被选中的个体进行交配的概率。

例如,交配概率为0.8,则80%的“夫妻”个体会生育后代。每两个个体通过交配产生两个新个体,代替原来的“老”个体,而不交配的个体则保持不变。交配过程,父母的染色体相互交換,从而产生两个新的染色体,第一个个体前半段是父亲的染色体,后半段是母亲的,第二个个体则正好相反。不过这里指的半段並不是真正的一半,这个位置叫做交配点,也是随机产生的,可以是染色体的任意位置。

突变过程,表示通过突变产生新的下一代个体。一般遗传算法都有一个固定的突变常数,又称为变异概率,通常是0.1或者更小,这代表变异发生的概率。根据这个概率,新个体的染色体随机的突变,通常就是改变染色体的一个字节(0变到1,或者1变到0)。

遗传算法实现将不断的重复这个过程:每个个体被评价,计算出适应度,两个个体交配,然后突变,产生下一代,直到终止条件满足为止。一般终止条件有以下几种:

  • 进化次数限制
  • 计算耗费的资源限制,如计算时间、计算占用的CPU,内存等
  • 个体已经满足最优值的条件,即最优值已经找到
  • 当适应度已经达到饱和,继续进化不会产生适应度更好的个体
  • 人为干预

算法实现原理:

ga

  • 1. 创建初始种群
  • 2. 循环:产生下一代
  • 3. 评价种群中的个体适应度
  • 4. 定义选择的适应度函数
  • 5. 改变该种群(交配和变异)
  • 6. 返回第二步
  • 7. 满足终止条件结束

3. 遗传算法R语言实现

本节的系统环境

  • Win7 64bit
  • R: 3.1.1 x86_64-w64-mingw32/x64 (64-bit)

一个典型的遗传算法要求:一个基因表示的求解域, 一个适应度函数来评价解决方案。

遗传算法的参数通常包括以下几个:

  • 种群规模(Population),即种群中染色体个体的数目。
  • 染色体的基因个数(Size),即变量的数目。
  • 交配概率(Crossover),用于控制交叉计算的使用频率。交叉操作可以加快收敛,使解达到最有希望的最优解区域,因此一般取较大的交叉概率,但交叉概率太高也可能导致过早收敛。
  • 变异概率(Mutation),用于控制变异计算的使用频率,决定了遗传算法的局部搜索能力。
  • 中止条件(Termination),结束的标志。

在R语言中,有一些现成的第三方包已经实现的遗传算法,我们可以直接进行使用。

  • mcga包,多变量的遗传算法,用于求解多维函数的最小值。
  • genalg包,多变量的遗传算法,用于求解多维函数的最小值。
  • rgenoud包,复杂的遗传算法,将遗传算法和衍生的拟牛顿算法结合起来,可以求解复杂函数的最优化化问题。
  • gafit包,利用遗传算法求解一维函数的最小值。不支持R 3.1.1的版本。
  • GALGO包,利用遗传算法求解多维函数的最优化解。不支持R 3.1.1的版本。

本文将介绍mcga包和genalg包的遗传算法的使用。

3.1 mcga包

我们使用mcga包的mcga()函数,可以实现多变量的遗传算法。

mcga包是一个遗传算法快速的工具包,主要解决实值优化的问题。它使用的变量值表示基因序列,而不是字节码,因此不需要编解码的处理。mcga实现了遗传算法的交配和突变的操作,并且可以进行大范围和高精度的搜索空间的计算,算法的主要缺点是使用了256位的一元字母表。

首先,安装mcga包。


> install.packages("mcga")
> library(mcga)

查看一下mcga()函数的定义。


> mcga
function (popsize, chsize, crossprob = 1, mutateprob = 0.01, elitism = 1, minval, maxval, maxiter = 10, evalFunc) 

参数说明:

  • popsize,个体数量,即染色体数目
  • chsize,基因数量,限参数的数量
  • crossprob,交配概率,默认为1.0
  • mutateprob,突变概率,默认为0.01
  • elitism,精英数量,直接复制到下一代的染色体数目,默认为1
  • minval,随机生成初始种群的下边界值
  • maxval,随机生成初始种群的上边界值
  • maxiter,繁殖次数,即循环次数,默认为10
  • evalFunc,适应度函数,用于给个体进行评价

接下来,我们给定一个优化的问题,通过mcga()函数,计算最优化的解。

题目1:设fx=(x1-5)^2 + (x2-55)^2 +(x3-555)^2 +(x4-5555)^2 +(x5-55555)^2,计算fx的最小值,其中x1,x2,x3,x4,x5为5个不同的变量。

从直观上看,如果想得到fx的最小值,其实当x1=5,x2=55,x3=555,x4=5555,x5=55555时,fx=0为最小值。如果使用穷举法,通过循环的方法找到这5个变量,估计会很费时的,我就不做测试了。下面我们看一下遗传算法的运行情况。


# 定义适应度函数
> f<-function(x){} # 代码省略

# 运行遗传算法
> m <- mcga(  popsize=200, 
+             chsize=5, 
+             minval=0.0, 
+             maxval=999999, 
+             maxiter=2500, 
+             crossprob=1.0, 
+             mutateprob=0.01, 
+             evalFunc=f)

# 最优化的个体结果
> print(m$population[1,])
[1]  5.000317    54.997099   554.999873  5555.003120 55554.218695

# 执行时间
> m$costs[1]
[1] 3.6104556

我们得到的最优化的结果为x1=5.000317, x2=54.997099, x3=554.999873, x4=5555.003120, x5=55554.218695,和我们预期的结果非常接近,而且耗时只有3.6秒。这个结果是非常令人满意地,不是么!如果使用穷举法,时间复杂度为O(n^5),估计没有5分钟肯定算不出来。

当然,算法执行时间和精度,都是通过参数进行配置的。如果增大个体数目或循环次数,一方面会增加算法的计算时间,另一方面结果也可能变得更精准。所以,在实际的使用过程中,需要根据一定的经验调整这几个参数。

3.2 genalg包

我们使用genalg包的rbga()函数,也可以实现多变量的遗传算法。

genalg包不仅实现了遗传算法,还提供了遗传算法的数据可视化,给用户更直观的角度理解算法。默认图显示的最小和平均评价值,表示遗传算法的计算进度。直方图显出了基因选择的频率,即基因在当前个体中被选择的次数。参数图表示评价函数和变量值,非常方便地看到评价函数和变量值的相关关系。

首先,安装genalg包。


> install.packages("genalg")
> library(genalg)

查看一下rbga()函数的定义。


> rbga(stringMin=c(), stringMax=c(),
     suggestions=NULL,
     popSize=200, iters=100,
     mutationChance=NA,
     elitism=NA,
     monitorFunc=NULL, evalFunc=NULL,
     showSettings=FALSE, verbose=FALSE)

参数说明:

  • stringMin,设置每个基因的最小值
  • stringMax,设置每个基因的最大值
  • suggestions,建议染色体的可选列表
  • popSize,个体数量,即染色体数目,默认为200
  • iters,迭代次数,默认为100
  • mutationChance,突变机会,默认为1/(size+1),它影响收敛速度和搜索空间的探测,低机率导致更快收敛,高机率增加了搜索空间的跨度。
  • elitism,精英数量,默认为20%,直接复制到下一代的染色体数目
  • monitorFunc,监控函数,每产生一代后运行
  • evalFunc,适应度函数,用于给个体进行评价
  • showSettings,打印设置,默认为false
  • verbose,打印算法运行日志,默认为false

接下来,我们给定一个优化的问题,通过rbga()函数,计算最优化的解。

题目2:设fx=abs(x1-sqrt(exp(1)))+abs(x2-log(pi)),计算fx的最小值,其中x1,x2为2个不同的变量。

从直观上看,如果想得到fx的最小值,其实当x1=sqrt(exp(1))=1.648721, x2=log(pi)=1.14473时,fx=0为最小值。同样地,如果使用穷举法,通过循环的方法找到这2个变量,估计会很费时的,我也不做测试了。下面我们看一下rbga()函数的遗传算法的运行情况。


# 定义适应度函数
> f<-function(x){}  #代码省略

# 定义监控函数
> monitor <- function(obj){}  #代码省略

# 运行遗传算法
> m2 = rbga(c(1,1),
+           c(3,3),
+           popSize=100,
+           iters=1000,
+           evalFunc=f,
+           mutationChance=0.01,
+           verbose=TRUE,
+           monitorFunc=monitor
+           )

Testing the sanity of parameters...
Not showing GA settings...
Starting with random values in the given domains...
Starting iteration 1
Calucating evaluation values... ....................................................... done.
Sending current state to rgba.monitor()...
Creating next generation...
  sorting results...
  applying elitism...
  applying crossover...
  applying mutations... 2 mutations applied
Starting iteration 2
Calucating evaluation values... ...................... done.
Sending current state to rgba.monitor()...
Creating next generation...
  sorting results...
  applying elitism...
  applying crossover...
  applying mutations... 4 mutations applied
Starting iteration 3
Calucating evaluation values... ....................... done.

# 省略输出...

程序运行截图

rbga-run

需要注意的是,程序在要命令行界面运行,如果在RStudio中运行,会出现下面的错误提示。


Creating next generation...
  sorting results...
  applying elitism...
  applying crossover...
  applying mutations... 1 mutations applied
Starting iteration 10 
Calucating evaluation values... .......................... done.
Sending current state to rgba.monitor()...
Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
  object 'rversion' not found
Graphics error: Error in get(name, envir = asNamespace(pkg), inherits = FALSE) : 
  object 'rversion' not found

我们迭代1000次后,查看计算结果。


# 计算结果
> m2$population[1,]
[1] 1.650571 1.145784

我们得到的最优化的结果为x1=1.650571, x2=1.145784,非常接近最终的结果。另外,我们可以通过genalg包的可视化功能,看到迭代过程的每次的计算结果。下面截图分为对应1次迭代,10次迭代,200次迭代和1000次迭代的计算结果。从图中可以看出,随着迭代次数的增加,优选出的结果集变得越来越少,而且越来越精准。

rbga-all

默认图输出,用于描述遗传过程的进展,X轴为迭代次数,Y轴评价值,评价值越接近于0越好。在1000迭代1000次后,基本找到了精确的结果。


> plot(m2)

rbga-ev

直方图输出,用于描述对染色体的基因选择频率,即一个基因在染色体中的当前人口被选择的次数。当x1在1.65区域时,被选择超过80次;当x2在1.146区域时,被选择超过了80次。通过直方图,我们可以理解为更优秀的基因被留给了后代。


> plot(m2,type='hist')

rbga-hist

参数图输出,用于描述评价函数和变量的值的相关关系。对于x1,评价值越小,变量值越准确,但相关关系不明显。对于x2,看不出相关关系。


> plot(m2,type='vars')

rbga-vars

对比mcga包和genalg包,mcga包适合计算大范围取值空间的最优解,而用genalg包对于大范围取值空间的计算就表现就不太好了。从另一个方面讲,genalg包提供了可视化工具,可以让我们直观的看遗传算法的收敛过程,对于算法的理解和调优是非常有帮助的。

在掌握了遗传算法后,我们就可以快度地处理一些优化的问题了,比如接下来我会介绍的金融回测系统的参数寻优方案。让我们远离穷举法,珍惜CPU的每一秒时间。

参考文章
http://zh.wikipedia.org/zh/遗传算法

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

打赏作者

桶排序的Nodejs实现

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

转载请注明出处:
http://blog.fens.me/algorithm-bucketsort-nodejs/

bucket-nodejs

前言

一个好的程序算法,可以提升百倍的程序性能。但并没有一种通用的算法,适用于所有场景。选择合适的算法,用在正确的地方,是一个好算法的开始。

本文将用Nodejs实现桶排序算法。

目录

  1. 桶排序介绍
  2. 桶排序算法演示
  3. Nodejs程序实现
  4. 案例:桶排序统计高考分数

1. 桶排序介绍

桶排序(Bucket sort)是一种基于计数的排序算法,工作的原理是将数据分到有限数量的桶子里,然后每个桶再分别排序(有可能再使用别的排序算法或是以递回方式继续使用桶排序进行排序)。当要被排序的数据内的数值是均匀分配的时候,桶排序时间复杂度为Θ(n)。桶排序不同于快速排序,并不是比较排序,不受到时间复杂度 O(nlogn) 下限的影响。关于快速排序,请参考文章:快速排序的Nodejs实现

桶排序按下面4步进行:

  • 1. 设置固定数量的空桶。
  • 2. 把数据放到对应的桶中。
  • 3. 对每个不为空的桶中数据进行排序。
  • 4. 拼接从不为空的桶中数据,得到结果。

桶排序,主要适用于小范围整数数据,且独立均匀分布,可以计算的数据量很大,而且符合线性期望时间。

2. 桶排序算法演示

举例来说,现在有一组数据[7, 36, 65, 56, 33, 60, 110, 42, 42, 94, 59, 22, 83, 84, 63, 77, 67, 101],怎么对其按从小到大顺序排序呢?

bucketsort

操作步骤说明:

  • 1. 设置桶的数量为5个空桶,找到最大值110,最小值7,每个桶的范围20.8=(110-7+1)/5 。
  • 2. 遍历原始数据,以链表结构,放到对应的桶中。数字7,桶索引值为0,计算公式为floor((7 – 7) / 20.8), 数字36,桶索引值为1,计算公式floor((36 – 7) / 20.8)。
  • 3. 当向同一个索引的桶,第二次插入数据时,判断桶中已存在的数字与新插入数字的大小,按照左到右,从小到大的顺序插入。如:索引为2的桶,在插入63时,桶中已存在4个数字56,59,60,65,则数字63,插入到65的左边。
  • 4. 合并非空的桶,按从左到右的顺序合并0,1,2,3,4桶。
  • 5. 得到桶排序的结构

3. Nodejs程序实现

像桶排序这种成熟的算法,自己实现一下并不难,按照上文的思路,我写了一个简单的程序实现。我感觉其中最麻烦的部分,是用Javascript操作链表。

现实代码如下:


'use strict';

/////////////////////////////////////////////////
// 桶排序
/////////////////////////////////////////////////
var _this = this
    , L = require('linklist');//链表

/**
 * 普通数组桶排序,同步
 *
 * @param arr Array 整数数组
 * @param num 桶的个数
 *
 * @example:
 * sort([1,4,1,5,3,2,3,3,2,5,2,8,9,2,1],5)
 * sort([1,4,1,5,3,2,3,3,2,5,2,8,9,2,1],5,0,5)
 */
exports.sort = function (arr, count) {
    if (arr.length == 0) return [];
    count = count || (count > 1 ? count : 10);

    // 判断最大值、最小值
    var min = arr[0], max = arr[0];
    for (var i = 1; i < arr.length; i++) {
        min = min < arr[i] ? min : arr[i];
        max = max > arr[i] ? max : arr[i];
    }
    var delta = (max - min + 1) / count;
    // console.log(min+","+max+","+delta);

    //初始化桶
    var buckets = [];

    //存储数据到桶
    for (var i = 0; i < arr.length; i++) {
        var idx = Math.floor((arr[i] - min) / delta); // 桶索引

        if (buckets[idx]) {//非空桶
            var bucket = buckets[idx];
            var insert = false;//插入标石
            L.reTraversal(bucket, function (item, done) {
                if (arr[i] <= item.v) {//小于,左边插入
                    L.append(item, _val(arr[i]));
                    insert = true;
                    done();//退出遍历
                }
            });
            if (!insert) { //大于,右边插入
                L.append(bucket, _val(arr[i]));
            }
        } else {//空桶
            var bucket = L.init();
            L.append(bucket, _val(arr[i]));
            buckets[idx] = bucket; //链表实现
        }
    }

    var result = [];
    for (var i = 0, j = 0; i < count; i++) {
        L.reTraversal(buckets[i], function (item) {
            // console.log(i+":"+item.v);
            result[j++] = item.v;
        });
    }
    return result;
}

//链表存储对象
function _val(v) {
    return {v: v}
}

运行程序:


var algo = require('./index.js');
var data = [ 7, 36, 65, 56, 33, 60, 110, 42, 42, 94,  59, 22, 83, 84, 63, 77, 67,101 ];
console.log(data);
console.log(algo.bucketsort.sort(data,5));//5个桶
console.log(algo.bucketsort.sort(data,10));//10个桶

//output
[ 7, 36, 65, 56, 33, 60, 110, 42, 42, 94, 59, 22, 83, 84, 63, 77, 67, 101 ]
[ 7, 22, 33, 36, 42, 42, 56, 59, 60, 63, 65, 67, 77, 83, 84, 94, 101, 110 ]
[ 7, 22, 33, 36, 42, 42, 56, 59, 60, 63, 65, 67, 77, 83, 84, 94, 101, 110 ]

需要说明的是:

  • 1. 桶内排序,可以像程序中所描述的,在插入过程中实现;也可以插入不排序,在合并过程中,再进行排序,可以调用快度排序。
  • 2. 链表,在Node的底层API中,有一个链表的实现< href="https://github.com/joyent/node/blob/master/lib/_linklist.js" target="_blank" ref="nofollow">_linklist.js,我没有直接使用,而是通过linklist包调用的。

程序源代码已上传到github, https://github.com/bsspirit/ape-algorithm。下载后,可以参考example.js文件中,来调用桶排序程序。

同时也在npm发布了,安装命令:

npm install ape-algorithm

4. 案例:桶排序统计高考分数

桶排序最出名的一个应用场景,就是统计高考的分数。一年的全国高考考生人数为900万人,分数使用标准分,最低200 ,最高900 ,没有小数,如果把这900万数字进行排序,应该如何做呢?

算法分析:

  • 1. 如果使用基于比较的排序,快速排序,平均时间复杂度为O(nlogn) = O(9000000*log9000000)=144114616=1.44亿次比较。
  • 2. 如果使用基于计数的排序,桶排序,平均的时候复杂度,可以控制在线性复杂度,当创建700桶时从200分到900分各一个桶,O(N)=O(9000000),就相当于扫描一次900W条数据。

我们跑一个程序,对比一次快速排序和桶排序。


//产生100W条,[200,900]闭区间的数据
var data = algo.data.randomData(1000*1000,200,900);
var s1 = new Date().getTime();
algo.quicksort.sort(data);//快速排序
var s2 = new Date().getTime();
algo.bucketsort.sort(data,700);//装到700个桶
var s3 = new Date().getTime();

console.log("quicksort time: %sms",s2-s1);
console.log("bucket time: %sms",s3-s2);

//output
quicksort time: 14768ms
bucket time: 1089ms

所以,对于高考计分的案例来说,桶排序是更适合的!我们把合适的算法,用在适合的场景,会给程序带来超越硬件的性能提升。

转载请注明出处:
http://blog.fens.me/algorithm-bucketsort-nodejs/

打赏作者

快速排序的Nodejs实现

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

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

转载请注明出处:
http://blog.fens.me/algorithm-quicksort-nodejs/

quicksort-nodejs

前言

排序算法,是程序开发中最基本的一项任务。教科课上讲的排序算法大概有7-8种,快速排序(QuickSort)是使用最广泛的一种基于比较排序的算法。网上已经有了各种语言的实现。本文则是利用快速排序的原理,实现对数组对象中的属性进行排序,利用Nodejs来实现。

目录

  1. 快速排序介绍
  2. 快速排序算法演示
  3. Nodejs程序实现

1. 快速排序算法介绍

快速排序是由东尼·霍尔所发展的一种排序算法。在平均状况下,排序 n 个项目要Ο(n log n)次比较。在最坏状况下则需要Ο(n^2)次比较,但这种状况并不常见。事实上,快速排序通常明显比其他Ο(n log n) 算法更快,因为它的内部循环(inner loop)可以在大部分的架构上很有效率地被实现出来。快速排序使用分治法(Divide and conquer)策略来把一个串行(list)分为两个子串行(sub-lists)。

快速排序的思想很简单,排序过程只需要三步:

  • 1. 从数列中挑出一个元素,称为 “基准”(pivot)
  • 2. 重新排序数列,所有元素比基准值小的摆放在基准前面,所有元素比基准值大的摆在基准的后面(相同的数可以到任一边)。在这个分区退出之后,该基准就处于数列的中间位置。这个称为分区(partition)操作。
  • 3. 递归地(recursive)把小于基准值元素的子数列和大于基准值元素的子数列排序。

递归的最底部情形,是数列的大小是零或一,也就是永远都已经被排序好了。虽然一直递归下去,但是这个算法总会退出,因为在每次的迭代(iteration)中,它至少会把一个元素摆到它最后的位置去。

摘自: http://zh.wikipedia.org/wiki/快速排序

2. 快速排序算法演示

举例来说,现在有一组数据[9,23,12,11,43,54,43,2,12,66],怎么对其按从小到大顺序排序呢?

quicksort

  • Step1,选择第一个元素9作为分隔值,把小于9的数字放左边,大于等于9的数字放右边。原来的一个数组就被分成了3个部分,左边+分隔值+右边=[2]+9+[23 12 11 43 54 43 12 66]
  • Step2, 选择Step1数组左边的第一元素2作为分隔值,左边只有一个元素,左边排序完成。选择Step1数组右边的第一个元素23作为分隔值,把小于23的数字放左边,大于等于23的数字放右边,得到[ 12 11 12 ] 23 [ 43 54 43 66 ]。
  • Step3,选择Step2数组以23分隔的左边第一个元素12作为分隔值,把小于12的数字放左边,大于等于12的数字放右边,得到[ 11 ] 12 [ 12 ]。选择Step2数组以23分隔的左边第一个元素43作为分割值,把小于43的数字放左边,大于等于43的数字右边,得到 [] 43 [ 54 43 66 ]。
  • Step4,对Step3数组以12分隔的左右两边进行排序,都只有一个元素,排序完成。对Step3数组以43分隔左右两进行排序,左边无元素,右边第一个元素54作为分隔值,得到[ 43 ] 54 [ 66 ]。
  • Step5,对Step4数组以54分隔的左右两边进行排序,都只有一个元素,排序完成。
  • 最后,后整个数组拼接在一起就得到排序结果:[2 9 11 12 12 23 43 43 54 66 ]

我们看到快速排序,其实就是一种分而治之的办法。网上还有很多快速排序的优化方法,可以尽一步减少时间复杂度和空间复杂度的开销。

3. Nodejs程序实现

像快速排序这种成熟的算法,我以为在Nodejs中早就有人已经封装好了程序包。等我用到的时候竟然找不到,所以自己花了点时间,动手实现了快速排序的Node模块。

当然,除了支持对普通数组的快速排序,还支持按数组对象中的某个属性进行快速排序,已满足自己的功能需求。

3.1 普通数组的快速排序

这部分的实现就和上文中演示的思路是一样的,我另外加了方向的判断。可以从小到大排序,也可以从大到小排序。

现实代码如下:


/**
 * 普通数组快速排序
 *
 * @param arr Array 数字数组
 * @param dir asc升序、desc降序
 *
 * @example:
 * sort([1,4,2,5])
 * sort([1,4,2,5],'asc')
 * sort([1,4,2,5],'desc')
 */
exports.sort=function(arr,dir){
    dir=dir||'asc';
    if (arr.length == 0) return [];

    var left = new Array();
    var right = new Array();
    var pivot = arr[0];

    if(dir==='asc'){//升序
        for (var i = 1; i < arr.length; i++) {
            arr[i] < pivot ? left.push(arr[i]): right.push(arr[i]);
        }
    }else{//降序
        for (var i = 1; i < arr.length; i++) {
            arr[i] > pivot ? left.push(arr[i]): right.push(arr[i]);
        }
    }
    return _this.sort(left,dir).concat(pivot, _this.sort(right,dir));
}

运行程序:


var algo = require('./index.js');
var arr = [23,12,11,43,54,43,2,12,66];
console.log(arr);
console.log(algo.quicksort.sort(arr));
console.log(algo.quicksort.sort(arr,'desc'));

//output
[ 23, 12, 11, 43, 54, 43, 2, 12, 66 ]
[ 2, 11, 12, 12, 23, 43, 43, 54, 66 ]
[ 66, 54, 43, 43, 23, 12, 12, 11, 2 ]

3.2 按数组对象中的属性进行快速排序

按数组对象中的属性快速排序,是指在数组中存储的不是数字类型,而是对象类型,如[{name:’小B’,id:12},{name:’小C’,id:21},{name:’小A’,id:2}],然后想要按对象中的id属性把对象排序。这样需求,其实在程序开发中更常见。那我们怎么做呢?其实,原理是一样的,只要把分隔值(pivot)和存储值(pivotObj)分成2个变量就行了。

实现代码如下:


/**
 * 对象数组快速排序
 *
 * @param arr Array 对象数组
 * @param key string 用于排序的属性
 * @param dir asc升序、desc降序
 *
 * @example:
 * sort([{name:'b',id:12},{name:'c',id:21},{name:'a',id:2}],'id')
 * sort([{name:'b',id:12},{name:'c',id:21},{name:'a',id:2}],'id','asc')
 * sort([{name:'b',id:12},{name:'c',id:21},{name:'a',id:2}],'id','desc')
 */
exports.sortObj=function(arr,key,dir){
    key=key||'id';
    dir=dir||'asc';
    if (arr.length == 0) return [];

    var left = new Array();
    var right = new Array();
    var pivot = arr[0][key];//分割值
    var pivotObj = arr[0];//存储值

    if(dir==='asc'){//升序
        for (var i = 1; i < arr.length; i++) {
            arr[i][key] < pivot ? left.push(arr[i]): right.push(arr[i]);
        }
    }else{//降序
        for (var i = 1; i < arr.length; i++) {
            arr[i][key] > pivot ? left.push(arr[i]): right.push(arr[i]);
        }
    }
    return _this.sortObj(left,key,dir).concat(pivotObj, _this.sortObj(right,key,dir));
}

运行程序:


var algo = require('./index.js');
var arrObj = [{name:'b',id:12},{name:'c',id:21},{name:'a',id:2}];
console.log(arrObj);
console.log(algo.quicksort.sortObj(arrObj,'id','asc'));
console.log(algo.quicksort.sortObj(arrObj,'id','desc'));

//output
[ { name: 'b', id: 12 },{ name: 'c', id: 21 },{ name: 'a', id: 2 } ]
[ { name: 'a', id: 2 },{ name: 'b', id: 12 },{ name: 'c', id: 21 } ]
[ { name: 'c', id: 21 },{ name: 'b', id: 12 },{ name: 'a', id: 2 } ]

不用多少代码,就能实现快速排序。程序源代码已上传到github, https://github.com/bsspirit/ape-algorithm。下载后,可以参考example.js文件中,来调用快速排序程序。

同时也在npm发布了,安装命令:

npm install ape-algorithm

有时候写写算法,感觉又回到了学生时代,还是挺有意思的。

转载请注明出处:
http://blog.fens.me/algorithm-quicksort-nodejs/

打赏作者

PeopleRank从社交网络中发现个体价值

Hadoop家族系列文章,主要介绍Hadoop家族产品,常用的项目包括Hadoop, Hive, Pig, HBase, Sqoop, Mahout, Zookeeper, Avro, Ambari, Chukwa,新增加的项目包括,YARN, Hcatalog, Oozie, Cassandra, Hama, Whirr, Flume, Bigtop, Crunch, Hue等。

从2011年开始,中国进入大数据风起云涌的时代,以Hadoop为代表的家族软件,占据了大数据处理的广阔地盘。开源界及厂商,所有数据软件,无一不向Hadoop靠拢。Hadoop也从小众的高富帅领域,变成了大数据开发的标准。在Hadoop原有技术基础之上,出现了Hadoop家族产品,通过“大数据”概念不断创新,推出科技进步。

作为IT界的开发人员,我们也要跟上节奏,抓住机遇,跟着Hadoop一起雄起!

关于作者:

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

转载请注明出处:
http://blog.fens.me/hadoop-social-peoplerank/

peoplerank

前言

如果说Google改变了互联网,那么社交网络就改变人们的生活方式。通过社交网络,我们每个个体,都是成为了网络的中心。我们的生活半径,被无限放大,通过6个朋友关系,就可以认识世界上任何一个人。

未来的互联网将是属于我们每一个人。

目录

  1. PeopleRank和PageRank
  2. 需求分析:从社交网络中发现个体价值
  3. 算法模型:PeopleRank算法
  4. 架构设计:PeopleRank计算引擎系统架构
  5. 程序开发:PeopleRank算法实现

1. PeopleRank和PageRank

PageRank让Google成为搜索领域的No.1,也是当今最有影响力的互联网公司之一,用技术创新改变人们的生活。PageRank主要用于网页评分计算,把互联网上的所有网页都进行打分,给网页价值的体现。

自2012以来,中国开始进入社交网络的时代,开心网,人人网,新浪微博,腾讯微博,微信等社交网络应用,开始进入大家的生活。最早是由“抢车位”,“偷菜”等社交游戏带动的社交网络的兴起,如今人们会更多的利用社交网络,获取信息和分享信息。我们的互联网,正在从以网页信息为核心的网络,向着以人为核心的网络转变着。

于是有人就提出了,把PageRank模型应用于社交网络,定义以人为核心的个体价值。这样PageRank模型就有了新的应用领域,同时也有了一个新的名字PeopleRank。

关于PageRank的介绍,请参考文章:PageRank算法R语言实现

注:PeopleRank网上还有不同的解释,我这里仅仅表示用来解释“PageRank模型”。

下面我们将从一个PeopleRank的案例来解释,如何从社交网络中发现个体价值。

2. 需求分析:从社交网络中发现个体价值

案例介绍:
以新浪微博为例,给微博中每个用户进行评分!
从新浪微博上,把我们的关注和粉丝的关系都找到。

如下图:我在微博上,随便找了几个微博账号。

weibo-logo

我们的任务是,需要给这些账号评分!

  • 方法一,简单求和:评分=关注数+粉丝数+微博数
  • 方法二,加权求和:评分=a*关注数+b*粉丝数+c*微博数

新建数据文件:weibo.csv


~ vi weibo.csv

A,314,1091,1488
B,196,10455,327
C,557,51635228,13793
D,55,14655464,1681
E,318,547,4899
F,166,145,170
G,17,890,169
H,54,946759,17229

R语言读入数据文件


weibo<-read.csv(file="weibo.csv",header=FALSE)
names(weibo)<-c("id","follow","fans","tweet")

1). 方法一,简单相加法


> data.frame(weibo[1],rank=rowSums(weibo[2:4]))
  id     rank
1  A     2893
2  B    10978
3  C 51649578
4  D 14657200
5  E     5764
6  F      481
7  G     1076
8  H   964042

这种方法简单粗暴的方式,是否能代码公平的打分呢?!

2). 方法二,加权求和

通过a,b,c的3个参数,分别设置权重求和。与方法一存在同样的问题,a,b,c的权值都是人为指定的,也是不能表示公平的打分的。

除了上面的两个方法,你能否想到不一样的思路呢!

3. 算法模型:PeopleRank算法

基于PageRank的理论,我们以每个微博账户的“关注”为链出链接,“粉丝”为链入链接,我们把这种以人为核心的关系,叫PeopleRank。

关于PageRank的介绍,请参考文章:PageRank算法R语言实现

PeopleRank假设条件:

  • 数量假设:如果一个用户节点接收到的其他用户“关注”的数量越多,那么这个用户越重要。
  • 质量假设:用户A的“粉丝”质量不同,质量高的“粉丝”会通“关注”接向其他用户传递更多的权重。所以越是质量高的“用户”关注用户A,则用户A越重要。

衡量PeopleRank的3个指标:

  • 粉丝数
  • 粉丝是否有较高PeopleRank值
  • 粉丝关注了多少人

我们以下的数据为例,构造基于微博的数据模型:
(由于微博数据已增加访问权限,我无法拿到当前的实际数据,我用以前@晒粉丝应用,收集到的微博数据为例,这里ID已经过处理)

测试数据集:people.csv

  • 25个用户
  • 66个关系,关注和粉丝的关系

数据集: people.csv


1,19
1,21
2,11
2,17
2,21
3,1
3,20
3,2
3,7
3,6
3,10
4,3
4,6
5,19
5,11
5,2
6,4
6,12
6,18
6,15
6,10
6,5
7,9
7,18
7,10
8,3
8,11
8,7
8,16
8,14
9,6
10,8
10,18
11,13
11,3
12,9
12,4
12,16
12,5
13,19
13,1
13,6
14,7
14,17
14,19
14,1
14,5
14,2
15,11
15,14
15,12
16,20
17,4
17,6
18,10
18,11
18,15
18,14
19,18
20,10
20,5
21,24
22,11
23,17
24,15
25,24

第一列为用户ID,第二列也是用户ID。第一列用户,关注了第二用户。

以R语言可视化输出

igraph-refs

R语言程序


library(igraph)
people<-read.csv(file="people.csv",header=FALSE)
drawGraph<-function(data){
  names(data)<-c("from","to") 
  g <- graph.data.frame(data, directed=TRUE)
  V(g)$label <- V(g)$name
  V(g)$size <- 15
  E(g)$color <- grey(0.5)
  g2 <- simplify(g)
  plot(g2,layout=layout.circle)
}
drawGraph(people)

用R语言构建PeopleRank的算法原型

  • 构建邻接矩阵
  • 变换概率矩阵
  • 递归计算矩阵特征值
  • 标准化结果
  • 对结果排序输出

R语言算法模型


#构建邻接矩阵
adjacencyMatrix<-function(pages){
  n<-max(apply(pages,2,max))
  A <- matrix(0,n,n)
  for(i in 1:nrow(pages)) A[pages[i,]$dist,pages[i,]$src]<-1
  A
}

#变换概率矩阵
dProbabilityMatrix<-function(G,d=0.85){
  cs <- colSums(G)
  cs[cs==0] <- 1
  n <- nrow(G)
  delta <- (1-d)/n
  A <- matrix(delta,n,n)
  for (i in 1:n) A[i,] <- A[i,] + d*G[i,]/cs
  A
}

#递归计算矩阵特征值
eigenMatrix<-function(G,iter=100){
  n<-nrow(G)
  x <- rep(1,n)
  for (i in 1:iter) x <- G %*% x
  x/sum(x)
}

#直接计算矩阵特征值
calcEigenMatrix<-function(G){
  x <- Re(eigen(G)$vectors[,1])
  x/sum(x)
}

PeopleRank计算,带入数据集people.csv


people<-read.csv(file="people.csv",header=FALSE)
names(people)<-c("src","dist");people
A<-adjacencyMatrix(people);A
G<-dProbabilityMatrix(A);G
q<-calcEigenMatrix(G);

q
[1] 0.03274732 0.03404052 0.05983465 0.03527074 0.04366519 0.07042752 0.02741232
 [8] 0.03378595 0.02118713 0.06537870 0.07788465 0.03491910 0.03910097 0.05076803
[15] 0.06685364 0.01916392 0.02793695 0.09450614 0.05056016 0.03076591 0.02956243
[22] 0.00600000 0.00600000 0.03622806 0.00600000

我们给这25用户进行打分,从高到低进行排序。

对结果排序输出:


result<-data.frame(userid=userid,PR=q[userid])
result
   userid          PR
1      18 0.09450614
2      11 0.07788465
3       6 0.07042752
4      15 0.06685364
5      10 0.06537870
6       3 0.05983465
7      14 0.05076803
8      19 0.05056016
9       5 0.04366519
10     13 0.03910097
11     24 0.03622806
12      4 0.03527074
13     12 0.03491910
14      2 0.03404052
15      8 0.03378595
16      1 0.03274732
17     20 0.03076591
18     21 0.02956243
19     17 0.02793695
20      7 0.02741232
21      9 0.02118713
22     16 0.01916392
23     22 0.00600000
24     23 0.00600000
25     25 0.00600000

查看评分最高的用户18的关系数据:


people[c(which(people$src==18), which(people$dist==18)),]

   src dist
55  18   10
56  18   11
57  18   15
58  18   14
19   6   18
24   7   18
33  10   18
59  19   18

粉丝的PeopleRank排名:


which(result$userid %in% people$src[which(people$dist==18)])

[1]  3  5  8 20

粉丝的关注数:


table(people$src)[people$src[which(people$dist==18)]]

 6  7 10 19 
 6  3  2  1 

数据解释:用户18

  • 有4个粉丝为别是6,7,10,19。(粉丝数)
  • 4个粉丝的PeopleRank排名,是3,5,8,20。(粉丝是否有较高PeopleRank值)
  • 粉丝的关注数量,是6,3,2,1。(粉丝关注了多少人)

因此,通过对上面3个指标的综合打分,用户18是评分最高的用户。

通过R语言实现的计算模型,已经比较符合我们的评分标准了,下面我们把PeopleRank用MapReduce实现,以满足对海量数据的计算需求。

4. 架构设计:PeopleRank计算引擎系统架构

pagerank-spider

上图中,左边是数据爬虫系统,右边是Hadoop的HDFS, MapReduce。

  • 数据爬虫系统实时爬取微博数据
  • 设置系统定时器CRON,每xx小时,增量向HDFS导入数据(userid1,userid2)
  • 完成导入后,设置系统定时器,启动MapReduce程序,运行推荐算法。
  • 完成计算后,设置系统定时器,从HDFS导出推荐结果数据到数据库,方便以后的及时查询。

5. 程序开发:PeopleRank算法实现

win7的开发环境 和 Hadoop的运行环境 ,请参考文章:用Maven构建Hadoop项目

开发步骤:

  • 微博好友的关系数据: people.csv
  • 出始的PR数据:peoplerank.csv
  • 邻接矩阵: AdjacencyMatrix.java
  • PeopleRank计算: PageRank.java
  • PR标准化: Normal.java
  • 启动程序: PageRankJob.java

1). 微博好友的关系数据: people.csv,在上文中已列出

2). 出始的PR数据:peoplerank.csv


1,1
2,1
3,1
4,1
5,1
6,1
7,1
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
25,1

3). 邻接矩阵

矩阵解释:

  • 阻尼系数为0.85
  • 用户数为25
  • reduce以行输出矩阵的列,输出列主要用于分步式存储,下一步需要转成行

部分数据输出:


~ hadoop fs -cat /user/hdfs/pagerank/tmp1/part-r-00000|head -n 4

1	0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.43100002,0.005999999,0.43100002,0.005999999,0.005999999,0.005999999,0.005999999
10	0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.43100002,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.43100002,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999
11	0.005999999,0.005999999,0.43100002,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.43100002,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999
12	0.005999999,0.005999999,0.005999999,0.2185,0.2185,0.005999999,0.005999999,0.005999999,0.2185,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.2185,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999,0.005999999

4). PeopleRank计算: PageRank.java

迭代一次的PeopleRank值


~ hadoop fs -cat /user/hdfs/pagerank/pr/part-r-00000

1       0.716666
10      1.354167
11      2.232500
12      0.575000
13      0.575000
14      0.815833
15      1.354167
16      0.532500
17      1.425000
18      1.850000
19      1.283334
2       0.716667
20      1.141667
21      0.858333
22      0.150000
23      0.150000
24      1.850000
25      0.150000
3       1.170001
4       0.929167
5       1.070833
6       2.275001
7       0.603333
8       0.575000
9       0.645833

5). PR标准化: Normal.java

迭代10次,并标准化的结果:


~ hadoop fs -cat /user/hdfs/pagerank/result/part-r-00000

1       0.032842
10      0.065405
11      0.077670
12      0.034864
13      0.039175
14      0.050574
15      0.066614
16      0.019167
17      0.027990
18      0.094460
19      0.050673
2       0.034054
20      0.030835
21      0.029657
22      0.006000
23      0.006000
24      0.036111
25      0.006000
3       0.059864
4       0.035314
5       0.043805
6       0.070516
7       0.027444
8       0.033715
9       0.021251

我们对结果进行排序


   id       pr
10 18 0.094460
3  11 0.077670
22  6 0.070516
7  15 0.066614
2  10 0.065405
19  3 0.059864
11 19 0.050673
6  14 0.050574
21  5 0.043805
5  13 0.039175
17 24 0.036111
20  4 0.035314
4  12 0.034864
12  2 0.034054
24  8 0.033715
1   1 0.032842
13 20 0.030835
14 21 0.029657
9  17 0.027990
23  7 0.027444
25  9 0.021251
8  16 0.019167
15 22 0.006000
16 23 0.006000
18 25 0.006000

第一名是用户18,第二名是用户11,第三名是用户6,第三名与之前R语言单机计算的结果有些不一样,而且PR值也稍有不同,这是因为我们迭代10次时,特征值还没有完全的收敛,需要更多次的迭代计算,才能得矩阵的特征值。

程序API的实现,请参考文章:PageRank算法并行实现

我们通过PageRank的模型,成功地应用到了社交网络,实现了PeopleRank的计算,通过设计数据挖掘算法,来取代不成熟的人脑思想。算法模型将更客观,更精准。

最后,大家可以利用这个案例的设计思路,认真地了解社交网络,做出属于的自己的算法。

由于时间仓促,代码可能存在bug。请有能力同学,自行发现问题,解决问题!!

######################################################
看文字不过瘾,作者视频讲解,请访问网站:http://onbook.me/video
######################################################

转载请注明出处:
http://blog.fens.me/hadoop-social-peoplerank/

打赏作者