xgboost入门与实战(原理篇)

来源:互联网 发布:中央数据交换平台 编辑:程序博客网 时间:2024/06/08 13:15

1.序

xgboost是大规模并行boosted tree的工具,它是目前最快最好的开源boosted tree工具包,比常见的工具包快10倍以上。在数据科学方面,有大量kaggle选手选用它进行数据挖掘比赛,其中包括两个以上kaggle比赛的夺冠方案。在工业界规模方面,xgboost的分布式版本有广泛的可移植性,支持在YARN, MPI, Sungrid Engine等各个平台上面运行,并且保留了单机并行版本的各种优化,使得它可以很好地解决于工业界规模的问题。

花了几天时间粗略地看完了xgboost原论文和作者的slide讲解,仅仅是入门入门入门笔记。给我的感觉就是xgboost算法比较复杂,针对传统GBDT算法做了很多细节改进,包括损失函数、正则化、切分点查找算法优化、稀疏感知算法、并行化算法设计等等。本文主要介绍xgboost基本原理以及与传统gbdt算法对比总结,后续会基于Python版本做了一些实战调参试验。想详细学习xgboost算法原理建议通读作者原始论文与slide讲解。

相关文献资料: 
Xgboost Slides 
XGBoost中文版原理介绍 
原始论文XGBoost: A Scalable Tree Boosting System 
XGBoost Parameters (official guide)

精彩博文: 
XGBoost浅入浅出——wepon 
xgboost: 速度快效果好的boosting模型 
Complete Guide to Parameter Tuning in XGBoost (with codes in Python)

XGBoost Plotting API以及GBDT组合特征实践

补充!LightGBM!:

微软出了个LightGBM,号称性能更强劲,速度更快。简单实践了一波,发现收敛速度要快一些,不过调参还不6 ,没有权威。看了GitHub上的介绍以及知乎上的一些回答,大致理解了性能提升的原因。 
主要是两个:①histogram算法替换了传统的Pre-Sorted,某种意义上是牺牲了精度(但是作者声明实验发现精度影响不大)换取速度,直方图作差构建叶子直方图挺有创造力的。(xgboost的分布式实现也是基于直方图的,利于并行)②带有深度限制的按叶子生长 (leaf-wise) 算法代替了传统的(level-wise) 决策树生长策略,提升精度,同时避免过拟合危险。

