13518219792
可以
我们提供的服务有：成都做网站、成都网站设计、成都外贸网站建设、微信公众号开发、网站优化、网站认证、申扎ssl等。为成百上千企事业单位解决了网站和推广的问题。提供周到的售前咨询和贴心的售后服务，是有科学管理、有技术的申扎网站制作公司
最近项目中涉及基于Gradient Boosting Regression 算法拟合时间序列曲线的内容，利用python机器学习包 scikit-learn 中的GradientBoostingRegressor完成
因此就学习了下Gradient Boosting算法，在这里分享下我的理解
Boosting 算法简介
Boosting算法，我理解的就是两个思想：
1）“三个臭皮匠顶个诸葛亮”，一堆弱分类器的组合就可以成为一个强分类器；
2）“知错能改，善莫大焉”，不断地在错误中学习，迭代来降低犯错概率
当然，要理解好Boosting的思想，首先还是从弱学习算法和强学习算法来引入：
1）强学习算法：存在一个多项式时间的学习算法以识别一组概念，且识别的正确率很高；
2）弱学习算法：识别一组概念的正确率仅比随机猜测略好；
Kearns Valiant证明了弱学习算法与强学习算法的等价问题，如果两者等价，只需找到一个比随机猜测略好的学习算法，就可以将其提升为强学习算法。
那么是怎么实现“知错就改”的呢？
Boosting算法，通过一系列的迭代来优化分类结果，每迭代一次引入一个弱分类器，来克服现在已经存在的弱分类器组合的shortcomings
在Adaboost算法中，这个shortcomings的表征就是权值高的样本点
而在Gradient Boosting算法中,这个shortcomings的表征就是梯度
无论是Adaboost还是Gradient Boosting，都是通过这个shortcomings来告诉学习器怎么去提升模型，也就是“Boosting”这个名字的由来吧
Adaboost算法
Adaboost是由Freund 和 Schapire在1997年提出的，在整个训练集上维护一个分布权值向量W,用赋予权重的训练集通过弱分类算法产生分类假设（基学习器）y(x),然后计算错误率,用得到的错误率去更新分布权值向量w,对错误分类的样本分配更大的权值,正确分类的样本赋予更小的权值。每次更新后用相同的弱分类算法产生新的分类假设,这些分类假设的序列构成多分类器。对这些多分类器用加权的方法进行联合,最后得到决策结果。
其结构如下图所示：
前一个学习器改变权重w，然后再经过下一个学习器，最终所有的学习器共同组成最后的学习器。
如果一个样本在前一个学习器中被误分，那么它所对应的权重会被加重，相应地，被正确分类的样本的权重会降低。
这里主要涉及到两个权重的计算问题：
1）样本的权值
1 没有先验知识的情况下,初始的分布应为等概分布,样本数目为n,权值为1/n
2 每一次的迭代更新权值，提高分错样本的权重
2）弱学习器的权值
1 最后的强学习器是通过多个基学习器通过权值组合得到的。
2 通过权值体现不同基学习器的影响,正确率高的基学习器权重高。实际上是分类误差的一个函数
Gradient Boosting
和Adaboost不同，Gradient Boosting 在迭代的时候选择梯度下降的方向来保证最后的结果最好。
损失函数用来描述模型的“靠谱”程度，假设模型没有过拟合，损失函数越大，模型的错误率越高
如果我们的模型能够让损失函数持续的下降，则说明我们的模型在不停的改进，而最好的方式就是让损失函数在其梯度方向上下降。
下面这个流程图是Gradient Boosting的经典图了，数学推导并不复杂，只要理解了Boosting的思想，不难看懂
这里是直接对模型的函数进行更新，利用了参数可加性推广到函数空间。
训练F0-Fm一共m个基学习器，沿着梯度下降的方向不断更新ρm和am
GradientBoostingRegressor实现
python中的scikit-learn包提供了很方便的GradientBoostingRegressor和GBDT的函数接口，可以很方便的调用函数就可以完成模型的训练和预测
GradientBoostingRegressor函数的参数如下：
class sklearn.ensemble.GradientBoostingRegressor(loss='ls', learning_rate=0.1, n_estimators=100, subsample=1.0, min_samples_split=2, min_samples_leaf=1, min_weight_fraction_leaf=0.0, max_depth=3, init=None, random_state=None, max_features=None, alpha=0.9, verbose=0, max_leaf_nodes=None, warm_start=False, presort='auto')[source]¶
loss: 选择损失函数，默认值为ls(least squres)
learning_rate: 学习率，模型是0.1
n_estimators: 弱学习器的数目，默认值100
max_depth: 每一个学习器的最大深度，限制回归树的节点数目，默认为3
min_samples_split: 可以划分为内部节点的最小样本数，默认为2
min_samples_leaf: 叶节点所需的最小样本数，默认为1
……
可以参考
官方文档里带了一个很好的例子，以500个弱学习器，最小平方误差的梯度提升模型，做波士顿房价预测，代码和结果如下：
1 import numpy as np 2 import matplotlib.pyplot as plt 3 4 from sklearn import ensemble 5 from sklearn import datasets 6 from sklearn.utils import shuffle 7 from sklearn.metrics import mean_squared_error 8 9 ###############################################################################10 # Load data11 boston = datasets.load_boston()12 X, y = shuffle(boston.data, boston.target, random_state=13)13 X = X.astype(np.float32)14 offset = int(X.shape[0] * 0.9)15 X_train, y_train = X[:offset], y[:offset]16 X_test, y_test = X[offset:], y[offset:]17 18 ###############################################################################19 # Fit regression model20 params = {'n_estimators': 500, 'max_depth': 4, 'min_samples_split': 1,21 'learning_rate': 0.01, 'loss': 'ls'}22 clf = ensemble.GradientBoostingRegressor(**params)23 24 clf.fit(X_train, y_train)25 mse = mean_squared_error(y_test, clf.predict(X_test))26 print("MSE: %.4f" % mse)27 28 ###############################################################################29 # Plot training deviance30 31 # compute test set deviance32 test_score = np.zeros((params['n_estimators'],), dtype=np.float64)33 34 for i, y_pred in enumerate(clf.staged_predict(X_test)):35 test_score[i] = clf.loss_(y_test, y_pred)36 37 plt.figure(figsize=(12, 6))38 plt.subplot(1, 2, 1)39 plt.title('Deviance')40 plt.plot(np.arange(params['n_estimators']) + 1, clf.train_score_, 'b-',41 label='Training Set Deviance')42 plt.plot(np.arange(params['n_estimators']) + 1, test_score, 'r-',43 label='Test Set Deviance')44 plt.legend(loc='upper right')45 plt.xlabel('Boosting Iterations')46 plt.ylabel('Deviance')47 48 ###############################################################################49 # Plot feature importance50 feature_importance = clf.feature_importances_51 # make importances relative to max importance52 feature_importance = 100.0 * (feature_importance / feature_importance.max())53 sorted_idx = np.argsort(feature_importance)54 pos = np.arange(sorted_idx.shape[0]) + .555 plt.subplot(1, 2, 2)56 plt.barh(pos, feature_importance[sorted_idx], align='center')57 plt.yticks(pos, boston.feature_names[sorted_idx])58 plt.xlabel('Relative Importance')59 plt.title('Variable Importance')60 plt.show()
可以发现，如果要用Gradient Boosting 算法的话，在sklearn包里调用还是非常方便的，几行代码即可完成，大部分的工作应该是在特征提取上。
感觉目前做数据挖掘的工作，特征设计是最重要的，据说现在kaggle竞赛基本是GBDT的天下，优劣其实还是特征上，感觉做项目也是，不断的在研究数据中培养对数据的敏感度。
上一期提到的图像阈值处理，不仅可以实现获取你想要的目标区域(作为mask使用)，还可以帮你获取图像的边缘信息，那关于图像边缘，本期将从另外的角度来处理。
对边缘信息与背景差异较大的场景，你也可以使用threshold分割，不过若阈值不好选取，Laplacian梯度算子就不失为一直尝试方案，而且上网看看，关于Laplacian算子还可以用来判断图像的模糊程度，这个在相机的自动对焦当中，是否可以尝试判断下？
不过处理的效果并不理想，图像低灰阶部分边缘信息丢失严重。
对于sobel，laplacian算子我们可以使用cv2.filter2D()来实现，配置相应的核模板即可，如实现提取水平方向边缘信息：
你可以依据实际的应用需求来配置提取边缘的角度信息，这里以45度角(垂直向下逆时针旋转45度)为例：
对此，你可以采用下面的方式来解决：
常用形式
odeint(func, y0, t,args,Dfun)
一般这种形式就够用了。
下面是官方的例子，求解的是
D(D(y1))-t*y1=0
为了方便，采取D=d/dt。如果我们令初值
y1(0) = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
D(y1)(0) = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
这个微分方程的解y1=airy(t)。
令D(y1)=y0，就有这个常微分方程组。
D(y0)=t*y1
D(y1)=y0
Python求解该微分方程。
from scipy.integrate import odeint
from scipy.special import gamma, airy
y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
y0 = [y0_0, y1_0]
def func(y, t):
... return [t*y[1],y[0]]
def gradient(y,t):
... return [[0,t],[1,0]]
x = arange(0,4.0, 0.01)
t = x
ychk = airy(x)[0]
y = odeint(func, y0, t)
y2 = odeint(func, y0, t, Dfun=gradient)
print ychk[:36:6]
[ 0.355028 0.339511 0.324068 0.308763 0.293658 0.278806]
print y[:36:6,1]
[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]
print y2[:36:6,1]
[ 0.355028 0.339511 0.324067 0.308763 0.293658 0.278806]
得到的解与精确值相比，误差相当小。
=======================================================================================================
args是额外的参数。
用法请参看下面的例子。这是一个洛仑兹曲线的求解，并且用matplotlib绘出空间曲线图。（来自《python科学计算》）
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
# 给出位置矢量w，和三个参数p, r, b 计算出
# dx/dt, dy/dt, dz/dt 的值
x, y, z = w
# 直接与lorenz 的计算公式对应
return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.01) # 创建时间点
# 调用ode 对lorenz 进行求解, 用两个不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
# 绘图
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()
===========================================================================
scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
计算常微分方程（组）
使用 FORTRAN库odepack中的lsoda解常微分方程。这个函数一般求解初值问题。
参数：
func : callable(y, t0, ...) 计算y在t0 处的导数。
y0 : 数组 y的初值条件（可以是矢量）
t : 数组 为求出y，这是一个时间点的序列。初值点应该是这个序列的第一个元素。
args : 元组 func的额外参数
Dfun : callable(y, t0, ...) 函数的梯度(Jacobian)。即雅可比多项式。
col_deriv : boolean. True，Dfun定义列向导数（更快），否则Dfun会定义横排导数
full_output : boolean 可选输出，如果为True 则返回一个字典，作为第二输出。
printmessg : boolean 是否打印convergence 消息。
返回： y : array, shape (len(y0), len(t))
数组，包含y值，每一个对应于时间序列中的t。初值y0 在第一排。
infodict : 字典，只有full_output == True 时，才会返回。
字典包含额为的输出信息。
键值：
‘hu’ vector of step sizes successfully used for each time step.
‘tcur’ vector with the value of t reached for each time step. (will always be at least as large as the input times).
‘tolsf’ vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected.
‘tsw’ value of t at the time of the last method switch (given for each time step)
‘nst’ cumulative number of time steps
‘nfe’ cumulative number of function evaluations for each time step
‘nje’ cumulative number of jacobian evaluations for each time step
‘nqu’ a vector of method orders for each successful step.
‘imxer’index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise.
‘lenrw’ the length of the double work array required.
‘leniw’ the length of integer work array required.
‘mused’a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff)
其他参数，官方网站和文档都没有明确说明。相关的资料，暂时也找不到。
一、概观
scipy中的optimize子包中提供了常用的最优化算法函数实现。我们可以直接调用这些函数完成我们的优化问题。optimize中函数最典型的特点就是能够从函数名称上看出是使用了什么算法。下面optimize包中函数的概览：
1.非线性最优化
fmin -- 简单Nelder-Mead算法
fmin_powell -- 改进型Powell法
fmin_bfgs -- 拟Newton法
fmin_cg -- 非线性共轭梯度法
fmin_ncg -- 线性搜索Newton共轭梯度法
leastsq -- 最小二乘
2.有约束的多元函数问题
fmin_l_bfgs_b ---使用L-BFGS-B算法
fmin_tnc ---梯度信息
fmin_cobyla ---线性逼近
fmin_slsqp ---序列最小二乘法
nnls ---解|| Ax - b ||_2 for x=0
3.全局优化
anneal ---模拟退火算法
brute --强力法
4.标量函数
fminbound
brent
golden
bracket
5.拟合
curve_fit-- 使用非线性最小二乘法拟合
6.标量函数求根
brentq ---classic Brent (1973)
brenth ---A variation on the classic Brent（1980）ridder ---Ridder是提出这个算法的人名
bisect ---二分法
newton ---牛顿法
fixed_point
7.多维函数求根
fsolve ---通用
broyden1 ---Broyden’s first Jacobian approximation.
broyden2 ---Broyden’s second Jacobian approximationnewton_krylov ---Krylov approximation for inverse Jacobiananderson ---extended Anderson mixing
excitingmixing ---tuned diagonal Jacobian approximationlinearmixing ---scalar Jacobian approximationdiagbroyden ---diagonal Broyden Jacobian approximation8.实用函数
line_search ---找到满足强Wolfe的alpha值
check_grad ---通过和前向有限差分逼近比较检查梯度函数的正确性二、实战非线性最优化
fmin完整的调用形式是：
fmin(func, x0, args=(), xtol=0.0001, ftol=0.0001, maxiter=None, maxfun=None, full_output=0, disp=1, retall=0, callback=None)不过我们最常使用的就是前两个参数。一个描述优化问题的函数以及初值。后面的那些参数我们也很容易理解。如果您能用到，请自己研究。下面研究一个最简单的问题，来感受这个函数的使用方法：f(x)=x**2-4*x+8，我们知道，这个函数的最小值是4，在x=2的时候取到。
from scipy.optimize import fmin #引入优化包def myfunc(x):
return x**2-4*x+8 #定义函数
x0 = [1.3] #猜一个初值
xopt = fmin(myfunc, x0) #求解
print xopt #打印结果
运行之后，给出的结果是：
Optimization terminated successfully.
Current function value: 4.000000
Iterations: 16
Function evaluations: 32
[ 2.00001953]
程序准确的计算得出了最小值，不过最小值点并不是严格的2，这应该是由二进制机器编码误差造成的。
除了fmin_ncg必须提供梯度信息外，其他几个函数的调用大同小异，完全类似。我们不妨做一个对比：
from scipy.optimize import fmin,fmin_powell,fmin_bfgs,fmin_cgdef myfunc(x):
return x**2-4*x+8
x0 = [1.3]
xopt1 = fmin(myfunc, x0)
print xopt1
xopt2 = fmin_powell(myfunc, x0)
print xopt2
xopt3 = fmin_bfgs(myfunc, x0)
print xopt3
xopt4 = fmin_cg(myfunc,x0)
print xopt4
给出的结果是：
Optimization terminated successfully.
Current function value: 4.000000
Iterations: 16
Function evaluations: 32
[ 2.00001953]
Optimization terminated successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 53
1.99999999997
Optimization terminated successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 12
Gradient evaluations: 4
[ 2.00000001]
Optimization terminated successfully.
Current function value: 4.000000
Iterations: 2
Function evaluations: 15
Gradient evaluations: 5
[ 2.]
我们可以根据给出的消息直观的判断算法的执行情况。每一种算法数学上的问题，请自己看书学习。个人感觉，如果不是纯研究数学的工作，没必要搞清楚那些推导以及定理云云。不过，必须了解每一种算法的优劣以及能力所及。在使用的时候，不妨多种算法都使用一下，看看效果分别如何，同时，还可以互相印证算法失效的问题。
在from scipy.optimize import fmin之后，就可以使用help(fmin)来查看fmin的帮助信息了。帮助信息中没有例子，但是给出了每一个参数的含义说明，这是调用函数时候的最有价值参考。
有源码研究癖好的，或者当你需要改进这些已经实现的算法的时候，可能需要查看optimize中的每种算法的源代码。在这里：https:/ / github. com/scipy/scipy/blob/master/scipy/optimize/optimize.py聪明的你肯定发现了，顺着这个链接往上一级、再往上一级，你会找到scipy的几乎所有源码！
