回归问题的一般模型如下:$$y=\sumw[i]*x[i]+b$$如下图所示,对于一维数据,线性回归就是根据给定的点$(x_i,y_i)$拟合出一条直线$$y=ax+b$$即求出系数a、b。
importwarningswarnings.filterwarnings('ignore')importpandasaspdimportnumpyasnpimportmatplotlib.pyplotasplt%matplotlibinlineimportseabornassnsimportmglearnmglearn.plots.plot_linear_regression_wave()w[0]:0.393906b:-0.031804
推广到多维数据,线性回归模型的训练过程就是寻找参数向量$\overrightarrow{w}$的过程,只是拟合的目标变为了高维平面,线性回归最常用的两种方法是最小二乘法(OLS)和梯度下降法。
最小二乘法基于这样一个目标,使得数据的实际值$y_i$与预测值$\hat{y_i}$之间的偏差最小,即损失函数最小,OLS使用均方误差(MeanSquareError,MSE)作为损失函数,优化目标为$$min\quadMSE=\frac1n\sum_{i=1}^n(y_i-\hat{y_i})^2$$对于一维数据而言$$MSE=\frac1n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)^2$$为求最小值,需要求偏导$$\frac{\partialMSE}{\partial\theta_0}=-\frac2n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)=0$$,$$\frac{\partialMSE}{\partial\theta_1}=-\frac2n\sum_{i=1}^n(y_i-\theta_0-\theta_1x_i)x_i=0$$联立可得$$\theta_1=\frac{\sum(y_i-\bary)(x_i-\barx)}{(x_i-\barx)^2}$$,$$\theta_0=\bary-\theta_1\barx$$
对于多元回归同理,下面是矩阵解法,损失函数定义为$$J(\mathbf\theta)=\frac{1}{2}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})$$求导$$\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta)=\mathbf{X}^T(\mathbf{X\theta}-\mathbf{Y})=0$$最终可得$$\mathbf{\theta}=(\mathbf{X^{T}X})^{-1}\mathbf{X^{T}Y}$$
最小二乘原理是通过求导的方式最小化MSE以求得参数$\theta$,下面我们介绍另一种方法梯度下降法。
梯度下降法是一个比较纯粹的计算机编程方法。
如图所示,我们知道,损失函数是系数的函数,一元线性回归有两个参数,组成了损失函数平面,我们首先随机制定一组系数,即在上图平面上随机选取一个初始点,然后同时进行以下变换$$\theta_0=\theta_0-\alpha\frac{\partialJ(\theta)}{\partial\theta_0}$$$$\theta_1=\theta_1-\alpha\frac{\partialJ(\theta)}{\partial\theta_1}$$重复该步骤直到终止。
我们来分析一下发生了什么。首先,偏导数的系数$\alpha$是正数。对于偏导数而言,当偏导大于零时候,$J(\theta)$随$\theta_i$增大而增大,同时,新的$\theta_i$小于旧的$\theta_i$,因此,$J(\theta)$减小;当偏导数小于零的时候,$J(\theta)$随着$\theta_i$增大而减小,同时,新的$\theta_i$大于旧的$\theta_i$,因此,$J(\theta)$还是减小,即每次循环,损失函数都会减小,最终到达一个局部的最小值,即极小值。
但是,我们的损失函数是凸函数,并不是有多个极小值的图形,其真实图形类似下图所示,极小值即为最小值。
算法步骤:
梯度下降法家族
python中可以使用sklearn或者statsmodel进行线性回归,sklearn偏向于机器学习,statsmodel本就是scipy.stats的补充,主要用于线性回归以及方差分析。
#使用sklearn进行线性回归fromsklearn.linear_modelimportLinearRegressionfromsklearn.model_selectionimporttrain_test_splitX,y=mglearn.datasets.make_wave(n_samples=60)#导入数据print("数据集的第一个样本点:X:{}y:{}".format(X[0],y[0]))#数据集划分,同一random_state表示对数据集进行相同的划分X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=42)#标准的sklearn风格API,Model().fit(X,y)lr=LinearRegression().fit(X_train,y_train)print('系数:{}'.format(lr.coef_))print('截距:{}'.format(lr.intercept_))print('训练精度:{}'.format(lr.score(X_train,y_train)))print('测试精度:{}'.format(lr.score(X_test,y_test)))数据集的第一个样本点:X:[-0.75275929]y:-1.1807331091906834系数:[0.39390555]截距:-0.031804343026759746训练精度:0.6700890315075756测试精度:0.65933685968637sklearn中使用OLS拟合模型,score是可决系数,可以看出,测试集可决系数只有0.65左右,可以说效果并不好,这是因为原数据为一维数据,当数据维度增加时,线性模型可以变得十分强大。下面,我们再使用statsmodel来训练模型:
对于一元线性回归,当因变量y与x并不成线性关系时,无法直接使用线性回归。根据泰勒定理:令a=0可知,y可以由$x,x^2,x^3...$线性表示,因此,可以将$x^n$看作额外的变量,将一元线性回归转化为多元线性回归,以此来增加模型的准确性。
即通过取对数将原本无线性关系的变量转化为近似线性关系以应用线性回归,如对数线性回归:$$ln\mathbf{Y}=\mathbf{X\theta}$$
岭回归使用L2正则化处理回归模型,其惩罚项为L2范数,惩罚项系数为正数,对应sklearn.Ridge中参数alpha,增大alpha会导致系数趋向于0,从而降低训练集性能,是解决过拟合的一种方法,同样的,sklearn.ridge使用OLS。$$J(\mathbf\theta)=\frac{1}{2}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\frac{1}{2}\alpha||\theta||_2^2$$
$$\frac{\partial}{\partial\mathbf\theta}J(\mathbf\theta)=\mathbf{X}^T(\mathbf{X\theta}-\mathbf{Y})+\alpha\theta=0$$
$$\theta=(X^TX+\alphaE)^{-1}X^TY$$
#岭回归在sklearn中的实现fromsklearn.linear_modelimportRidgeX,y=mglearn.datasets.load_extended_boston()print('数据规模:{}'.format(X.shape))X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=0)ridge=Ridge(alpha=1).fit(X_train,y_train)LR=LinearRegression().fit(X_train,y_train)print('线性回归精度([训练,测试]):{}'.format([LR.score(X_train,y_train),LR.score(X_test,y_test)]))print('岭回归精度([训练,测试]):{}'.format([ridge.score(X_train,y_train),ridge.score(X_test,y_test)]))数据规模:(506,104)线性回归精度([训练,测试]):[0.952051960903273,0.6074721959665708]岭回归精度([训练,测试]):[0.885796658517094,0.7527683481744755]可以看出,由于boston数据集拥有104个特征,但只有506条数据。线性回归具有十分明显的过拟合,岭回归模型训练精度低于线性回归,但是测试精度高于线性回归。
增加数据量可以缓解过拟合问题,如下图所示,当数据量增大时,线性回归的测试精度与ridge相似。
mglearn.plots.plot_ridge_n_samples()
lasso使用L1正则化,惩罚项是L1范数,但是可以使某特征系数为0,模型更容易解释,也可以呈现模型的重要特征,由于使用绝对值,存在不可导点,因此OLS、梯度下降都不可用。$$J(\mathbf\theta)=\frac{1}{2m}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\alpha||\theta||_1$$
#lasso实现fromsklearn.linear_modelimportLassolasso=Lasso(alpha=0.1,max_iter=1000000).fit(X_train,y_train)print('训练精度:{}'.format(lasso.score(X_train,y_train)))print('测试精度:{}'.format(lasso.score(X_test,y_test)))print('模型所用特征数:{}'.format(np.sum(lasso.coef_!=0)))训练精度:0.7709955157630054测试精度:0.6302009976110041模型所用特征数:8ElasticNetElasticNet同时使用L1、L2范数进行正则化,$$J(\mathbf\theta)=\frac{1}{2m}(\mathbf{X\theta}-\mathbf{Y})^T(\mathbf{X\theta}-\mathbf{Y})+\alpha\rho||\theta||_1+\frac{\alpha(1-\rho)}{2}||\theta||_2^2$$
前面提到广义线性回归:$$f(y)=X\theta$$使用线性模型进行二分类的想法就自然产生了:只需找一个单调可微函数将分类任务的真实标记y与线性回归模型的预测值联系起来。
如用于二分类的线性模型可以用以下公式预测:$$y=\sumw[i]*x[i]+b>0$$
常用的二分类模型有Logistic回归(LogisticRegression)和线性支持向量机(LinearSupportVectorMachine,线性SVM)。
LogisticRegression是将线性回归所得因变量y进行非线性(Sigmoid函数)变换映射到[0,1]之内,作为分类样本点分属0,1两类的概率。
Sigmoid函数$$g(z)=y=\frac{1}{1+e^{-z}}$$其图像如下:
X=np.linspace(-10,10)y=[]foriinX:y.append(1/(1+np.exp(-i)))plt.plot(X,y)plt.xlabel('z')plt.ylabel('g(z)')在x=0处函数值为0.5,x趋向于无穷时,函数值分别趋向0和1。如果$${z=X\theta}$$那么就把线性回归所得的函数值映射到了0-1之间,$g(z)$可以看作分类为1的概率,越靠近1,被分类为1的概率越大,在临界值0.5处最容易被误分类。
对y进行变形可得$$z=X\theta=ln\frac{y}{1-y}$$
我们将y看作分类为正例的概率,则z就可以看作正反概率之比(称为几率)的对数,因此,模型称为“对数几率回归”。
对于每个样本点$(x_i,y_i)$,$y_i=1,y_i=0$的概率分别为$$P(y_i=1|x_i,\theta)=h_\theta(x_i)$$$$P(y_i=0|x_i,\theta)=1-h_\theta(x_i)$$将其合并为$$P(y_i|x_i,\theta)=h_\theta(x_i)^{y_i}(1-h_\theta(x_i))^{1-y_i}$$假设每个样本点独立同分布,样本数为n,由最大似然法(MLE)构造似然函数得$$L(\theta)=\prod_{i=1}^nP(y_i|x_i,\theta)$$由于似然函数表示的是取得现有样本的概率,应当予以最大化,因此,取似然函数的对数的相反数作为损失函数$$J(\theta)=-lnL(\theta)=-\sum\limits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+(1-y_i)ln(1-h_{\theta}(x_i)))$$
求偏导得$$\frac{\partial}{\partial\theta}J(\theta)=X^T(h_{\theta}(X)-Y)$$
使用梯度下降法$$\theta=\theta-\alphaX^T(h_{\theta}(X)-Y)$$
#乳腺癌数据上使用LogisticRegressionfromsklearn.datasetsimportload_breast_cancerfromsklearn.linear_modelimportLogisticRegressioncancer=load_breast_cancer()X_train,X_test,y_train,y_test=train_test_split(cancer.data,cancer.target,stratify=cancer.target,random_state=42)forC,makerinzip([0.001,1,100],['o','^','v']):logistic=LogisticRegression(C=C,penalty='l2').fit(X_train,y_train)print('训练精度(C={}):{}'.format(C,logistic.score(X_train,y_train)))print('训练精度(C={}):{}'.format(C,logistic.score(X_test,y_test)))plt.plot(logistic.coef_.T,maker,label='C={}'.format(C))plt.xticks(range(cancer.data.shape[1]),cancer.feature_names,rotation=90)plt.xlabel('CoefficientIndex')plt.ylabel('Coefficient')plt.legend()训练精度(C=0.001):0.9507042253521126训练精度(C=0.001):0.9440559440559441训练精度(C=1):0.9413145539906104训练精度(C=1):0.965034965034965训练精度(C=100):0.9507042253521126训练精度(C=100):0.958041958041958
LogisticRegression也可以使用正则化,方法同样是在损失函数后增加正则化项。sklearn中LogisticRegression默认使用L2正则化,参数penalty可修改正则化方式。上图是不同正则化参数训练所得模型系数,可以看出skleran中正则化项C越小,正则化程度越强,参数的变换范围越小。这是由于取损失函数时取了相反数,-C相当于lasso中的$\alpha$,C越小,-C越大。$$J(\mathbf\theta)=-\sum\limits_{i=1}^{m}(y_iln(h_{\theta}(x_i))+(1-y_i)ln(1-h_{\theta}(x_i)))-\frac{1}{2}C||\theta||_2^2$$
在使用逻辑回归进行预测时,通常要进行特征筛选,以便减少特征数量,增强模型可解释性并加快训练速度。
#导入数据data=pd.read_excel(r'G:\MYFiles\学习\Python数据分析\Python数据分析与挖掘实战\Python数据分析与挖掘实战\chapter5\demo\data\bankloan.xls')data.head(3)年龄教育工龄地址收入负债率信用卡负债其他负债违约041317121769.311.3593925.008608112711063117.31.3622024.000798024011514555.50.8560752.1689250X=data.iloc[:,:-1]y=data.iloc[:,-1]X_train,X_test,y_train,y_test=train_test_split(X,y,stratify=y)fromsklearn.linear_modelimportLogisticRegressionCVlr=LogisticRegressionCV(Cs=[0.0001,0.001,0.01,0.1,1,10],cv=5).fit(X_train,y_train)print("训练精度:{}\n测试精度:{}\n总精度:{}".format(lr.score(X_train,y_train),lr.score(X_test,y_test),lr.score(X,y)))print("正则化系数:%s"%lr.C_)print("变量系数:%s"%lr.coef_)训练精度:0.8285714285714286测试精度:0.7771428571428571总精度:0.8157142857142857正则化系数:[10.]变量系数:[[0.029383150.08359976-0.24002552-0.08184292-0.004432810.100485780.58023866-0.02504046]]线性SVC见SVM
线性判别分析(LinearDiscriminantAnalysis,简称LDA)是一种监督学习的降维方法。
LDA的思想非常朴素。给定训练集,将样例投影到一条直线上,使得同类样例的投影点尽可能接近、异类样例的投影点尽可能远离;在对新样本进行分类时,将其投影到同样的这条直线上,再根据投影点的位置来确定新样本的类别。
一般来说,如果我们的数据是有类别标签的,那么优先选择LDA去尝试降维;当然也可以使用PCA做很小幅度的降维去消去噪声,然后再使用LDA降维。如果没有类别标签,那么肯定PCA是最先考虑的一个选择了。以下是一个使用LDA进行降维的例子:
frommpl_toolkits.mplot3dimportAxes3Dfromsklearn.datasets.samples_generatorimportmake_classificationX,y=make_classification(n_samples=1000,n_features=3,n_redundant=0,n_classes=3,n_informative=2,n_clusters_per_class=1,class_sep=0.5,random_state=10)fig=plt.figure()ax=Axes3D(fig,rect=[0,0,1,1],elev=30,azim=20)ax.scatter(X[:,0],X[:,1],X[:,2],marker='o',c=y)plt.show()
fromsklearn.discriminant_analysisimportLinearDiscriminantAnalysislda=LinearDiscriminantAnalysis(n_components=2).fit(X,y)X_trans=lda.transform(X)plt.scatter(X_trans[:,0],X_trans[:,1],marker='o',c=y)plt.show()
许多线性模型不适用于多分类问题,通常,我们基于一些基本策略,利用二分类学习器来解决多分类问题。多分类学习的基本思路是“拆解法”,即将多分类任务拆为若干个二分类任务求解。具体来说,先对问题进行拆分,然后为拆出的每个二分类任务训练一个分类器;在测试时,对这些分类器的预测结果进行集成以获得最终的多分类结果。有两点关键:
一对一(OvO):对样本进行N分类,将类别两两组合,共$C_N^2=\frac{N(N-1)}2$个二分类器,使用每个分类器对二分类器预测,对预测结果使用多数表决的方法选出最终结果。
一对其余(OVR):对样本进行N分类,对每个类别构造一个二分类器,共N个二分类器,二分类器的结果是“是”或者“不是”该类别。使用N个分类器对每个样本进行分类,被分类为正例的则为该样例的标签,若被预测为多个标签则选择置信度较高的标签作为预测结果。
多对多(MvM):每次将若干个类作为正类,若干个其他类作为反类。下面是一个多对多的常用方法:
纠错输出码(ErrorCorrectingOutputCodes,简称ECOC)工作分为两步:
如图所示,是两种编码方式。对于左侧的二元码,构造五个二分类器,也就是对数据集进行五次划分,如$f_1$将$C_2$类数据划分为正例,其他划分为反例。对测试样例使用五个分类器分别分类,得到五个分类标签,作为分类向量,计算分类向量与各个类别分类向量的距离,即可衡量测试样例与五个类别的差异程度,选择差异最小的类别作为最终分类类别。右侧分类器只是划分时分为三个类别。左边分类结果是$C_3$,右侧是$C_2$。
ECOC编码对分类器的错误有一定的容忍和修正能力。例如图(a)中对测试示例的正确预测编码是(-1,+1,+1,-1,+1),在预测时某个分类器出错了,例如f2出错从而导致了错误编码(-1,-1,+1,-1,+1),但基于这个编码仍能产生正确的最终分类结果$C_3$.一般来说对同一个学习任务,ECOC编码越长,纠错能力越强。但是,编码越长,意味着所需训练的分类器越多,计算、存储开销都会增大;另一方面,对有限类别数,可能的组合数目是有限的,码长超过一定范围后就失去了意义。
sklearn中,线性SVC使用的就是一对其余方法进行多分类,数据如下所示:
fromsklearn.datasetsimportmake_blobsX,y=make_blobs(random_state=42)mglearn.discrete_scatter(X[:,0],X[:,1],y)plt.xlabel('Feature0')plt.ylabel('Featrue1')plt.legend(['Class0','Class1','Class2'])plt.show()
对于以上数据使用线性SVC进行分类:
fromsklearn.svmimportLinearSVCLSVC=LinearSVC().fit(X,y)print("系数矩阵:{}\n截距:{}\n评分:{}".format(LSVC.coef_,LSVC.intercept_,LSVC.score(X,y)))系数矩阵:[[-0.174924820.23141171][0.47621845-0.06937277][-0.18913971-0.20400358]]截距:[-1.077455270.13140594-0.08604909]评分:1.0可以看出,LinearSVC输出了三条直线,下面将其可视化,三条直线将其分为7个区域,交叉区域平均分配。
mglearn.plots.plot_2d_classification(LSVC,X,fill=True,alpha=.7)mglearn.discrete_scatter(X[:,0],X[:,1],y)line=np.linspace(-10,10)forcoef,intercept,colorinzip(LSVC.coef_,LSVC.intercept_,['b','r','g']):plt.plot(line,-(line*coef[0]+intercept)/coef[1])plt.xlabel('Feature0')plt.ylabel('Featrue1')plt.legend(['Class0','Class1','Class2','Lineclass0','Lineclass1','Lineclass2'],loc=(1.01,0))
类别不平衡(classimbalance)就是指分类任务中不同类别的训练样例数目差别很大的情况。
比如,一个样本,正例2个,反例998个,只要分类器预测全部为反例,就能得到99.8%的精度,但这个精度是没有意义的。
再缩放:判定正例的临界条件改为$$\frac{y}{1-y}\times\frac{m^-}{m^+}>1$$
在使用分类器进行分类时,分类器给出一个预测值y,将y与阈值进行比较,如通常在$y>0.5$的时候预测正例,否则反例,y反映的是正例的可能性,设置阈值为0.5,则是认为正例与反例的可能性相等。
然而,当训练集中正、反例的数目不同时,令$m^+$表示正例数目,$m^-$表示反例数目,则观测几率是$\frac{m^+}{m^-}$;,由于我们通常假设训练集是真实样本总体的无偏采样,因此观测几率就代表了真实几率。
在实际中,针对类别不平衡问题通常有以下几种处理方法:
欠采样法若随机丢弃反例,可能丢失一些重要信息。欠采样法的代表性算法EasyEnsemble则是利用集成学习机制,将反例划分为若干个集合供不同学习器使用,这样对每个学习器来看都进行了欠采样,但在全局来看却不会丢失重要信息。
过采样法不能简单地对初始正例样本进行重复采样,否则会招致严重的过拟合,过采样法的代表性算法SMOTE是通过对训练集里的正例进行插值来产生额外的正例。