惩罚线性回归---Python机器学习:预测分析核心算法

惩罚线性回归

参考教材:Python机器学习预测分析核心算法，书中代码较为过时，借用sklearn等工具包进行了重写。

实践中遇到的绝大多数预测分析（函数逼近）问题，惩罚线性回归和集成方法都具有最优或接近最优的性能。这些问题包含：大数据集、小数据集、宽数据集（wide data sets）、高瘦数据集（tall skinny data sets）、复杂问题、简单问题

惩罚线性回归模型一个重要优势就是它训练所需时间。当面对大规模的数据时，训练所需时间就成为一个需要考量的因素

---

1. 声呐数据集

声纳信号又叫作啁啾信号（chirpedsignal），即信号在一个脉冲期间频率会增加或降低。此数据集的测量值代表声纳接收器在不同地点接收到的返回信号，其中在大约一半的例子中，返回的声纳信号反映的是岩石的形状，而另一半是金属圆筒的形状（水雷）

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from random import uniform
from math import exp

%matplotlib inline

# 显示中文
plt.rcParams['font.sans-serif'] = [u'SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 加载数据
rocksVMines = pd.read_csv('./sonar.txt',header=None, prefix="V")
rocksVMines.head()

	V0	V1	V2	V3	V4	V5	V6	V7	V8	V9	...	V51	V52	V53	V54	V55	V56	V57	V58	V59	V60
0	0.0200	0.0371	0.0428	0.0207	0.0954	0.0986	0.1539	0.1601	0.3109	0.2111	...	0.0027	0.0065	0.0159	0.0072	0.0167	0.0180	0.0084	0.0090	0.0032	R
1	0.0453	0.0523	0.0843	0.0689	0.1183	0.2583	0.2156	0.3481	0.3337	0.2872	...	0.0084	0.0089	0.0048	0.0094	0.0191	0.0140	0.0049	0.0052	0.0044	R
2	0.0262	0.0582	0.1099	0.1083	0.0974	0.2280	0.2431	0.3771	0.5598	0.6194	...	0.0232	0.0166	0.0095	0.0180	0.0244	0.0316	0.0164	0.0095	0.0078	R
3	0.0100	0.0171	0.0623	0.0205	0.0205	0.0368	0.1098	0.1276	0.0598	0.1264	...	0.0121	0.0036	0.0150	0.0085	0.0073	0.0050	0.0044	0.0040	0.0117	R
4	0.0762	0.0666	0.0481	0.0394	0.0590	0.0649	0.1209	0.2467	0.3564	0.4459	...	0.0031	0.0054	0.0105	0.0110	0.0015	0.0072	0.0048	0.0107	0.0094	R

5 rows × 61 columns

用分位数图展示异常点

分位数图展示了数据的百分位边界与高斯分布的同样百分位的边界对比。如果此数据服从高斯分布，则画出来的点应该是一条直线。

from scipy import stats

stats.probplot(rocksVMines.iloc[:,3].values,dist='norm',plot=plt)
plt.xlabel(u'分位数')
plt.title(u'概率图')
plt.ylabel('已经排序数值')
plt.show()

观察可知:声呐数据集尾部的数据要多于高斯分布尾部的数据

利用平行坐标图进行可视化

对于具有多个属性问题的一种可视化方法叫平行坐标图，即绘制每行数据由属性值构成的行向量。

整个数据集的平行坐标图对于数据集中的每一行属性都有对应的一条折线。基于标签对折线标示不同的颜色，更有利于观测到属性值与标签之间的关系。

for i in range(rocksVMines.shape[0]):
    if rocksVMines.loc[i,'V60'] == 'M':
        pcolor = 'red'
    else:
        pcolor = 'blue'
    rocksVMines.drop('V60',axis=1).iloc[i].plot(color=pcolor)

plt.xlabel("属性")
plt.ylabel("属性值")
plt.show()

看不出明显的差别

1. 沿图的底部，蓝色的线要突出一点；
2. 在30-40之间蓝色的线多少要比红色的线要高点。

这些观察有助于解释和确认某些预测的结果

---

属性和标签的关系可视化

1. 绘制属性之间的散点图（以2-3， 2-21属性对为例）

下面绘制的两张散点图直观显示，频率离的近的信号相关性要更高，后面的相关系数计算更加说明了这一点。

plt.scatter(rocksVMines.iloc[:,1],rocksVMines.iloc[:,2])
plt.xlabel("第2个属性")
plt.ylabel("第3个属性")
plt.show()


plt.scatter(rocksVMines.iloc[:,1],rocksVMines.iloc[:,21])
plt.xlabel("第2个属性")
plt.ylabel("第21个属性")
plt.show()

print("第2属性和第3属性的相关系数：",rocksVMines.iloc[:,[1,2]].corr().values[0,1])
print("第2属性和第21属性的相关系数：",rocksVMines.iloc[:,[1,20]].corr().values[0,1])

第2属性和第3属性的相关系数： 0.7799158719104263
第2属性和第21属性的相关系数： 0.07525528679117005

2） 绘制区分度最大的属性（猜测）和标签的散点图

根据前面的平行坐标图，在30-40之间蓝色的线多少要比红色的线要高点，并且大概在35索引值附近有所分离。

于是我们猜测第35个属性区分度最大，并绘制这个属性和标签的散点图。

下面第一张图的散点挤在一起看不清楚，于是故意将标签值上下随机做点扰动，并且进行半透明处理，于是有了第二张图。

plt.scatter(rocksVMines.iloc[:,35],rocksVMines.iloc[:,60].map(lambda x: 1.0 if x == 'M' else 0.0))
plt.xlabel("属性值")
plt.ylabel("标签值")
plt.show()

# 扰动
plt.scatter(rocksVMines.iloc[:,35],rocksVMines.iloc[:,60].map(lambda x: (1.0 if x == 'M' else 0.0)+uniform(-0.1, 0.1)), alpha=0.5, s=60)
plt.xlabel("属性值")
plt.ylabel("标签值")
plt.show()

注意到第35个属性在左上方的点更加集中一些，然而在下面的数据从右到左分布得更加均匀些。上方的数据对应水雷的数据。下面的数据对应岩石的数据。由图观察可知，可以因此建立一个分类器，判断第35个属性是否大于或小于0.5。如果大于0.5，就判断为岩石，如果小于0.5，就判断为水雷。这样就可以获得一个比随机猜测好些的结果。

---

用热度图展示属性和标签的相关性

多个属性间的相关性很高（相关系数＞0.7），即多重共线性（multicollinearity），往往会导致预测结果不稳定

rocksVMines['V60'] = rocksVMines['V60'].map(lambda x: 1.0 if x == 'M' else 0.0)

plt.pcolor(rocksVMines.corr())
plt.xlim(0, 61)
plt.ylim(0, 61)
plt.show()

ROC曲线

ROC曲线就是使用一个图来展示不同的混淆矩阵。实际上，是绘制真正率(TPR)随假正率(FPR)的变化情况。

$$TRP=\frac{TP}{TP+FN}$$

$$FRP=\frac{FP}{TN+FP}$$

ROC曲线越接近左上角，效果越好。 如果ROC曲线掉到对角线下面，很可能是把预测方向搞反了。

from sklearn import datasets, linear_model
from sklearn.metrics import roc_curve, auc

# 混淆矩阵
def confusion_matrix(predicted, actual, threshold):
    if len(predicted) != len(actual):
        return -1
    tp, fp, tn, fn = 0.0, 0.0, 0.0, 0.0
    for i in range(len(actual)):
        if actual[i] > 0.5:
            if predicted[i] > threshold:
                tp += 1.0
            else:
                fn += 1.0
        else:
            if predicted[i] < threshold:
                tn += 1.0
            else:
                fp += 1.0
    rtn = [tp, fn, fp, tn]
    return rtn

X = rocksVMines.iloc[:,:-1]
labels =rocksVMines.iloc[:,-1]
X.shape,labels.shape

((208, 60), (208,))

# 划分数据集
indices = range(X.shape[0])
X_test = X[X.index %3 ==0].values
X_train= X[X.index %3 !=0].values

y_test = labels[labels.index%3 == 0].values
y_train= labels[labels.index%3 != 0].values

X_test.shape,y_test.shape,X_train.shape,y_train.shape

((70, 60), (70,), (138, 60), (138,))

# 建模
model = linear_model.LinearRegression()
model.fit(X_train,y_train)

pre_X = model.predict(X_train)
pre_X[:5]

array([-0.10240253,  0.42090698,  0.38593034,  0.36094537,  0.31520494])

# 计算训练集的混淆矩阵
threshold = 0.5
confusion_mat_train = confusion_matrix(pre_X,y_train,threshold)

confusion_mat_train

[68.0, 6.0, 7.0, 57.0]

# 计算测试集的混淆矩阵
pre_test = model.predict(X_test)
confusion_mat_test = confusion_matrix(pre_test,y_test,threshold)
confusion_mat_test

[28.0, 9.0, 9.0, 24.0]

# 计算训练集auc
fpr,tpr,thresholds = roc_curve(y_train,pre_X)
roc_auc = auc(fpr,tpr)
roc_auc

0.9795185810810811

# Plot ROC curve
plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('In sample ROC rocks versus mines')
plt.legend(loc="lower right")
plt.show()

# 计算测试集auc
fpr, tpr, thresholds = roc_curve(y_test,pre_test)
roc_auc = auc(fpr, tpr)
print( 'AUC for out-of-sample ROC curve: %f' % roc_auc)

AUC for out-of-sample ROC curve: 0.848485

plt.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Out-of-sample ROC rocks versus mines')
plt.legend(loc="lower right")
plt.show()

最优alpha

alphas = [-3, -2, -1, 0,1, 2, 3, 4, 5]
alphaList = [0.1**i for i in alphas]

aucList = []
for alph in alphaList:
    rocksVMinesRidgeModel = linear_model.Ridge(alpha=alph)
    rocksVMinesRidgeModel.fit(X_train,y_train)
    fpr, tpr, thresholds = roc_curve(y_test,rocksVMinesRidgeModel.predict(X_test))
    roc_auc = auc(fpr, tpr)
    aucList.append(roc_auc)

aucList

[0.8411138411138411,
 0.864045864045864,
 0.9074529074529074,
 0.9180999180999181,
 0.8828828828828829,
 0.8615888615888616,
 0.8517608517608517,
 0.8509418509418509,
 0.8493038493038493]

plt.plot(alphas, aucList)
plt.xlabel('-log(alpha)')
plt.ylabel('AUC')
plt.show()

indexBest = aucList.index(max(aucList))
alph = alphaList[indexBest]
rocksVMinesRidgeModel = linear_model.Ridge(alpha=alph)
rocksVMinesRidgeModel.fit(X_train, y_train)

#scatter plot of actual vs predicted
plt.scatter(rocksVMinesRidgeModel.predict(X_test), y_test, s=100, alpha=0.25)
plt.xlabel("Predicted Value")
plt.ylabel("Actual Value")
plt.show()

2. 鲍鱼数据集

根据某些测量值预测鲍鱼年龄

abalone = pd.read_csv('./abalone.txt',header=None,prefix='V')
# abalone.columns =  ['Sex', 'Length', 'Diameter', 'Height',
#                    'Whole Wt', 'Shucked Wt','Viscera Wt', 
#                    'Shell Wt', 'Rings']


abalone.columns = ['性别','长度','直径','高度','整体重量', '去壳重量', '脏器重量', '壳的重量', '环数']

abalone

	性别	长度	直径	高度	整体重量	去壳重量	脏器重量	壳的重量	环数
0	M	0.455	0.365	0.095	0.5140	0.2245	0.1010	0.1500	15
1	M	0.350	0.265	0.090	0.2255	0.0995	0.0485	0.0700	7
2	F	0.530	0.420	0.135	0.6770	0.2565	0.1415	0.2100	9
3	M	0.440	0.365	0.125	0.5160	0.2155	0.1140	0.1550	10
4	I	0.330	0.255	0.080	0.2050	0.0895	0.0395	0.0550	7
...	...	...	...	...	...	...	...	...	...
4172	F	0.565	0.450	0.165	0.8870	0.3700	0.2390	0.2490	11
4173	M	0.590	0.440	0.135	0.9660	0.4390	0.2145	0.2605	10
4174	M	0.600	0.475	0.205	1.1760	0.5255	0.2875	0.3080	9
4175	F	0.625	0.485	0.150	1.0945	0.5310	0.2610	0.2960	10
4176	M	0.710	0.555	0.195	1.9485	0.9455	0.3765	0.4950	12

4177 rows × 9 columns

查看热力图

多个属性间的相关性很高（相关系数＞0.7），即多重共线性（multicollinearity），往往会导致预测结果不稳定

corMat = abalone.corr()
plt.pcolor(corMat)
plt.show()

summary = abalone.describe()
summary

	长度	直径	高度	整体重量	去壳重量	脏器重量	壳的重量	环数
count	4177.000000	4177.000000	4177.000000	4177.000000	4177.000000	4177.000000	4177.000000	4177.000000
mean	0.523992	0.407881	0.139516	0.828742	0.359367	0.180594	0.238831	9.933684
std	0.120093	0.099240	0.041827	0.490389	0.221963	0.109614	0.139203	3.224169
min	0.075000	0.055000	0.000000	0.002000	0.001000	0.000500	0.001500	1.000000
25%	0.450000	0.350000	0.115000	0.441500	0.186000	0.093500	0.130000	8.000000
50%	0.545000	0.425000	0.140000	0.799500	0.336000	0.171000	0.234000	9.000000
75%	0.615000	0.480000	0.165000	1.153000	0.502000	0.253000	0.329000	11.000000
max	0.815000	0.650000	1.130000	2.825500	1.488000	0.760000	1.005000	29.000000

minRings = summary.loc['min']['环数']
maxRings = summary.loc['max']['环数']
nrows = abalone.shape[0]

minRings,maxRings,nrows

(1.0, 29.0, 4177)

abalone.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 4177 entries, 0 to 4176
Data columns (total 9 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   性别      4177 non-null   object 
 1   长度      4177 non-null   float64
 2   直径      4177 non-null   float64
 3   高度      4177 non-null   float64
 4   整体重量    4177 non-null   float64
 5   去壳重量    4177 non-null   float64
 6   脏器重量    4177 non-null   float64
 7   壳的重量    4177 non-null   float64
 8   环数      4177 non-null   int64  
dtypes: float64(7), int64(1), object(1)
memory usage: 293.8+ KB

平行坐标图

折线代表了一行数据，折线的颜色表明了其所属的类别。

鲍鱼问题是一个回归问题，应该用不同的颜色来对应标签值的高低。也就是实现由标签的实数值到颜色值的映射，需要将标签的实数值压缩到[0.0,1.0]区间。

可视化结果发现，属性值相近的地方，折线的颜色也很接近。这可以暗示可以构建相当准确的预测模型。

for i in range(nrows):
    dataRow = abalone.iloc[i].drop(['性别','环数'],axis=0)
    labelColor = (abalone.iloc[i]['环数']-minRings)/(maxRings-minRings)
    # 通过控制每一次循环计算的labelColor的值控制线条显示的色差，选择色环RdYlBu上面的不同颜色
    dataRow.plot(color=plt.cm.RdYlBu(labelColor),alpha=0.5)

plt.ylabel(("属性值"))
plt.show()

meanRings = summary.loc['mean']['环数']
stdRings = summary.loc['std']['环数']
meanRings,stdRings

(9.933684462532918, 3.2241690320681133)

for i in range(nrows):
    dataRow = abalone.iloc[i].drop(['性别','环数'],axis=0)
    normTarget = (abalone.iloc[i]['环数']-meanRings)/stdRings
    # 让取负值的数据与取正值的数据基本上一样多
    labelColor = 1.0/(1.0 + exp(-normTarget))
    dataRow.plot(color=plt.cm.RdYlBu(labelColor),alpha=0.5)

plt.ylabel(("属性值"))
plt.show()

箱线图

箱线图是一种识别异常点的强大工具：

• 中位数：中间红线；
• 25%分位数：箱顶；
• 75%分位数：箱底；
• 上边缘水平线： 到箱顶的距离是箱高的1.4倍（默认参数）, 即4分位间距的1.4倍
• 下边缘水平线： 到箱底的距离是箱高的1.4倍（默认参数）, 即4分位间距的1.4倍
• 异常值： 上下边缘以外的数据。

第一张图由于环数远大于其他指标，看不清楚，剔除后得到第二张图，将第二张图进行标准化于是有第三张图。

abalone.boxplot()
plt.ylabel(("四分位间距"))
plt.show()

abalone.drop('环数',axis=1).boxplot()
plt.ylabel(("四分位间距"))
plt.show()

new_abalone = abalone.drop(['性别','环数'],axis=1)
((new_abalone-new_abalone.mean())/new_abalone.std()).boxplot()

plt.ylabel("四分位间距 - 标准化")
plt.show()

最佳子集选择算法

在列的个数上增加一个约束（假设为nCol），然后从X的所有列中抽取特定个数的列构成数据集，在上边执行最小二乘法，遍历所有列的组合（列数为nCol），找到在测试集上取得最佳效果的nCol值；增加nCol值，重复上述过程。以上过程产生最佳的一列子集、两列子集一直到所有列子集（对应矩阵X）。对于每个子集同样有一个性能与之对应。下一步决定在部署时是使用一列子集版本、两列子集版本，还是其他版本，直接选择错误率最低的版本。

最佳子集选择存在的一个问题是该算法需要大量计算，即使属性不多的情况下（属性数对应X的列数），计算量也非常巨大。例如，10个属性对应于210=1000个子集。有几种方法可以避免这种情况，其中一种就是前向逐步回归。

前向逐步回归算法

前向逐步回归的想法是从1列子集开始，找到效果最佳的那一列属性，接着寻找与其组合与效果最佳的第2列属性，而不是评估所有的2列子集

基本步骤:

• 将β的所有值初始化为0。

在每一步中

• 使用已经选择的变量找到残差值。
• 确定哪个未使用的变量能够最佳的解释残差，将该变量加入选择变量中。

LARS算法与前向逐步回归算法非常类似。LARS与前向逐步回归算法的主要差异是LARS在引入新属性时只是部分引入，引入属性过程并非不可逆

from sklearn import datasets, linear_model
from sklearn.model_selection import train_test_split
from sklearn import datasets, linear_model
from math import sqrt          

wine_data = pd.read_csv('./winequality-red.csv',sep=';')
wine_data.head()

	fixed acidity	volatile acidity	citric acid	residual sugar	chlorides	free sulfur dioxide	total sulfur dioxide	density	pH	sulphates	alcohol	quality
0	7.4	0.70	0.00	1.9	0.076	11.0	34.0	0.9978	3.51	0.56	9.4	5
1	7.8	0.88	0.00	2.6	0.098	25.0	67.0	0.9968	3.20	0.68	9.8	5
2	7.8	0.76	0.04	2.3	0.092	15.0	54.0	0.9970	3.26	0.65	9.8	5
3	11.2	0.28	0.56	1.9	0.075	17.0	60.0	0.9980	3.16	0.58	9.8	6
4	7.4	0.70	0.00	1.9	0.076	11.0	34.0	0.9978	3.51	0.56	9.4	5

namelist = np.array(['非挥发性酸','挥发性酸'