Python时间序列数据分析–以示例说明

导读

本文主要分为四个部分:

  1. 用pandas处理时序数据
  2. 怎样检查时序数据的稳定性
  3. 怎样让时序数据具有稳定性
  4. 时序数据的预测

1. 用pandas导入和处理时序数据

第一步:导入常用的库

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from matplotlib.pylab import rcParams
#rcParams设定好画布的大小
rcParams['figure.figsize'] = 15, 6

第二步:导入时序数据
数据文件可在github:
http://github.com/aarshayj/Analytics_Vidhya/tree/master/Articles/Time_Series_Analysis 中下载

data = pd.read_csv(path+"AirPassengers.csv")
print data.head()
print '\n Data types:'
print data.dtypes

运行结果如下:数据包括每个月对应的passenger的数目。
可以看到data已经是一个DataFrame,包含两列Month和#Passengers,其中Month的类型是object,而index是0,1,2…

1

第三步:处理时序数据
我们需要将Month的类型变为datetime,同时作为index。

dateparse = lambda dates: pd.datetime.strptime(dates, '%Y-%m')
#---其中parse_dates 表明选择数据中的哪个column作为date-time信息,
#---index_col 告诉pandas以哪个column作为 index
#--- date_parser 使用一个function(本文用lambda表达式代替),使一个string转换为一个datetime变量
data = pd.read_csv('AirPassengers.csv', parse_dates=['Month'], index_col='Month',date_parser=dateparse)
print data.head()
print data.index

结果如下:可以看到data的index已经变成datetime类型的Month了。
2
3

2.怎样检查时序数据的稳定性(Stationarity)

因为ARIMA模型要求数据是稳定的,所以这一步至关重要。

1. 判断数据是稳定的常基于对于时间是常量的几个统计量:

  1. 常量的均值
  2. 常量的方差
  3. 与时间独立的自协方差

用图像说明如下:

  1. 均值
    Mean_nonstationary
    X是时序数据的值,t是时间。可以看到左图,数据的均值对于时间轴来说是常量,即数据的均值不是时间的函数,所有它是稳定的;右图随着时间的推移,数据的值整体趋势是增加的,所有均值是时间的函数,数据具有趋势,所以是非稳定的。
  2. 方差
    Var_nonstationary
    可以看到左图,数据的方差对于时间是常量,即数据的值域围绕着均值上下波动的振幅是固定的,所以左图数据是稳定的。而右图,数据的振幅在不同时间点不同,所以方差对于时间不是独立的,数据是非稳定的。但是左、右图的均值是一致的。
  3. 自协方差
    Cov_nonstationary
    一个时序数据的自协方差,就是它在不同两个时刻i,j的值的协方差。可以看到左图的自协方差于时间无关;而右图,随着时间的不同,数据的波动频率明显不同,导致它i,j取值不同,就会得到不同的协方差,因此是非稳定的。虽然右图在均值和方差上都是与时间无关的,但仍是非稳定数据。

2. python判断时序数据稳定性

有两种方法:
1.Rolling statistic– 即每个时间段内的平均的数据均值和标准差情况。

  1. Dickey-Fuller Test — 这个比较复杂,大致意思就是在一定置信水平下,对于时序数据假设 Null hypothesis: 非稳定。
    if 通过检验值(statistic)< 临界值(critical value),则拒绝null hypothesis,即数据是稳定的;反之则是非稳定的。