细节大家直接看作者的解释以及GitHub上的介绍吧,还是挺好理解的~ 
链接: 
https://www.zhihu.com/question/51644470/answer/130946285 
https://github.com/Microsoft/LightGBM/wiki/Features

  关于xgboost的原理网络上的资源很少,大多数还停留在应用层面,本文通过学习陈天奇博士的PPT、论文、一些网络资源,希望对xgboost原理进行深入理解。(笔者在最后的参考文献中会给出地址

2.xgboost vs gbdt

  说到xgboost,不得不说gbdt,两者都是boosting方法(如图1所示),了解gbdt可以看我这篇文章 地址。

 
图1

  如果不考虑工程实现、解决问题上的一些差异,xgboost与gbdt比较大的不同就是目标函数的定义。

  注:红色箭头指向的l即为损失函数;红色方框为正则项,包括L1、L2;红色圆圈为常数项。xgboost利用泰勒展开三项,做一个近似,我们可以很清晰地看到,最终的目标函数只依赖于每个数据点的在误差函数上的一阶导数和二阶导数。

xgboost相比传统gbdt有何不同?xgboost为什么快?xgboost如何支持并行?

看了陈天奇大神的文章和slides,略抒己见,没有面面俱到,不恰当的地方欢迎讨论:

  • 传统GBDT以CART作为基分类器,xgboost还支持线性分类器,这个时候xgboost相当于带L1和L2正则化项的逻辑斯蒂回归(分类问题)或者线性回归(回归问题)。
  • 传统GBDT在优化时只用到一阶导数信息,xgboost则对代价函数进行了二阶泰勒展开,同时用到了一阶和二阶导数。顺便提一下,xgboost工具支持自定义代价函数,只要函数可一阶和二阶求导。
  • xgboost在代价函数里加入了正则项,用于控制模型的复杂度。正则项里包含了树的叶子节点个数、每个叶子节点上输出的score的L2模的平方和。从Bias-variance tradeoff角度来讲,正则项降低了模型的variance,使学习出来的模型更加简单,防止过拟合,这也是xgboost优于传统GBDT的一个特性。
  • Shrinkage(缩减),相当于学习速率(xgboost中的eta)。xgboost在进行完一次迭代后,会将叶子节点的权重乘上该系数,主要是为了削弱每棵树的影响,让后面有更大的学习空间。实际应用中,一般把eta设置得小一点,然后迭代次数设置得大一点。(补充:传统GBDT的实现也有学习速率)
  • 列抽样(column subsampling)。xgboost借鉴了随机森林的做法,支持列抽样,不仅能降低过拟合,还能减少计算,这也是xgboost异于传统gbdt的一个特性。
  • 对缺失值的处理。对于特征的值有缺失的样本,xgboost可以自动学习出它的分裂方向。
  • xgboost工具支持并行。boosting不是一种串行的结构吗?怎么并行的?注意xgboost的并行不是tree粒度的并行,xgboost也是一次迭代完才能进行下一次迭代的(第t次迭代的代价函数里包含了前面t-1次迭代的预测值)。xgboost的并行是在特征粒度上的。我们知道,决策树的学习最耗时的一个步骤就是对特征的值进行排序(因为要确定最佳分割点),xgboost在训练之前,预先对数据进行了排序,然后保存为block结构,后面的迭代中重复地使用这个结构,大大减小计算量。这个block结构也使得并行成为了可能,在进行节点的分裂时,需要计算每个特征的增益,最终选增益最大的那个特征去做分裂,那么各个特征的增益计算就可以开多线程进行。
  • 可并行的近似直方图算法。树节点在进行分裂时,我们需要计算每个特征的每个分割点对应的增益,即用贪心法枚举所有可能的分割点。当数据无法一次载入内存或者在分布式情况下,贪心算法效率就会变得很低,所以xgboost还提出了一种可并行的近似直方图算法,用于高效地生成候选的分割点。

=============

回复 @肖岩在评论里的问题,因为有些公式放正文比较好。评论里讨论的问题的大意是 “xgboost代价函数里加入正则项,是否优于cart的剪枝”。其实陈天奇大神的slides里面也是有提到的,我当一下搬运工。
决策树的学习过程就是为了找出最优的决策树,然而从函数空间里所有的决策树中找出最优的决策树是NP-C问题,所以常采用启发式(Heuristic)的方法,如CART里面的优化GINI指数、剪枝、控制树的深度。这些启发式方法的背后往往隐含了一个目标函数,这也是大部分人经常忽视掉的。xgboost的目标函数如下:

其中正则项控制着模型的复杂度,包括了叶子节点数目T和leaf score的L2模的平方:

那这个跟剪枝有什么关系呢???
跳过一系列推导,我们直接来看xgboost中树节点分裂时所采用的公式:

这个公式形式上跟ID3算法(采用entropy计算增益) 、CART算法(采用gini指数计算增益) 是一致的,都是用分裂后的某种值 减去 分裂前的某种值,从而得到增益。为了限制树的生长,我们可以加入阈值,当增益大于阈值时才让节点分裂,上式中的gamma即阈值,它是正则项里叶子节点数T的系数,所以xgboost在优化目标函数的同时相当于做了预剪枝。另外,上式中还有一个系数lambda,是正则项里leaf score的L2模平方的系数,对leaf score做了平滑,也起到了防止过拟合的作用,这个是传统GBDT里不具备的特性。


3.原理

对于上面给出的目标函数,我们可以进一步化简

(1)定义树的复杂度

对于f的定义做一下细化,把树拆分成结构部分q叶子权重部分w。下图是一个具体的例子。结构函数q把输入映射到叶子的索引号上面去,而w给定了每个索引号对应的叶子分数是什么。

定义这个复杂度包含了一棵树里面节点的个数,以及每个树叶子节点上面输出分数的L2模平方。当然这不是唯一的一种定义方式,不过这一定义方式学习出的树效果一般都比较不错。下图还给出了复杂度计算的一个例子。

注:方框部分在最终的模型公式中控制这部分的比重,对应模型参数中的lambda ,gamma

在这种新的定义下,我们可以把目标函数进行如下改写,其中I被定义为每个叶子上面样本集合 ,g是一阶导数,h是二阶导数

这一个目标包含了T个相互独立的单变量二次函数。我们可以定义

最终公式可以化简为

通过对求导等于0,可以得到

然后把最优解代入得到:

(2)打分函数计算示例

Obj代表了当我们指定一个树的结构的时候,我们在目标上面最多减少多少。我们可以把它叫做结构分数(structure score)

(3)分裂节点

论文中给出了两种分裂节点的方法

(1)贪心法:

每一次尝试去对已有的叶子加入一个分割

对于每次扩展,我们还是要枚举所有可能的分割方案,如何高效地枚举所有的分割呢?我假设我们要枚举所有x < a 这样的条件,对于某个特定的分割a我们要计算a左边和右边的导数和。

我们可以发现对于所有的a,我们只要做一遍从左到右的扫描就可以枚举出所有分割的梯度和GL和GR。然后用上面的公式计算每个分割方案的分数就可以了。

观察这个目标函数,大家会发现第二个值得注意的事情就是引入分割不一定会使得情况变好,因为我们有一个引入新叶子的惩罚项。优化这个目标对应了树的剪枝, 当引入的分割带来的增益小于一个阀值的时候,我们可以剪掉这个分割。大家可以发现,当我们正式地推导目标的时候,像计算分数和剪枝这样的策略都会自然地出现,而不再是一种因为heuristic(启发式)而进行的操作了。

下面是论文中的算法

(2)近似算法:

主要针对数据太大,不能直接进行计算

4.自定义损失函数(指定grad、hess)

(1)损失函数

(2)grad、hess推导

(3)官方代码

#!/usr/bin/pythonimport numpy as npimport xgboost as xgb#### advanced: customized loss function#print ('start running example to used customized objective function')dtrain = xgb.DMatrix('../data/agaricus.txt.train')dtest = xgb.DMatrix('../data/agaricus.txt.test')# note: for customized objective function, we leave objective as default# note: what we are getting is margin value in prediction# you must know what you are doingparam = {'max_depth': 2, 'eta': 1, 'silent': 1}watchlist = [(dtest, 'eval'), (dtrain, 'train')]num_round = 2# user define objective function, given prediction, return gradient and second order gradient# this is log likelihood lossdef logregobj(preds, dtrain):    labels = dtrain.get_label()    preds = 1.0 / (1.0 + np.exp(-preds))    grad = preds - labels    hess = preds * (1.0-preds)    return grad, hess# user defined evaluation function, return a pair metric_name, result# NOTE: when you do customized loss function, the default prediction value is margin# this may make builtin evaluation metric not function properly# for example, we are doing logistic loss, the prediction is score before logistic transformation# the builtin evaluation error assumes input is after logistic transformation# Take this in mind when you use the customization, and maybe you need write customized evaluation functiondef evalerror(preds, dtrain):    labels = dtrain.get_label()    # return a pair metric_name, result    # since preds are margin(before logistic transformation, cutoff at 0)    return 'error', float(sum(labels != (preds > 0.0))) / len(labels)# training with customized objective, we can also do step by step training# simply look at xgboost.py's implementation of trainbst = xgb.train(param, dtrain, num_round, watchlist, logregobj, evalerror)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42

5.Xgboost调参

由于xgboost的参数过多,这里介绍三种思路

(1)GridSearch

(2)Hyperopt

(3)老外写的一篇文章,操作性比较强,推荐学习一下。地址

xgboost使用经验总结

  • 多类别分类时,类别需要从0开始编码
  • Watchlist不会影响模型训练。
  • 类别特征必须编码,因为xgboost把特征默认都当成数值型的
  • 调参:Notes on Parameter Tuning 以及 Complete Guide to Parameter Tuning in XGBoost (with codes in Python)
  • 训练的时候,为了结果可复现,记得设置随机数种子。
  • XGBoost的特征重要性是如何得到的?某个特征的重要性(feature score),等于它被选中为树节点分裂特征的次数的和,比如特征A在第一次迭代中(即第一棵树)被选中了1次去分裂树节点,在第二次迭代被选中2次…..那么最终特征A的feature score就是 1+2+….

6.工程实现优化

(1)Column Blocks and Parallelization

(2)Cache Aware Access

  • A thread pre-fetches data from non-continuous memory into a continuous bu ffer.
  • The main thread accumulates gradients statistics in the continuous buff er.

(3)System Tricks

  • Block pre-fetching.
  • Utilize multiple disks to parallelize disk operations.
  • LZ4 compression(popular recent years for outstanding performance).
  • Unrolling loops.
  • OpenMP

7.代码走读

这块非常感谢杨军老师的无私奉献【4】

个人看代码用的是SourceInsight,由于xgboost有些文件是cc后缀名,可以通过以下命令修改下(默认的识别不了)

find ./ -name "*.cc" | awk -F "." '{print $2}' | xargs -i -t mv ./{}.cc  ./{}.cpp
  • 1
  • 1

实际上,对XGBoost的源码进行走读分析之后,能够看到下面的主流程:

cli_main.cc:main()     -> CLIRunTask()          -> CLITrain()               -> DMatrix::Load()               -> learner = Learner::Create()               -> learner->Configure()               -> learner->InitModel()               -> for (i = 0; i < param.num_round; ++i)                    -> learner->UpdateOneIter()                    -> learner->Save()    learner.cc:Create()      -> new LearnerImpl()Configure()InitModel()     -> LazyInitModel()          -> obj_ = ObjFunction::Create()               -> objective.cc                    Create()                         -> SoftmaxMultiClassObj(multiclass_obj.cc)/                              LambdaRankObj(rank_obj.cc)/                              RegLossObj(regression_obj.cc)/                              PoissonRegression(regression_obj.cc)          -> gbm_ = GradientBooster::Create()               -> gbm.cc                    Create()                         -> GBTree(gbtree.cc)/                              GBLinear(gblinear.cc)          -> obj_->Configure()          -> gbm_->Configure()UpdateOneIter()      -> PredictRaw()      -> obj_->GetGradient()      -> gbm_->DoBoost()         gbtree.cc:Configure()      -> for (up in updaters)           -> up->Init()DoBoost()      -> BoostNewTrees()           -> new_tree = new RegTree()           -> for (up in updaters)                -> up->Update(new_tree)    tree_updater.cc:Create()     -> ColMaker/DistColMaker(updater_colmaker.cc)/        SketchMaker(updater_skmaker.cc)/        TreeRefresher(updater_refresh.cc)/        TreePruner(updater_prune.cc)/        HistMaker/CQHistMaker/                  GlobalProposalHistMaker/                  QuantileHistMaker(updater_histmaker.cc)/        TreeSyncher(updater_sync.cc)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56

从上面的代码主流程可以看到,在XGBoost的实现中,对算法进行了模块化的拆解,几个重要的部分分别是:

I. ObjFunction:对应于不同的Loss Function,可以完成一阶和二阶导数的计算。 
II. GradientBooster:用于管理Boost方法生成的Model,注意,这里的Booster Model既可以对应于线性Booster Model,也可以对应于Tree Booster Model。 
III. Updater:用于建树,根据具体的建树策略不同,也会有多种Updater。比如,在XGBoost里为了性能优化,既提供了单机多线程并行加速,也支持多机分布式加速。也就提供了若干种不同的并行建树的updater实现,按并行策略的不同,包括: 
  I). inter-feature exact parallelism (特征级精确并行) 
  II). inter-feature approximate parallelism(特征级近似并行,基于特征分bin计算,减少了枚举所有特征分裂点的开销) 
  III). intra-feature parallelism (特征内并行)

