机器学习（5）——多项式回归与模型泛化

什么是多项式回归

import numpy as np
import matplotlib.pyplot as plt

x=np.random.uniform(-3,3,size=100)
X=x.reshape(-1,1)

y=0.5* x**2 +x+ 2 + np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

用线性回归试试：

from sklearn.linear_model import LinearRegression
lin_reg=LinearRegression()
lin_reg.fit(X,y)

y_predict=lin_reg.predict(X)

plt.scatter(x,y)
plt.plot(x,y_predict,color='r')
plt.show()

发现根本没有拟合，我们可以把x**2这一项设为另一个特征，和原来的特征合并后一起训练，这就是多项式回归：

X2=np.hstack([X,X**2])

lin_reg2=LinearRegression()
lin_reg2.fit(X2,y)
y_predict2=lin_reg2.predict(X2)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict2[np.argsort(x)],color='r') #要对x排序，否则是混乱的折线图
plt.show()

看看预测的特征值和截距：

lin_reg2.coef_

lin_reg2.intercept_

发现都很接近我们开始设置的。

scikit-learn中的多项式回归和Pipeline

scikit-learn中的多项式回归

from sklearn.preprocessing import PolynomialFeatures

poly=PolynomialFeatures(degree=2)
poly.fit(X)
X2=poly.transform(X)

X2.shape

X2[:5,:] #分别是x的0次方,x的1次方,x的二次方

from sklearn.linear_model import LinearRegression
lin_reg2=LinearRegression()
lin_reg2.fit(X2,y)
y_predict2=lin_reg2.predict(X2)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict2[np.argsort(x)],color='r')
plt.show()

x^0的系数是0，这是正确的，因为常数是截距。

X=np.arange(1,11).reshape(-1,2)
poly=PolynomialFeatures(degree=2)
poly.fit(X)
X2=poly.transform(X)

X2.shape
X2#如果原来有两个特征a、b，那么平方后是三项：a^2,a*b,b^2

那么degree=3时应该有10项。

Pipeline

Pipeline是scikit-learn提供的一大利器，可以顺序处理一些流程：

x=np.random.uniform(-3,3,size=100)
X=x.reshape(-1,1)
y=0.5* x**2 + x + 2 +np.random.normal(0,1,100)

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
poly_reg=Pipeline([
    ("poly",PolynomialFeatures(degree=2)),#最大为2次方
    ("std_scaler",StandardScaler()),#归一化
    ("lin_reg",LinearRegression()) #线性回归
])

poly_reg.fit(X,y)
y_predict=poly_reg.predict(X)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict[np.argsort(x)],color='r')
plt.show()

过拟合与欠拟合

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x=np.random.uniform(-3.0,3.0,size=100)
X=x.reshape(-1,1)
y=0.5* x**2 + x +2 +np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

from sklearn.linear_model import LinearRegression
lin_reg=LinearRegression()
lin_reg.fit(X,y)
lin_reg.score(X,y) #R^2值

下面我们使用均方误差：

直线：

from sklearn.metrics import mean_squared_error

y_predict=lin_reg.predict(X)
mean_squared_error(y,y_predict)

多项式：

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomialRegression(degree):
    return Pipeline([
        ("poly",PolynomialFeatures(degree=degree)),
        ("std_scaler",StandardScaler()),
        ("lin_reg",LinearRegression())
    ])

poly2_reg=PolynomialRegression(degree=2)
poly2_reg.fit(X,y)

y2_predict=poly2_reg.predict(X)
mean_squared_error(y,y2_predict)

测试一下扩展到100次方：

poly100_reg=PolynomialRegression(degree=100)
poly100_reg.fit(X,y)

y100_predict=poly100_reg.predict(X)
mean_squared_error(y,y100_predict)

plt.scatter(x,y)
plt.plot(np.sort(x),y100_predict[np.argsort(x)],'r')
plt.show()

用-3~3的等差数列来测试：