from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries):
    
    #这里以一年为一个窗口,每一个时间t的值由它前面12个月(包括自己)的均值代替,标准差同理。
    rolmean = pd.rolling_mean(timeseries,window=12)
    rolstd = pd.rolling_std(timeseries, window=12)
    
    #plot rolling statistics:
    fig = plt.figure()
    fig.add_subplot()
    orig = plt.plot(timeseries, color = 'blue',label='Original')
    mean = plt.plot(rolmean , color = 'red',label = 'rolling mean')
    std = plt.plot(rolstd, color = 'black', label= 'Rolling standard deviation')
    
    plt.legend(loc = 'best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    
    #Dickey-Fuller test:
    
    print 'Results of Dickey-Fuller Test:'
    dftest = adfuller(timeseries,autolag = 'AIC')
    #dftest的输出前一项依次为检测值,p值,滞后数,使用的观测数,各个置信度下的临界值
    dfoutput = pd.Series(dftest[0:4],index = ['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dftest[4].items():
        dfoutput['Critical value (%s)' %key] = value
    
    print dfoutput
    
ts = data['#Passengers']
test_stationarity(ts)

结果如下:
1037438-20170509173335269-842159812

8
可以看到,数据的rolling均值/标准差具有越来越大的趋势,是不稳定的。
且DF-test可以明确的指出,在任何置信度下,数据都不是稳定的。

3. 让时序数据变成稳定的方法

让数据变得不稳定的原因主要有俩:

  1. 趋势(trend)-数据随着时间变化。比如说升高或者降低。
  2. 季节性(seasonality)-数据在特定的时间段内变动。比如说节假日,或者活动导致数据的异常。

由于原数据值域范围比较大,为了缩小值域,同时保留其他信息,常用的方法是对数化,取log。

ts_log = np.log(ts)
  1. 检测和去除趋势
    通常有三种方法:

    • 聚合 : 将时间轴缩短,以一段时间内星期/月/年的均值作为数据值。使不同时间段内的值差距缩小。
    • 平滑: 以一个滑动窗口内的均值代替原来的值,为了使值之间的差距缩小
    • 多项式过滤:用一个回归模型来拟合现有数据,使得数据更平滑。

本文主要使用平滑方法

Moving Average–移动平均

moving_avg = pd.rolling_mean(ts_log,12)
plt.plot(ts_log ,color = 'blue')
plt.plot(moving_avg, color='red')

9
可以看出moving_average要比原值平滑许多。

然后作差:

ts_log_moving_avg_diff = ts_log-moving_avg
ts_log_moving_avg_diff.dropna(inplace = True)
test_stationarity(ts_log_moving_avg_diff)

1011

可以看到,做了处理之后的数据基本上没有了随时间变化的趋势,DFtest的结果告诉我们在95%的置信度下,数据是稳定的。

上面的方法是将所有的时间平等看待,而在许多情况下,可以认为越近的时刻越重要。所以引入指数加权移动平均– Exponentially-weighted moving average.(pandas中通过ewma()函数提供了此功能。)

# halflife的值决定了衰减因子alpha:  alpha = 1 - exp(log(0.5) / halflife)
expweighted_avg = pd.ewma(ts_log,halflife=12)
ts_log_ewma_diff = ts_log - expweighted_avg
test_stationarity(ts_log_ewma_diff)

1213

可以看到相比普通的Moving Average,新的数据平均标准差更小了。而且DFtest可以得到结论:数据在99%的置信度上是稳定的。

  1. 检测和去除季节性
    有两种方法:

    • 1 差分化: 以特定滞后数目的时刻的值的作差
    • 2 分解: 对趋势和季节性分别建模在移除它们

Differencing–差分

ts_log_diff = ts_log - ts_log.shift()
ts_log_diff.dropna(inplace=True)
test_stationarity(ts_log_diff)

14
如图,可以看出相比MA方法,Differencing方法处理后的数据的均值和方差的在时间轴上的振幅明显缩小了。DFtest的结论是在90%的置信度下,数据是稳定的。

3.Decomposing-分解

#分解(decomposing) 可以用来把时序数据中的趋势和周期性数据都分离出来:
from statsmodels.tsa.seasonal import seasonal_decompose
def decompose(timeseries):
    
    # 返回包含三个部分 trend(趋势部分) , seasonal(季节性部分) 和residual (残留部分)
    decomposition = seasonal_decompose(timeseries)
    
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    
    plt.subplot(411)
    plt.plot(ts_log, label='Original')
    plt.legend(loc='best')
    plt.subplot(412)
    plt.plot(trend, label='Trend')
    plt.legend(loc='best')
    plt.subplot(413)
    plt.plot(seasonal,label='Seasonality')
    plt.legend(loc='best')
    plt.subplot(414)
    plt.plot(residual, label='Residuals')
    plt.legend(loc='best')
    plt.tight_layout()
    
    return trend , seasonal, residual

16
如图可以明显的看到,将original数据 拆分成了三份。Trend数据具有明显的趋势性,Seasonality数据具有明显的周期性,Residuals是剩余的部分,可以认为是去除了趋势和季节性数据之后,稳定的数据,是我们所需要的。

#消除了trend 和seasonal之后,只对residual部分作为想要的时序数据进行处理
trend , seasonal, residual = decompose(ts_log)
residual.dropna(inplace=True)
test_stationarity(residual)

17
18
如图所示,数据的均值和方差趋于常数,几乎无波动(看上去比之前的陡峭,但是要注意他的值域只有[-0.05,0.05]之间),所以直观上可以认为是稳定的数据。另外DFtest的结果显示,Statistic值原小于1%时的Critical value,所以在99%的置信度下,数据是稳定的。

4. 对时序数据进行预测

假设经过处理,已经得到了稳定时序数据。接下来,我们使用ARIMA模型
对数据已经预测。ARIMA的介绍可以见本目录下的另一篇文章。

step1: 通过ACF,PACF进行ARIMA(p,d,q)的p,q参数估计

由前文Differencing部分已知,一阶差分后数据已经稳定,所以d=1。
所以用一阶差分化的ts_log_diff = ts_log – ts_log.shift() 作为输入。
等价于

yt=YtYt1yt=Yt−Yt−1

作为输入。先画出ACF,PACF的图像,代码如下:

#ACF and PACF plots:
from statsmodels.tsa.stattools import acf, pacf
lag_acf = acf(ts_log_diff, nlags=20)
lag_pacf = pacf(ts_log_diff, nlags=20, method='ols')
#Plot ACF: 
plt.subplot(121) 
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Autocorrelation Function')

#Plot PACF:
plt.subplot(122)
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='gray')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.axhline(y=1.96/np.sqrt(len(ts_log_diff)),linestyle='--',color='gray')
plt.title('Partial Autocorrelation Function')
plt.tight_layout()

19
图中,上下两条灰线之间是置信区间,p的值就是ACF第一次穿过上置信区间时的横轴值。q的值就是PACF第一次穿过上置信区间的横轴值。所以从图中可以得到p=2,q=2。

step2: 得到参数估计值p,d,q之后,生成模型ARIMA(p,d,q)
为了突出差别,用三种参数取值的三个模型作为对比。
模型1:AR模型(ARIMA(2,1,0))

from statsmodels.tsa.arima_model import ARIMA
model = ARIMA(ts_log, order=(2, 1, 0))  
results_AR = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_AR.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_AR.fittedvalues-ts_log_diff)**2))

20
图中,蓝线是输入值,红线是模型的拟合值,RSS的累计平方误差。

模型2:MA模型(ARIMA(0,1,2))

model = ARIMA(ts_log, order=(0, 1, 2))  
results_MA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_MA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_MA.fittedvalues-ts_log_diff)**2))

21

模型3:ARIMA模型(ARIMA(2,1,2))

