线性回归模型
目录：
1，最小二乘公式推导：
1.1，α、β推导
1.2，多项式回归
2，损失函数、正则化：
2.1，Ridge回归
2.2，LASSO回归
2.3，Elasitc Net算法
2.4，局部加权回归(Local Weight Regression)
3，模型评价指标
4，code：
4.1，波士顿房价预测：线性回归(无正则化项)、Ridge回归(L2正则)、LASSO回归(L1正则)、Elasitc Net算法(L1和L2正则)；
4.2，鲍鱼年龄预测：局部加权回归(样本残差平方加权)，不调包自建局部加权回归模型。
1，最小二乘公式推导：
线性模型的假设条件：
1），y的均值是x的线性组合（LinearFunction）；
2），残差e_i独立于x；
3），给定x, 残差e_i要服从正态分布（Normal Distribution）；
4），对于不同的xi, 残差e_i的方差variance应相同（qual Variance）。
1.1，α、β推导：
“二乘”指的是用平方来度量观测点与估计点的远近（在古汉语中“平方”称为“二乘”），“最小”指的是参数的估计值要保证各个观测点与估计点的距离的平方和达到最小。最小化残差平方和，即最小化SSE：
最终求得参数：
1.2，多项式回归：
对于多项式回归，即自变量多余一个，如下图所示：
此时可以以向量和矩阵形式表示多项式线性模型为：
此时，优化目标为最小化下面这个损失函数：
此时，对于多项式回归，最小二乘法的推导公式如下：
2，损失函数、正则化与Ridge回归、LASSO回归、Elasitc Net算法和局部加权回归：
为了克服模型过拟合的风险，加入正则化项，比如，L1正则、L2正则：
2.1，Ridge回归：
使用L2正则的线性回归模型就称为Ridge回归(岭回归)，即上图的第一个公式。
2.2，LASSO回归：
使用L1正则的线性回归模型就称为LASSO回归(Least Absolute Shrinkage and Selection Operator)，即上图的第二个公式。
上图直观地显示了Ridge回归于LASSO回归的区别。
2.3，Elasitc Net算法：
同时使用L1正则和L2正则的线性回归模型就称为Elasitc Net算法(弹性网络算法)，公式如下：
其中，λ>0; P∈[0,1].
2.4，局部加权回归：
局部加权回归(Local Weight Regression)损失函数：
w(i)是权重，它根据要预测的点与数据集中的点的距离来为数据集中的点赋权值。 当某点离要预测的点越远，其权重越小，否则越大。例如选择指数衰减函数(其中k为波长参数，它控制了权值随距离下降的速率)作为w公式：
使用该方式主要应用到样本之间的相似性考虑，实际上，LocalWeight Regression 与KNN都是为位置数据量身定制的，在局部进行训练。
局部加权回归（lowess）较好的解决了平滑问题。在做数据平滑的时候，会有遇到有趋势或者季节性的数据，对于这样的数据，我们不能使用简单的均值正负3倍标准差以外做异常值剔除，需要考虑到趋势性等条件。使用局部加权回归，可以拟合一条趋势线，将该线作为基线，偏离基线距离较远的则是真正的异常值点。实际上，局部加权回归(Lowess)主要还是处理平滑问题的多，因为预测问题，可以有更多模型做的更精确。但就平滑来说，Lowess很直观而且很有说服力。
3，模型评估：
回归算法的评价指标有：R-Squared、MSE，RMSE，MAE等。
4，code：
基于波士顿房价预测数据实现：1，线性回归(无正则化项)、2，Ridge回归(L2正则),3，LASSO回归(L1正则),4，Elasitc Net算法(L1和L2正则),5，局部加权回归(样本残差平方加权)，并不调包自建局部加权回归模型预测鲍鱼年龄。
代码语言：javascript复制import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.linear_model import LinearRegression, LassoCV, RidgeCV, ElasticNetCV
from sklearn.preprocessing import PolynomialFeatures #数据预处理，标准化
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV
import warnings
warnings.filterwarnings("ignore") # 拦截异常
# 1,加载数据：
def notEmpty(s):
return s != ''
names = ['CRIM','ZN', 'INDUS','CHAS','NOX','RM','AGE','DIS','RAD','TAX','PTRATIO','B','LSTAT']
df = pd.read_csv('../DataSets/boston_housing.data', header=None) # 数据文件格式不统一，所以读取时，先按照一行一个字段属性读取数据，然后再按照每行数据进行处理
data = np.empty((len(df), 14))
for i, d in enumerate(df.values): # enumerate生成一列索 引i,d为其元素
d = map(float, filter(notEmpty, d[0].split(' '))) # filter一个函数，一个list
data[i] = list(d) # 根据函数结果是否为真，来过滤list中的项
x, y = np.split(data, (13,), axis=1) # 分割数据
y = y.ravel() # 转换格式 拉直操作
ly=len(y)
print (y.shape, x[0:5], y[0:5])
print ("样本数据量:%d, 特征个数： {}, target样本数据量: {}".format(x.shape,y.shape[0]))
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=0) # 分割数据
# 2,Pipeline用于并行调参
models = [
Pipeline([
('ss', StandardScaler()),
('poly', PolynomialFeatures()),
('linear', LinearRegression())
]),
Pipeline([
('ss', StandardScaler()),
('poly', PolynomialFeatures()), # 特征工程：2次多项式特征（1,a,b,a^2,ab, b^2）
('linear', RidgeCV(alphas=np.logspace(-3,1,20))) # RidgeCV和Ridge的区别：前者可以进行交叉验证
]),
Pipeline([
('ss', StandardScaler()),
('poly', PolynomialFeatures()),
('linear', LassoCV(alphas=np.logspace(-3,1,20)))
]),
Pipeline([
('ss', StandardScaler()),
('poly', PolynomialFeatures()),
('linear', ElasticNetCV(alphas=np.logspace(-3,1,20), l1_ratio=[.1, .5, .7, .9, .95, 1]))
])
]
# 网格调参：参数字典，字典中的key是属性的名称，value是可选的参数列表
parameters = {
"poly__degree": [3,2,1],
"poly__include_bias": [True, False], #多项式幂为零的特征作为线性模型中的截距
"poly__interaction_only": [True, False], #不产生交互项，如X1*X1
"linear__fit_intercept": [True, False]
}
# 3,线性回归、Lasso回归、Ridge回归、ElasticNet比较：
colors = ['g-', 'b-', 'y-', 'k-']
titles = [u'Linear Regression', u'Ridge Regression', u'Lasso Regression', u'ElasticNet Regression']
plt.figure(figsize=(16,8), facecolor='w')
ln_x_test = range(len(x_test))
plt.plot(ln_x_test, y_test, 'r-', lw=1.5, label=u'fact')
import time
start = time.time()
# 训练模型：
for t in range(4):
model = GridSearchCV(models[t], param_grid=parameters,cv=5, n_jobs=1) # 网格搜索, 五折交叉验证
model.fit(x_train, y_train)
print ("%s算法:最优参数:" % titles[t],model.best_params_) # 模型效果值获取（最优参数）
print ("%s算法:R值=%.3f" % (titles[t], model.best_score_))
y_predict = model.predict(x_test) # 模型预测
# 画图
plt.plot(ln_x_test, y_predict, colors[t], lw = t*0.5 2, label=u'%s:,$R^2$=%.3f' % (titles[t], model.best_score_))
end = time.time()
print('cost time: ', end-start)
# 图形显示
plt.legend(loc = 'upper left')
plt.grid(True)
plt.title(u"boston_housing price-prediction")
plt.show()
# 4,局部加权回归(Local Weight Regression)：预测鲍鱼年龄
# 4.1，对象化封装功能函数：
import numpy as np
# =============================================================================
# process data
# =============================================================================
def loadDataSet(fileName): #general function to parse tab -delimited floats
numFeat = len(open(fileName).readline().split('t')) - 1 #get number of fields 2
dataMat = []; labelMat = []
fr = open(fileName)
for line in fr.readlines():
lineArr =[]
curLine = line.strip().split('t') #['1.000000', '0.116163', '3.129283']
for i in range(numFeat):
lineArr.append(float(curLine[i]))
dataMat.append(lineArr) #[[1.0, 0.52707], [1.0, 0.116163],...]
labelMat.append(float(curLine[-1])) #[4.225236, 4.231083,...]
return dataMat,labelMat #list
# =============================================================================
# linear regression
# =============================================================================
def standRegres(xArr,yArr):
xMat = np.mat(xArr) #转为矩阵
yMat = np.mat(yArr).T #转为矩阵再转置
xTx = xMat.T*xMat #xMat.T*xMat*w - xMat.T*yMat = 0
if np.linalg.det(xTx) == 0.0:
print ("This matrix is singular, cannot do inverse")
return
ws = xTx.I * (xMat.T*yMat)
return ws
# =============================================================================
# locally weighted linear regression
# =============================================================================
def lwlr(testPoint,xArr,yArr,k=1.0):
xMat = np.mat(xArr); yMat = np.mat(yArr).T
m = np.shape(xMat)[0]
weights = np.mat(np.eye((m)))
for j in range(m): #next 2 lines create weights matrix
diffMat = testPoint - xMat[j,:] #difference matrix
weights[j,j] = np.exp(diffMat*diffMat.T/(-2.0*k**2)) #weighted matrix
xTx = xMat.T * (weights * xMat)
if np.linalg.det(xTx) == 0.0:
print ("This matrix is singular, cannot do inverse")
return
ws = xTx.I * (xMat.T * (weights * yMat)) #normal equation
#ws #7 feature,and 1
return testPoint * ws
def lwlrTest(testArr,xArr,yArr,k=1.0): #loops over all the data points and applies lwlr to each one
m = np.shape(testArr)[0] # #trianing sample
yHat = np.zeros(m) #array
for i in range(m):
yHat[i] = lwlr(testArr[i],xArr,yArr,k)
return yHat
#others
def rssError(yArr,yHatArr):
return ((yArr-yHatArr)**2).sum()
# 4.2加载数据，并开始训练：
abX,abY = loadDataSet('../DataSets/abalone.txt')
#locally weighted linear regression
yHat01 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],0.1) #training set #0-99
yHat1 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],1)
yHat10 = lwlrTest(abX[0:99],abX[0:99],abY[0:99],10)
#error
error01 = rssError(abY[0:99],yHat01.T) #56.820227823572182
error1 = rssError(abY[0:99],yHat1.T) #429.89056187016683
error10 = rssError(abY[0:99],yHat10.T) #549.1181708825128
#generalization
yHat01g = lwlrTest(abX[100:199],abX[100:199],abY[100:199],0.1) #test set #100-199
yHat1g = lwlrTest(abX[100:199],abX[100:199],abY[100:199],1)
yHat10g = lwlrTest(abX[100:199],abX[100:199],abY[100:199],10)
#error
error01g = rssError(abY[100:199],yHat01g.T) #36199.797699875046
error1g = rssError(abY[100:199],yHat1g.T) #231.81344796874004
error10g = rssError(abY[100:199],yHat10g.T) #291.87996390562728 #compare：k = 0.1 overfitting
ws = standRegres(abX[0:99],abY[0:99]) # 线性回归
yHat = np.mat(abX[100:199])*ws
errorlr = rssError(yHat.T.A,abY[100:199]) #518.63631532510897
print('error10g < errorlr: ', error10g, errorlr) #compare: lwlr is better than lr