X_plot=np.linspace(-3,3,100).reshape(100,1)
y_plot=poly100_reg.predict(X_plot)
plt.scatter(x,y)
plt.plot(X_plot[:,0],y_plot,color="r")
plt.axis([-3,3,-1,10])
plt.show()

曲线变得很复杂，这就是过拟合；而直线太过简单，是欠拟合。

train_test_split的意义

可以发现这条过拟合的曲线面对新数据预测是非常弱的，我们称模型的泛化能力弱。

所以测试数据集的意义就是把数据拆分成训练数据和测试数据，如果训练得到的模型对测试数据准确率也比较高，那么泛化能力强。

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=666)

lin_reg=LinearRegression()
lin_reg.fit(X_train,y_train)
y_predict=lin_reg.predict(X_test)
mean_squared_error(y_test,y_predict)

对于线性回归，拆分成训练和测试数据集后：

对于多项式回归：

poly2_reg=PolynomialRegression(degree=2)
poly2_reg.fit(X_train,y_train)
y2_predict=poly2_reg.predict(X_test)
mean_squared_error(y_test,y2_predict)

传参为10的时候：

发现刚才把整个数据当训练数据的时候误差要比传参为2的时候还小，但拆分数据集后误差是比2次多项式要大的。

学习曲线

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x=np.random.uniform(-3.0,3.0,size=100)
X=x.reshape(-1,1)
y=0.5* x**2 + x + 2 + np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=10)

X_train.shape

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

train_score=[]
test_score=[]
for i in range (1,76):
    lin_reg=LinearRegression()
    lin_reg.fit(X_train[:i],y_train[:i])
    y_train_predict=lin_reg.predict(X_train[:i])
    train_score.append(mean_squared_error(y_train[:i],y_train_predict))
    
    y_test_predict=lin_reg.predict(X_test)
    test_score.append(mean_squared_error(y_test,y_test_predict))
    
plt.plot([i for i in range(1,76)],np.sqrt(train_score),label="train")
plt.plot([i for i in range(1,76)],np.sqrt(test_score),label="test")
plt.legend()
plt.show()

把绘制学习曲线封装成函数：

def plot_learning_curve(algo,X_train,X_test,y_train,y_test):
    train_score=[]
    test_score=[]
    for i in range (1,len(X_train)+1):
        algo.fit(X_train[:i],y_train[:i])
        y_train_predict=algo.predict(X_train[:i])
        train_score.append(mean_squared_error(y_train[:i],y_train_predict))
    
        y_test_predict=algo.predict(X_test)
        test_score.append(mean_squared_error(y_test,y_test_predict))
        
    plt.plot([i for i in range(1,len(X_train)+1)],np.sqrt(train_score),label="train")
    plt.plot([i for i in range(1,len(X_train)+1)],np.sqrt(test_score),label="test")
    plt.legend()
    plt.axis([0,len(X_train)+1,0,4])
    plt.show()

线性回归的学习曲线（欠拟合）：

plot_learning_curve(LinearRegression(),X_train,X_test,y_train,y_test)

二次多项式回归的学习曲线（最佳）：

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomialRegression(degree):
    return Pipeline([
        ("poly",PolynomialFeatures(degree=degree)),
        ("std_scaler",StandardScaler()),
        ("lin_reg",LinearRegression())
    ])

poly2_reg=PolynomialRegression(degree=2)

plot_learning_curve(poly2_reg,X_train,X_test,y_train,y_test)

20次多项式回归的学习曲线：

poly20_reg=PolynomialRegression(degree=20)

plot_learning_curve(poly20_reg,X_train,X_test,y_train,y_test)

可以发现，过拟合和欠拟合最终趋近于的高度（误差）是比最佳曲线要高的。

验证数据集与交叉验证

我们采用train_test_split的方法每次根据测试数据集的准确率来修改参数重新拟合，那么最后可能会对测试数据集过拟合了，可以再拆分出一块验证数据集来承担原来测试数据集的任务，而测试数据集是最终评判模型的标准：

交叉验证：

将训练数据拆成k个模型，K个模型的均值作为结果调参。

import numpy as np
from sklearn import datasets