model = ARIMA(ts_log, order=(2, 1, 2))  
results_ARIMA = model.fit(disp=-1)  
plt.plot(ts_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')
plt.title('RSS: %.4f'% sum((results_ARIMA.fittedvalues-ts_log_diff)**2))

22
由RSS,可知模型3–ARIMA(2,1,2)的拟合度最好,所以我们确定了最终的预测模型。

step3: 将模型代入原数据进行预测
因为上面的模型的拟合值是对原数据进行稳定化之后的输入数据的拟合,所以需要对拟合值进行相应处理的逆操作,使得它回到与原数据一致的尺度。


#ARIMA拟合的其实是一阶差分ts_log_diff,predictions_ARIMA_diff[i]是第i个月与i-1个月的ts_log的差值。
#由于差分化有一阶滞后,所以第一个月的数据是空的,
predictions_ARIMA_diff = pd.Series(results_ARIMA.fittedvalues, copy=True)
print predictions_ARIMA_diff.head()
#累加现有的diff,得到每个值与第一个月的差分(同log底的情况下)。
#即predictions_ARIMA_diff_cumsum[i] 是第i个月与第1个月的ts_log的差值。
predictions_ARIMA_diff_cumsum = predictions_ARIMA_diff.cumsum()
#先ts_log_diff => ts_log=>ts_log => ts 
#先以ts_log的第一个值作为基数,复制给所有值,然后每个时刻的值累加与第一个月对应的差值(这样就解决了,第一个月diff数据为空的问题了)
#然后得到了predictions_ARIMA_log => predictions_ARIMA
predictions_ARIMA_log = pd.Series(ts_log.ix[0], index=ts_log.index)
predictions_ARIMA_log = predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)
predictions_ARIMA = np.exp(predictions_ARIMA_log)
plt.figure()
plt.plot(ts)
plt.plot(predictions_ARIMA)
plt.title('RMSE: %.4f'% np.sqrt(sum((predictions_ARIMA-ts)**2)/len(ts)))

23

5.总结

前面一篇文章,总结了ARIMA建模的步骤。
(1). 获取被观测系统时间序列数据;
(2). 对数据绘图,观测是否为平稳时间序列;对于非平稳时间序列要先进行d阶差分运算,化为平稳时间序列;
(3). 经过第二步处理,已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF,通过对自相关图和偏自相关图的分析,得到最佳的阶层 p 和阶数 q
(4). 由以上得到的d、q、p,得到ARIMA模型。然后开始对得到的模型进行模型检验。
具体例子会在另一篇文章中给出。

本文结合一个例子,说明python如何解决:
1.判断一个时序数据是否是稳定。对应步骤(1)

  1. 怎样让时序数据稳定化。对应步骤(2)
  2. 使用ARIMA模型进行时序数据预测。对应步骤(3,4)

另外对data science感兴趣的同学可以关注这个网站,干货还挺多的。
https://www.analyticsvidhya.com/blog/

GBDT(MART) 迭代决策树入门教程 | 简介

2012年11月29日 19:12:19

在网上看到一篇对从代码层面理解gbdt比较好的文章,转载记录一下:

GBDT(Gradient Boosting Decision Tree) 又叫 MART(Multiple Additive Regression Tree),是一种迭代的决策树算法,该算法由多棵决策树组成,所有树的结论累加起来做最终答案。它在被提出之初就和SVM一起被认为是泛化能力(generalization)较强的算法。近些年更因为被用于搜索排序的机器学习模型而引起大家关注。

后记:发现GBDT除了我描述的残差版本外还有另一种GBDT描述,两者大概相同,但求解方法(Gradient应用)不同。其区别和另一版本的介绍链接见这里。由于另一版本介绍博客中亦有不少错误,建议大家还是先看本篇,再跳到另一版本描述,这个顺序当能两版本都看懂。

第1~4节:GBDT算法内部究竟是如何工作的?

第5节:它可以用于解决哪些问题?

第6节:它又是怎样应用于搜索排序的呢? 

在此先给出我比较推荐的两篇英文文献,喜欢英文原版的同学可直接阅读:

【1】Boosting Decision Tree入门教程 http://www.schonlau.net/publication/05stata_boosting.pdf

【2】LambdaMART用于搜索排序入门教程 http://research.microsoft.com/pubs/132652/MSR-TR-2010-82.pdf

GBDT主要由三个概念组成:Regression Decistion Tree(即DT),Gradient Boosting(即GB),Shrinkage (算法的一个重要演进分枝,目前大部分源码都按该版本实现)。搞定这三个概念后就能明白GBDT是如何工作的,要继续理解它如何用于搜索排序则需要额外理解RankNet概念,之后便功德圆满。下文将逐个碎片介绍,最终把整张图拼出来。

一、 DT:回归树 Regression Decision Tree

提起决策树(DT, Decision Tree) 绝大部分人首先想到的就是C4.5分类决策树。但如果一开始就把GBDT中的树想成分类树,那就是一条歪路走到黑,一路各种坑,最终摔得都要咯血了还是一头雾水说的就是LZ自己啊有木有。咳嗯,所以说千万不要以为GBDT是很多棵分类树。决策树分为两大类,回归树和分类树。前者用于预测实数值,如明天的温度、用户的年龄、网页的相关程度;后者用于分类标签值,如晴天/阴天/雾/雨、用户性别、网页是否是垃圾页面。这里要强调的是,前者的结果加减是有意义的,如10岁+5岁-3岁=12岁,后者则无意义,如男+男+女=到底是男是女? GBDT的核心在于累加所有树的结果作为最终结果,就像前面对年龄的累加(-3是加负3),而分类树的结果显然是没办法累加的,所以GBDT中的树都是回归树,不是分类树,这点对理解GBDT相当重要(尽管GBDT调整后也可用于分类但不代表GBDT的树是分类树)。那么回归树是如何工作的呢?

下面我们以对人的性别判别/年龄预测为例来说明,每个instance都是一个我们已知性别/年龄的人,而feature则包括这个人上网的时长、上网的时段、网购所花的金额等。

作为对比,先说分类树,我们知道C4.5分类树在每次分枝时,是穷举每一个feature的每一个阈值,找到使得按照feature<=阈值,和feature>阈值分成的两个分枝的熵最大的feature和阈值(熵最大的概念可理解成尽可能每个分枝的男女比例都远离1:1),按照该标准分枝得到两个新节点,用同样方法继续分枝直到所有人都被分入性别唯一的叶子节点,或达到预设的终止条件,若最终叶子节点中的性别不唯一,则以多数人的性别作为该叶子节点的性别。