此外,为了避免overfit,还提供了一个用于对树进行剪枝的updater(TreePruner),以及一个用于在分布式场景下完成结点模型参数信息通信的updater(TreeSyncher),这样设计,关于建树的主要操作都可以通过Updater链的方式串接起来,比较一致干净,算是Decorator设计模式[4]的一种应用。

XGBoost的实现中,最重要的就是建树环节,而建树对应的代码中,最主要的也是Updater的实现。所以我们会以Updater的实现作为介绍的入手点。

以ColMaker(单机版的inter-feature parallelism,实现了精确建树的策略)为例,其建树操作大致如下:

updater_colmaker.cc:ColMaker::Update()     -> Builder builder;     -> builder.Update()          -> InitData()          -> InitNewNode() // 为可用于split的树结点(即叶子结点,初始情况下只有一个                           // 叶结点,也就是根结点) 计算统计量,包括gain/weight等          ->  for (depth = 0; depth < 树的最大深度; ++depth)               -> FindSplit()                    -> for (each feature) // 通过OpenMP获取                                          // inter-feature parallelism                         -> UpdateSolution()                                    -> EnumerateSplit()  // 每个执行线程处理一个特征,                                                   // 选出每个特征的                                                   // 最优split point                              -> ParallelFindSplit()                                      // 多个执行线程同时处理一个特征,选出该特征                                   //的最优split point;                                    // 在每个线程里汇总各个线程内分配到的数据样                                   //本的统计量(grad/hess);                                   // aggregate所有线程的样本统计(grad/hess),                                          //计算出每个线程分配到的样本集合的边界特征值作为                                   //split point的最优分割点;                                   // 在每个线程分配到的样本集合对应的特征值集合进                                   //行枚举作为split point,选出最优分割点                         -> SyncBestSolution()                                 // 上面的UpdateSolution()/ParallelFindSplit()                               //会为所有待扩展分割的叶结点找到特征维度的最优split                                //point,比如对于叶结点A,OpenMP线程1会找到特征F1                                //的最优split point,OpenMP线程2会找到特征F2的最                               //优split point,所以需要进行全局sync,找到叶结点A                               //的最优split point。                         -> 为需要进行分割的叶结点创建孩子结点                    -> ResetPosition()                       //根据上一步的分割动作,更新样本到树结点的映射关系                      // Missing Value(i.e. default)和非Missing Value(i.e.                       //non-default)分别处理               -> UpdateQueueExpand()                       // 将待扩展分割的叶子结点用于替换qexpand_,作为下一轮split的                      //起始基础               -> InitNewNode()  // 为可用于split的树结点计算统计量
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41