digits=datasets.load_digits()
X=digits.data
y=digits.target

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test=train_test_split(X,y,test_size=0.4,random_state=666) 

#采用一般方法寻找最优超参数
from sklearn.neighbors import KNeighborsClassifier
best_score,best_p,best_k=0,0,0
for k in range(2,11):
    for p in range(1,6):
        knn_clf=KNeighborsClassifier(weights="distance",n_neighbors=k,p=p)
        knn_clf.fit(X_train,y_train)
        score=knn_clf.score(X_test,y_test)
        if score>best_score:
            best_score,best_p,best_k=score,p,k
        
print("Best k=",best_k)
print("Best p=",best_p)
print("Best Score=",best_score)

调用交叉验证：

from sklearn.model_selection import cross_val_score
knn_clf=KNeighborsClassifier()
cross_val_score(knn_clf,X_train,y_train)

默认是把训练数据拆成三分，如上图所示，交叉验证三次。

#使用交叉验证寻找最优超参数
best_score,best_p,best_k=0,0,0
for k in range(2,11):
    for p in range(1,6):
        knn_clf=KNeighborsClassifier(weights="distance",n_neighbors=k,p=p)
        scores=cross_val_score(knn_clf,X_train,y_train)
        score=np.mean(scores)
        if score>best_score:
            best_score,best_p,best_k=score,p,k
        
print("Best k=",best_k)
print("Best p=",best_p)
print("Best Score=",best_score)

best_knn_clf=KNeighborsClassifier(weights="distance",n_neighbors=2,p=2)

best_knn_clf.fit(X_train,y_train)
best_knn_clf.score(X_test,y_test)

虽然正确率比上面用一般方法（train_test_split）低了一点，但是上面的方法可能对测试数据集过拟合了，交叉验证的更可信一些。

回顾网格搜索：

from sklearn.model_selection import GridSearchCV
param_grid=[
    {
        'weights':['distance'],
        'n_neighbors':[i for i in range(2,11)],
        'p':[i for i in range(1,6)]
    }
]

grid_search=GridSearchCV(knn_clf,param_grid,verbose=1)
grid_search.fit(X_train,y_train)

其中的CV就是交叉验证，返回的结果意思是把训练数据集分成三块，共有45中模型，一共要交叉验证3*45次。

grid_search.best_score_
grid_search.best_params_
best_knn_clf=grid_search.best_estimator_
best_knn_clf.score(X_test,y_test)

可以看出和我们上面用交叉验证得到的结果一样。

把数据集分成5份：

cross_val_score(knn_clf,X_train,y_train,cv=5)

GridSearchCV(knn_clf,param_grid,verbose=1,cv=5)

k-folds交叉验证：把训练数据集分成k份，称为k-fold cross validation。缺点是每次训练k个模型，相当于整体性能慢了k倍。

留一法（LOO-CV）：把训练数据集分成m分，称为留一法（Leave-One-Out Cross Validation）。完全不受随机的影响，最接近模型真正的性能指标。缺点是计算量巨大。

偏差方差平衡

偏差（bias）、方差（variance）

模型误差=偏差+方差+不可避免的误差

导致偏差的主要原因：对问题本身的假设不正确，如非线性数据使用线性回归，欠拟合

方差：数据的一点点扰动都会较大影响模型，通常原因是使用的模型太复杂，如高阶多项式回归，过拟合。

有一些算法天生是高方差的算法，如KNN。

非参数学习通常都是高方差，因为不对数据进行任何假设。

有一些算法天生是高偏差算法，如线性回归。

参数学习通常都是高偏差算法，因为对数据具有极强的假设。

大多数算法具有相应的参数，可以调整偏差和方差。如KNN的K，线性回归中使用多项式回归。

偏差和方差通常是矛盾的，降低偏差会提高方差，降低方差会提高偏差。

机器学习的主要挑战，来自于方差！

解决高方差的通常手段：

1. 降低模型复杂度
2. 减少数据维度、降噪
3. 增加样本数
4. 使用验证集
5. 模型正则化

模型泛化与岭回归