回归树总体流程也是类似,不过在每个节点(不一定是叶子节点)都会得一个预测值,以年龄为例,该预测值等于属于这个节点的所有人年龄的平均值。分枝时穷举每一个feature的每个阈值找最好的分割点,但衡量最好的标准不再是最大熵,而是最小化均方差–即(每个人的年龄-预测年龄)^2 的总和 / N,或者说是每个人的预测误差平方和 除以 N。这很好理解,被预测出错的人数越多,错的越离谱,均方差就越大,通过最小化均方差能够找到最靠谱的分枝依据。分枝直到每个叶子节点上人的年龄都唯一(这太难了)或者达到预设的终止条件(如叶子个数上限),若最终叶子节点上人的年龄不唯一,则以该节点上所有人的平均年龄做为该叶子节点的预测年龄。若还不明白可以Google “Regression Tree”,或阅读本文的第一篇论文中Regression Tree部分。

二、 GB:梯度迭代 Gradient Boosting

好吧,我起了一个很大的标题,但事实上我并不想多讲Gradient Boosting的原理,因为不明白原理并无碍于理解GBDT中的Gradient Boosting。喜欢打破砂锅问到底的同学可以阅读这篇英文wikihttp://en.wikipedia.org/wiki/Gradient_boosted_trees#Gradient_tree_boosting

Boosting,迭代,即通过迭代多棵树来共同决策。这怎么实现呢?难道是每棵树独立训练一遍,比如A这个人,第一棵树认为是10岁,第二棵树认为是0岁,第三棵树认为是20岁,我们就取平均值10岁做最终结论?–当然不是!且不说这是投票方法并不是GBDT,只要训练集不变,独立训练三次的三棵树必定完全相同,这样做完全没有意义。之前说过,GBDT是把所有树的结论累加起来做最终结论的,所以可以想到每棵树的结论并不是年龄本身,而是年龄的一个累加量。GBDT的核心就在于,每一棵树学的是之前所有树结论和的残差,这个残差就是一个加预测值后能得真实值的累加量。比如A的真实年龄是18岁,但第一棵树的预测年龄是12岁,差了6岁,即残差为6岁。那么在第二棵树里我们把A的年龄设为6岁去学习,如果第二棵树真的能把A分到6岁的叶子节点,那累加两棵树的结论就是A的真实年龄;如果第二棵树的结论是5岁,则A仍然存在1岁的残差,第三棵树里A的年龄就变成1岁,继续学。这就是Gradient Boosting在GBDT中的意义,简单吧。

三、 GBDT工作过程实例。

还是年龄预测,简单起见训练集只有4个人,A,B,C,D,他们的年龄分别是14,16,24,26。其中A、B分别是高一和高三学生;C,D分别是应届毕业生和工作两年的员工。如果是用一棵传统的回归决策树来训练,会得到如下图1所示结果:

现在我们使用GBDT来做这件事,由于数据太少,我们限定叶子节点做多有两个,即每棵树都只有一个分枝,并且限定只学两棵树。我们会得到如下图2所示结果:

在第一棵树分枝和图1一样,由于A,B年龄较为相近,C,D年龄较为相近,他们被分为两拨,每拨用平均年龄作为预测值。此时计算残差(残差的意思就是: A的预测值 + A的残差 = A的实际值),所以A的残差就是16-15=1(注意,A的预测值是指前面所有树累加的和,这里前面只有一棵树所以直接是15,如果还有树则需要都累加起来作为A的预测值)。进而得到A,B,C,D的残差分别为-1,1,-1,1。然后我们拿残差替代A,B,C,D的原值,到第二棵树去学习,如果我们的预测值和它们的残差相等,则只需把第二棵树的结论累加到第一棵树上就能得到真实年龄了。这里的数据显然是我可以做的,第二棵树只有两个值1和-1,直接分成两个节点。此时所有人的残差都是0,即每个人都得到了真实的预测值。

换句话说,现在A,B,C,D的预测值都和真实年龄一致了。Perfect!:

A: 14岁高一学生,购物较少,经常问学长问题;预测年龄A = 15 – 1 = 14

B: 16岁高三学生;购物较少,经常被学弟问问题;预测年龄B = 15 + 1 = 16

C: 24岁应届毕业生;购物较多,经常问师兄问题;预测年龄C = 25 – 1 = 24

D: 26岁工作两年员工;购物较多,经常被师弟问问题;预测年龄D = 25 + 1 = 26

那么哪里体现了Gradient呢?其实回到第一棵树结束时想一想,无论此时的cost function是什么,是均方差还是均差,只要它以误差作为衡量标准,残差向量(-1, 1, -1, 1)都是它的全局最优方向,这就是Gradient。

讲到这里我们已经把GBDT最核心的概念、运算过程讲完了!没错就是这么简单。不过讲到这里很容易发现三个问题:

1)既然图1和图2 最终效果相同,为何还需要GBDT呢?

答案是过拟合。过拟合是指为了让训练集精度更高,学到了很多”仅在训练集上成立的规律“,导致换一个数据集当前规律就不适用了。其实只要允许一棵树的叶子节点足够多,训练集总是能训练到100%准确率的(大不了最后一个叶子上只有一个instance)。在训练精度和实际精度(或测试精度)之间,后者才是我们想要真正得到的。

我们发现图1为了达到100%精度使用了3个feature(上网时长、时段、网购金额),其中分枝“上网时长>1.1h” 很显然已经过拟合了,这个数据集上A,B也许恰好A每天上网1.09h, B上网1.05小时,但用上网时间是不是>1.1小时来判断所有人的年龄很显然是有悖常识的;