8.python、R对于xgboost的简单使用

任务:二分类,存在样本不均衡问题(scale_pos_weight可以一定程度上解读此问题)

【Python】

【R】

9.xgboost中比较重要的参数介绍

(1)objective [ default=reg:linear ] 定义学习任务及相应的学习目标,可选的目标函数如下:

  • “reg:linear” –线性回归。
  • “reg:logistic” –逻辑回归。
  • “binary:logistic” –二分类的逻辑回归问题,输出为概率。
  • “binary:logitraw” –二分类的逻辑回归问题,输出的结果为wTx。
  • “count:poisson” –计数问题的poisson回归,输出结果为poisson分布。 在poisson回归中,max_delta_step的缺省值为0.7。(used to safeguard optimization)
  • “multi:softmax” –让XGBoost采用softmax目标函数处理多分类问题,同时需要设置参数num_class(类别个数)
  • “multi:softprob” –和softmax一样,但是输出的是ndata * nclass的向量,可以将该向量reshape成ndata行nclass列的矩阵。没行数据表示样本所属于每个类别的概率。
  • “rank:pairwise” –set XGBoost to do ranking task by minimizing the pairwise loss

(2)’eval_metric’ The choices are listed below,评估指标:

  • “rmse”: root mean square error
  • “logloss”: negative log-likelihood
  • “error”: Binary classification error rate. It is calculated as #(wrong cases)/#(all cases). For the predictions, the evaluation will regard the instances with prediction value larger than 0.5 as positive instances, and the others as negative instances.
  • “merror”: Multiclass classification error rate. It is calculated as #(wrong cases)/#(all cases).
  • “mlogloss”: Multiclass logloss
  • “auc”: Area under the curve for ranking evaluation.
  • “ndcg”:Normalized Discounted Cumulative Gain
  • “map”:Mean average precision
  • “ndcg@n”,”map@n”: n can be assigned as an integer to cut off the top positions in the lists for evaluation.
  • “ndcg-“,”map-“,”ndcg@n-“,”map@n-“: In XGBoost, NDCG and MAP will evaluate the score of a list without any positive samples as 1. By adding “-” in the evaluation metric XGBoost will evaluate these score as 0 to be consistent under some conditions.

(3)lambda [default=0] L2 正则的惩罚系数

(4)alpha [default=0] L1 正则的惩罚系数

(5)lambda_bias 在偏置上的L2正则。缺省值为0(在L1上没有偏置项的正则,因为L1时偏置不重要)

(6)eta [default=0.3] 
为了防止过拟合,更新过程中用到的收缩步长。在每次提升计算之后,算法会直接获得新特征的权重。 eta通过缩减特征的权重使提升计算过程更加保守。缺省值为0.3 
取值范围为:[0,1]

(7)max_depth [default=6] 数的最大深度。缺省值为6 ,取值范围为:[1,∞]

(8)min_child_weight [default=1] 
孩子节点中最小的样本权重和。如果一个叶子节点的样本权重和小于min_child_weight则拆分过程结束。在现行回归模型中,这个参数是指建立每个模型所需要的最小样本数。该成熟越大算法越conservative 
取值范围为: [0,∞]

10.参考文献

(1)xgboost导读和实战 
(2)xgboost 
(3)自定义目标函数 
(4)机器学习算法中GBDT和XGBOOST的区别有哪些? 
(5)XGBoost: Reliable Large-scale Tree Boosting System 
(6)XGBoost: A Scalable Tree Boosting System 

1 0
原创粉丝点击