模型正则化（Regularization）：限制参数的大小

过拟合的曲线非常陡峭，这是因为求得的参数大小规模相差较大

使用线性数据加噪音来测试一下：

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
x=np.random.uniform(-3.0,3.0,size=100)
X=x.reshape(-1,1)
y=0.5* x + 3 +np.random.normal(0,1,size=100)

plt.scatter(x,y)
plt.show()

多项式回归：

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

def PolynomialRegression(degree):
    return Pipeline([
        ("poly",PolynomialFeatures(degree=degree)),
        ("std_scaler",StandardScaler()),
        ("lin_reg",LinearRegression())
    ])

from sklearn.model_selection import train_test_split
np.random.seed(666)
X_train,X_test,y_train,y_test=train_test_split(X,y)

from sklearn.metrics import mean_squared_error
poly_reg=PolynomialRegression(degree=20)
poly_reg.fit(X_train,y_train)
y_poly_predict=poly_reg.predict(X_test)
mean_squared_error(y_test,y_poly_predict)
    

可见MSE是非常大的。

把模型的拟合曲线绘制封装成函数：

def plot_model(model):
    X_plot=np.linspace(-3,3,100).reshape(100,1)
    y_plot=model.predict(X_plot)
    plt.scatter(x,y)
    plt.plot(X_plot[:,0],y_plot,color="r")
    plt.axis([-3,3,0,6])
    plt.show()
    
plot_model(poly_reg)

很明显，20阶的多项式对训练数据过拟合了。

我们用岭回归来一发：

from sklearn.linear_model import Ridge
def RidgeRegression(degree,alpha):#alpha就是我们公式的阿尔法
    return Pipeline([
        ("poly",PolynomialFeatures(degree=degree)),
        ("std_scaler",StandardScaler()),
        ("ridge_reg",Ridge(alpha=alpha))
    ])
    
ridge1_reg=RidgeRegression(20,0.0001)
ridge1_reg.fit(X_train,y_train)
y1_predict=ridge1_reg.predict(X_test)
mean_squared_error(y_test,y1_predict)

可见MSE降了非常多！！

plot_model(ridge1_reg)

如果alpha=1：

ridge2_reg=RidgeRegression(20,1)
ridge2_reg.fit(X_train,y_train)
y2_predict=ridge2_reg.predict(X_test)
mean_squared_error(y_test,y2_predict)

plot_model(ridge2_reg)

alpha=100：

ridge3_reg=RidgeRegression(20,100)
ridge3_reg.fit(X_train,y_train)
y3_predict=ridge3_reg.predict(X_test)
mean_squared_error(y_test,y3_predict)

plot_model(ridge3_reg)

可以发现，随着alpha越来越大，损失函数中系数占比越来越大，那么重心就会转到降低每个系数的大小，最后系数趋近于0，就是直线了。

LASSO回归

继续使用上面岭回归的数据：

from sklearn.linear_model import Lasso

def LassoRegression(degree,alpha):
    return Pipeline([
        ("poly",PolynomialFeatures(degree=degree)),
        ("std_scaler",StandardScaler()),
        ("lasso_reg",Lasso(alpha=alpha))
    ])
    
lasso1_reg=LassoRegression(20,0.01)
lasso1_reg.fit(X_train,y_train)

y1_predict=lasso1_reg.predict(X_test)
mean_squared_error(y_test,y1_predict)

plot_model(lasso1_reg)

传入alpha=0.1：

lasso2_reg=LassoRegression(20,0.1)
lasso2_reg.fit(X_train,y_train)

y2_predict=lasso2_reg.predict(X_test)
mean_squared_error(y_test,y2_predict)

plot_model(lasso2_reg)

可见已经近似直线了，而岭回归是曲线。

比较Ridge和LASSO:

L1、L2正则项和弹性网络

L0正则就是让theta的数量尽量少。

弹性网综合了岭回归（计算相对是精准的，但特征比较多的时候因为不会特征选择所以计算量比较大）和lasso回归（特征选择，因为急于把特征化成0，所以可能会产生错误）的优点。