相对来说图2的boosting虽然用了两棵树 ,但其实只用了2个feature就搞定了,后一个feature是问答比例,显然图2的依据更靠谱。(当然,这里是LZ故意做的数据,所以才能靠谱得如此狗血。实际中靠谱不靠谱总是相对的) Boosting的最大好处在于,每一步的残差计算其实变相地增大了分错instance的权重,而已经分对的instance则都趋向于0。这样后面的树就能越来越专注那些前面被分错的instance。就像我们做互联网,总是先解决60%用户的需求凑合着,再解决35%用户的需求,最后才关注那5%人的需求,这样就能逐渐把产品做好,因为不同类型用户需求可能完全不同,需要分别独立分析。如果反过来做,或者刚上来就一定要做到尽善尽美,往往最终会竹篮打水一场空。

2)Gradient呢?不是“G”BDT么?

到目前为止,我们的确没有用到求导的Gradient。在当前版本GBDT描述中,的确没有用到Gradient,该版本用残差作为全局最优的绝对方向,并不需要Gradient求解.

3)这不是boosting吧?Adaboost可不是这么定义的。

这是boosting,但不是Adaboost。GBDT不是Adaboost Decistion Tree。就像提到决策树大家会想起C4.5,提到boost多数人也会想到Adaboost。Adaboost是另一种boost方法,它按分类对错,分配不同的weight,计算cost function时使用这些weight,从而让“错分的样本权重越来越大,使它们更被重视”。Bootstrap也有类似思想,它在每一步迭代时不改变模型本身,也不计算残差,而是从N个instance训练集中按一定概率重新抽取N个instance出来(单个instance可以被重复sample),对着这N个新的instance再训练一轮。由于数据集变了迭代模型训练结果也不一样,而一个instance被前面分错的越厉害,它的概率就被设的越高,这样就能同样达到逐步关注被分错的instance,逐步完善的效果。Adaboost的方法被实践证明是一种很好的防止过拟合的方法,但至于为什么则至今没从理论上被证明。GBDT也可以在使用残差的同时引入Bootstrap re-sampling,GBDT多数实现版本中也增加的这个选项,但是否一定使用则有不同看法。re-sampling一个缺点是它的随机性,即同样的数据集合训练两遍结果是不一样的,也就是模型不可稳定复现,这对评估是很大挑战,比如很难说一个模型变好是因为你选用了更好的feature,还是由于这次sample的随机因素。

四、Shrinkage 

Shrinkage(缩减)的思想认为,每次走一小步逐渐逼近结果的效果,要比每次迈一大步很快逼近结果的方式更容易避免过拟合。即它不完全信任每一个棵残差树,它认为每棵树只学到了真理的一小部分,累加的时候只累加一小部分,通过多学几棵树弥补不足。用方程来看更清晰,即

没用Shrinkage时:(yi表示第i棵树上y的预测值, y(1~i)表示前i棵树y的综合预测值)

y(i+1) = 残差(y1~yi), 其中: 残差(y1~yi) =  y真实值 – y(1 ~ i)

y(1 ~ i) = SUM(y1, …, yi)

Shrinkage不改变第一个方程,只把第二个方程改为:

y(1 ~ i) = y(1 ~ i-1) + step * yi

即Shrinkage仍然以残差作为学习目标,但对于残差学习出来的结果,只累加一小部分(step*残差)逐步逼近目标,step一般都比较小,如0.01~0.001(注意该step非gradient的step),导致各个树的残差是渐变的而不是陡变的。直觉上这也很好理解,不像直接用残差一步修复误差,而是只修复一点点,其实就是把大步切成了很多小步。本质上,Shrinkage为每棵树设置了一个weight,累加时要乘以这个weight,但和Gradient并没有关系。这个weight就是step。就像Adaboost一样,Shrinkage能减少过拟合发生也是经验证明的,目前还没有看到从理论的证明。

五、 GBDT的适用范围

该版本GBDT几乎可用于所有回归问题(线性/非线性),相对logistic regression仅能用于线性回归,GBDT的适用面非常广。亦可用于二分类问题(设定阈值,大于阈值为正例,反之为负例)。

六、 搜索引擎排序应用 RankNet

搜索排序关注各个doc的顺序而不是绝对值,所以需要一个新的cost function,而RankNet基本就是在定义这个cost function,它可以兼容不同的算法(GBDT、神经网络…)。

实际的搜索排序使用的是LambdaMART算法,必须指出的是由于这里要使用排序需要的cost function,LambdaMART迭代用的并不是残差。Lambda在这里充当替代残差的计算方法,它使用了一种类似Gradient*步长模拟残差的方法。这里的MART在求解方法上和之前说的残差略有不同,其区别描述见这里

就像所有的机器学习一样,搜索排序的学习也需要训练集,这里一般是用人工标注实现,即对每一个(query,doc) pair给定一个分值(如1,2,3,4),分值越高表示越相关,越应该排到前面。然而这些绝对的分值本身意义不大,例如你很难说1分和2分文档的相关程度差异是1分和3分文档差距的一半。相关度本身就是一个很主观的评判,标注人员无法做到这种定量标注,这种标准也无法制定。但标注人员很容易做到的是”AB都不错,但文档A比文档B更相关,所以A是4分,B是3分“。RankNet就是基于此制定了一个学习误差衡量方法,即cost function。具体而言,RankNet对任意两个文档A,B,通过它们的人工标注分差,用sigmoid函数估计两者顺序和逆序的概率P1。然后同理用机器学习到的分差计算概率P2(sigmoid的好处在于它允许机器学习得到的分值是任意实数值,只要它们的分差和标准分的分差一致,P2就趋近于P1)。这时利用P1和P2求的两者的交叉熵,该交叉熵就是cost function。它越低说明机器学得的当前排序越趋近于标注排序。为了体现NDCG的作用(NDCG是搜索排序业界最常用的评判标准),RankNet还在cost function中乘以了NDCG。

好,现在我们有了cost function,而且它是和各个文档的当前分值yi相关的,那么虽然我们不知道它的全局最优方向,但可以求导求Gradient,Gradient即每个文档得分的一个下降方向组成的N维向量,N为文档个数(应该说是query-doc pair个数)。这里仅仅是把”求残差“的逻辑替换为”求梯度“,可以这样想:梯度方向为每一步最优方向,累加的步数多了,总能走到局部最优点,若该点恰好为全局最优点,那和用残差的效果是一样的。这时套到之前讲的逻辑,GDBT就已经可以上了。那么最终排序怎么产生呢?很简单,每个样本通过Shrinkage累加都会得到一个最终得分,直接按分数从大到小排序就可以了(因为机器学习产生的是实数域的预测分,极少会出现在人工标注中常见的两文档分数相等的情况,几乎不同考虑同分文档的排序方式)

另外,如果feature个数太多,每一棵回归树都要耗费大量时间,这时每个分支时可以随机抽一部分feature来遍历求最优(ELF源码实现方式)。

xgboost: 速度快效果好的 boosting 模型

本文作者:何通,SupStat Inc(总部在纽约,中国分部为北京数博思达信息科技有限公司)数据科学家,加拿大 Simon Fraser University 计算机学院研究生,研究兴趣为数据挖掘和生物信息学。

主页:https://github.com/hetong007

引言

在数据分析的过程中,我们经常需要对数据建模并做预测。在众多的选择中,randomForestgbm 和 glmnet 是三个尤其流行的 R 包,它们在 Kaggle 的各大数据挖掘竞赛中的出现频率独占鳌头,被坊间人称为 R 数据挖掘包中的三驾马车。根据我的个人经验,gbm 包比同样是使用树模型的 randomForest 包占用的内存更少,同时训练速度较快,尤其受到大家的喜爱。在 python 的机器学习库 sklearn 里也有 GradientBoostingClassifier 的存在。

Boosting 分类器属于集成学习模型,它基本思想是把成百上千个分类准确率较低的树模型组合起来,成为一个准确率很高的模型。这个模型会不断地迭代,每次迭代就生成一颗新的树。对于如何在每一步生成合理的树,大家提出了很多的方法,我们这里简要介绍由 Friedman 提出的 Gradient Boosting Machine。它在生成每一棵树的时候采用梯度下降的思想,以之前生成的所有树为基础,向着最小化给定目标函数的方向多走一步。在合理的参数设置下,我们往往要生成一定数量的树才能达到令人满意的准确率。在数据集较大较复杂的时候,我们可能需要几千次迭代运算,如果生成一个树模型需要几秒钟,那么这么多迭代的运算耗时,应该能让你专心地想静静……

Rocket

现在,我们希望能通过 xgboost 工具更好地解决这个问题。xgboost 的全称是 eXtreme Gradient Boosting。正如其名,它是 Gradient Boosting Machine 的一个 c++ 实现,作者为正在华盛顿大学研究机器学习的大牛陈天奇。他在研究中深感自己受制于现有库的计算速度和精度,因此在一年前开始着手搭建 xgboost 项目,并在去年夏天逐渐成型。xgboost 最大的特点在于,它能够自动利用 CPU 的多线程进行并行,同时在算法上加以改进提高了精度。它的处女秀是 Kaggle 的希格斯子信号识别竞赛,因为出众的效率与较高的预测准确度在比赛论坛中引起了参赛选手的广泛关注,在 1700 多支队伍的激烈竞争中占有一席之地。随着它在 Kaggle 社区知名度的提高,最近也有队伍借助 xgboost 在比赛中夺得第一

为了方便大家使用,陈天奇将 xgboost 封装成了 python 库。我有幸和他合作,制作了 xgboost 工具的 R 语言接口,并将其提交到了 CRAN 上。也有用户将其封装成了 julia 库。python 和 R 接口的功能一直在不断更新,大家可以通过下文了解大致的功能,然后选择自己最熟悉的语言进行学习。

功能介绍

一、基础功能

首先,我们从 github 上安装这个包:

devtools::install_github('dmlc/xgboost',subdir='R-package')

动手时间到!第一步,运行下面的代码载入样例数据:

require(xgboost)
data(agaricus.train, package='xgboost')
data(agaricus.test, package='xgboost')
train <- agaricus.train
test <- agaricus.test

这份数据需要我们通过一些蘑菇的若干属性判断这个品种是否有毒。数据以 1 或 0 来标记某个属性存在与否,所以样例数据为稀疏矩阵类型:

class(train$data)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

不用担心,xgboost 支持稀疏矩阵作为输入。下面就是训练模型的命令

bst <- xgboost(data = train$data, label = train$label, max.depth = 2, eta = 1,nround = 2, objective = "binary:logistic")
[0] train-error:0.046522
[1] train-error:0.022263

我们迭代了两次,可以看到函数输出了每一次迭代模型的误差信息。这里的数据是稀疏矩阵,当然也支持普通的稠密矩阵。如果数据文件太大不希望读进 R 中,我们也可以通过设置参数data = 'path_to_file'使其直接从硬盘读取数据并分析。目前支持直接从硬盘读取 libsvm 格式的文件。

做预测只需要一句话:

pred <- predict(bst, test$data)

做交叉验证的函数参数与训练函数基本一致,只需要在原有参数的基础上设置nfold

cv.res <- xgb.cv(data = train$data, label = train$label, max.depth = 2, eta = 1, nround = 2, objective = "binary:logistic", 
nfold = 5)
[0] train-error:0.046522+0.001102   test-error:0.046523+0.004410
[1] train-error:0.022264+0.000864   test-error:0.022266+0.003450
cv.res
   train.error.mean train.error.std test.error.mean test.error.std
1:         0.046522        0.001102        0.046523       0.004410
2:         0.022264        0.000864        0.022266       0.003450

交叉验证的函数会返回一个 data.table 类型的结果,方便我们监控训练集和测试集上的表现,从而确定最优的迭代步数。

二、高速准确

上面的几行代码只是一个入门,使用的样例数据没法表现出 xgboost 高效准确的能力。xgboost 通过如下的优化使得效率大幅提高:

  1. xgboost 借助 OpenMP,能自动利用单机 CPU 的多核进行并行计算。需要注意的是,Mac 上的 Clang 对 OpenMP 的支持较差,所以默认情况下只能单核运行。
  2. xgboost 自定义了一个数据矩阵类 DMatrix,会在训练开始时进行一遍预处理,从而提高之后每次迭代的效率。

在尽量保证所有参数都一致的情况下,我们使用希格斯子竞赛的数据做了对照实验。

Model and Parameter gbm xgboost(1 thread) xgboost(2 threads) xgboost(4 threads) xgboost(8 threads)
Time (in secs) 761.48 450.22 102.41 44.18 34.04

以上实验使用的 CPU 是 i7-4700MQ。python 的 sklearn 速度与 gbm 相仿。如果想要自己对这个结果进行测试,可以在比赛的官方网站下载数据,并参考这份 demo 中的代码。

除了明显的速度提升外,xgboost 在比赛中的效果也非常好。在这个竞赛初期,大家惊讶地发现 R 和 python 中的 gbm 竟然难以突破组织者预设的 benchmark。而 xgboost 横空出世,用不到一分钟的训练时间便打入当时的 top 10,引起了大家的兴趣与关注。准确度提升的主要原因在于,xgboost 的模型和传统的 GBDT 相比加入了对于模型复杂度的控制以及后期的剪枝处理,使得学习出来的模型更加不容易过拟合。更多算法上的细节可以参考这份陈天奇给出的介绍性讲义

三、进阶特征

除了速度快精度高,xgboost 还有一些很有用的进阶特性。下面的 “demo” 链接对应着相应功能的简单样例代码。

  1. 只要能够求出目标函数的梯度和 Hessian 矩阵,用户就可以自定义训练模型时的目标函数。demo
  2. 允许用户在交叉验证时自定义误差衡量方法,例如回归中使用 RMSE 还是 RMSLE,分类中使用 AUC,分类错误率或是 F1-score。甚至是在希格斯子比赛中的 “奇葩” 衡量标准 AMSdemo
  3. 交叉验证时可以返回模型在每一折作为预测集时的预测结果,方便构建 ensemble 模型。demo
  4. 允许用户先迭代 1000 次,查看此时模型的预测效果,然后继续迭代 1000 次,最后模型等价于一次性迭代 2000 次。demo
  5. 可以知道每棵树将样本分类到哪片叶子上,facebook 介绍过如何利用这个信息提高模型的表现。demo
  6. 可以计算变量重要性并画出树状图。demo
  7. 可以选择使用线性模型替代树模型,从而得到带 L1+L2 惩罚的线性回归或者 logistic 回归。demo

这些丰富的功能来源于对日常使用场景的总结,数据挖掘比赛需求以及许多用户给出的精彩建议。

四、未来计划

现在,机器学习工具在实用中会不可避免地遇到 “单机性能不够” 的问题。目前,xgboost 的多机分布式版本正在开发当中。基础设施搭建完成之日,便是新一轮 R 包开始设计与升级之时。

结语

我为 xgboost 制作 R 接口的目的就是希望引进好的工具,让大家使用 R 的时候心情更愉悦。总结下来,xgboost 的特点有三个:速度快,效果好,功能多,希望它能受到大家的喜爱,成为一驾新的马车。

xgboost 功能较多,参数设置比较繁杂,希望在上手之后有更全面了解的读者可以参考项目 wiki。欢迎大家多多交流,在项目 issue 区提出疑问与建议。我们也邀请有兴趣的读者提交代码完善功能,让 xgboost 成为更好用的工具。

另外,在使用 github 开发的过程中,我深切地感受到了协作写代码带来的变化。一群人在一起的时候,可以写出更有效率的代码,在丰富的使用场景中发现新的需求,在极端情况发现隐藏很深的 bug,甚至在主代码手拖延症较为忙碌的时候有人挺身而出拿下一片 issue。这样的氛围,能让一个语言社区的交流丰富起来,从而充满生命力地活跃下去。

转:时间序列预测之–ARIMA模型

什么是 ARIMA模型

ARIMA模型的全称叫做自回归移动平均模型,全称是(ARIMA, Autoregressive Integrated Moving Average Model)。也记作ARIMA(p,d,q),是统计模型(statistic model)中最常见的一种用来进行时间序列 预测的模型。

1. ARIMA的优缺点

优点: 模型十分简单,只需要内生变量而不需要借助其他外生变量。

缺点:
1.要求时序数据是稳定的(stationary),或者是通过差分化(differencing)后是稳定的。
2.本质上只能捕捉线性关系,而不能捕捉非线性关系。
注意,采用ARIMA模型预测时序数据,必须是稳定的,如果不稳定的数据,是无法捕捉到规律的。比如股票数据用ARIMA无法预测的原因就是股票数据是非稳定的,常常受政策和新闻的影响而波动。

2. 判断是时序数据是稳定的方法。

严谨的定义: 一个时间序列的随机变量是稳定的,当且仅当它的所有统计特征都是独立于时间的(是关于时间的常量)。
判断的方法:

  1. 稳定的数据是没有趋势(trend),没有周期性(seasonality)的; 即它的均值,在时间轴上拥有常量的振幅,并且它的方差,在时间轴上是趋于同一个稳定的值的。
  2. 可以使用Dickey-Fuller Test进行假设检验。(另起文章介绍)

3. ARIMA的参数与数学形式

ARIMA模型有三个参数:p,d,q。

  • p–代表预测模型中采用的时序数据本身的滞后数(lags) ,也叫做AR/Auto-Regressive项
  • d–代表时序数据需要进行几阶差分化,才是稳定的,也叫Integrated项。
  • q–代表预测模型中采用的预测误差的滞后数(lags),也叫做MA/Moving Average项

先解释一下差分: 假设y表示t时刻的Y的差分。

if d=0, yt=Ytif d=1, yt=YtYt1if d=2, yt=(YtYt1)(Yt1Yt2)=Yt2Yt1+Yt2if d=0, yt=Ytif d=1, yt=Yt−Yt−1if d=2, yt=(Yt−Yt−1)−(Yt−1−Yt−2)=Yt−2Yt−1+Yt−2

ARIMA的预测模型可以表示为:

Y的预测值 = 常量c and/or 一个或多个最近时间的Y的加权和 and/or 一个或多个最近时间的预测误差。

假设p,q,d已知,
ARIMA用数学形式表示为:

ytˆ=μ+ϕ1yt1+...+ϕpytp+θ1et1+...+θqetqyt^=μ+ϕ1∗yt−1+…+ϕp∗yt−p+θ1∗et−1+…+θq∗et−q

 

,ϕARθMA其中,ϕ表示AR的系数,θ表示MA的系数

4.ARIMA模型的几个特例

1.ARIMA(0,1,0) = random walk:

当d=1,p和q为0时,叫做random walk,如图所示,每一个时刻的位置,只与上一时刻的位置有关。
file-list
预测公式如下:

Yˆt=μ+Yt1Y^t=μ+Yt−1

2. ARIMA(1,0,0) = first-order autoregressive model:

p=1, d=0,q=0。说明时序数据是稳定的和自相关的。一个时刻的Y值只与上一个时刻的Y值有关。

Yˆt=μ+ϕ1Yt1.where, ϕ[1,1],Y^t=μ+ϕ1∗Yt−1.where, ϕ∈[−1,1],是一个斜率系数

3. ARIMA(1,1,0) = differenced first-order autoregressive model:

p=1,d=1,q=0. 说明时序数据在一阶差分化之后是稳定的和自回归的。即一个时刻的差分(y)只与上一个时刻的差分有关。

yˆt=μ+ϕ1yt1YˆtYt1=μ+ϕ1(Yt1Yt2)Yˆt=μ+Yt1+ϕ1(Yt1Yt2)y^t=μ+ϕ1∗yt−1结合一阶差分的定义,也可以表示为:Y^t−Yt−1=μ+ϕ1∗(Yt−1−Yt−2)或者Y^t=μ+Yt−1+ϕ1∗(Yt−1−Yt−2)

4. ARIMA(0,1,1) = simple exponential smoothing with growth.

p=0, d=1 ,q=1.说明数据在一阶差分后市稳定的和移动平均的。即一个时刻的估计值的差分与上一个时刻的预测误差有关。

yˆt=μ+α1et1q=1ytp=1ytyˆt=YˆtYˆt1, et1=Yt1Yˆt1,θ1=1α1Yˆt=μ+Yˆt1+α1(Yt1Yˆt1)=μ+Yt1θ1et1y^t=μ+α1∗et−1注意q=1的差分yt与p=1的差分yt的是不一样的其中,y^t=Y^t−Y^t−1, et−1=Yt−1−Y^t−1,设θ1=1−α1则也可以写成:Y^t=μ+Y^t−1+α1(Yt−1−Y^t−1)=μ+Yt−1−θ1∗et−1

5. ARIMA(2,1,2)

在通过上面的例子,可以很轻松的写出它的预测模型:

yˆt=μ+ϕ1yt1+ϕ2yt2θ1et1θ2et2:Yˆt=μ+ϕ1(Yt1Yt2)+ϕ2(Yt2Yt3)θ1(Yt1Yˆt1)θ2(Yt2Yˆt2)y^t=μ+ϕ1∗yt−1+ϕ2∗yt−2−θ1∗et−1−θ2∗et−2也可以写成:Y^t=μ+ϕ1∗(Yt−1−Yt−2)+ϕ2∗(Yt−2−Yt−3)−θ1∗(Yt−1−Y^t−1)−θ2∗(Yt−2−Y^t−2)

6. ARIMA(2,2,2)

yˆt=μ+ϕ1yt1+ϕ2yt2θ1et1θ2et2Yˆt=μ+ϕ1(Yt12Yt2+Yt3)+ϕ2(Yt22Yt3+Yt4)θ1(Yt1Yˆt1)θ2(Yt2Yˆt2)y^t=μ+ϕ1∗yt−1+ϕ2∗yt−2−θ1∗et−1−θ2∗et−2Y^t=μ+ϕ1∗(Yt−1−2Yt−2+Yt−3)+ϕ2∗(Yt−2−2Yt−3+Yt−4)−θ1∗(Yt−1−Y^t−1)−θ2∗(Yt−2−Y^t−2)

7. ARIMA建模基本步骤

  1. 获取被观测系统时间序列数据;
  2. 对数据绘图,观测是否为平稳时间序列;对于非平稳时间序列要先进行d阶差分运算,化为平稳时间序列;
  3. 经过第二步处理,已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF,通过对自相关图和偏自相关图的分析,得到最佳的阶层 p 和阶数 q
  4. 由以上得到的d、q、p,得到ARIMA模型。然后开始对得到的模型进行模型检验。
    具体例子会在另一篇文章中给出。

转自geek精神   https://www.cnblogs.com/bradleon/p/6827109.html

机器学习的六大条件

1. 你的问题能否被自动化
    问题复杂度高、问题的数量级比较大的问题,适合进行机器学习。
2. 问题是否被定义清楚
3. 是否能够获取数据
    一般来说,数据越多,结果越好。
4. 数据是否有显著的规律
    机器学习可以学习规律和固定的模式,但是对于无规律的数据是无法学习的。
5. 是否能够提炼出量化的特征
    机器只能识别量化的特征,所有的特征都要转化为量化的才能被计算机识别。
6. 能否评估模型的好坏
    模型的结果需要满足业务目标