数据分析师-机器学习之线性回归

回归分析：
回归：统计学分析数据的方法，目的在于了解两个或多个变量间是否相关、研究其相关方向与强度，并建立数学模型以便观察特定变量来预测研究者感兴趣的变量，回归分析可以帮助人们了解在自变量变化时因变量的变化量。一般来说，通过回归分析我们可以由给出的自变量估计因变量的条件期望。

误差分布：
真实值和预测值之间的差异：

误差分布特点（独立、同分布、高斯分布）。

线性回归的最小二乘法：

最小二乘法回归存在问题：
矩阵必须可逆。
当特征数较多时，求逆运算时间开销较大。

线性回归的梯度下降法：

梯度下降法（Gradient descent）是使得代价函数达到最小的经典方法之一。

梯度下降的算法调优：
算法的步长选择。
算法参数的初始值选择。
归一化。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline
# y = kx + b,k is slope, b is y-intercept
def step_gradient(b_current, k_current, points, learningRate):b_gradient = 0k_gradient = 0for i in range(0, len(points)):x = points[i, 0]y = points[i, 1]k_gradient += (((k_current * x) + b_current)-y)*xb_gradient += (((k_current * x) + b_current)-y)*1new_b = b_current - (learningRate * b_gradient)new_k = k_current - (learningRate * k_gradient)return [new_b, new_k]
#误差函数
def compute_error_for_line_given_points(b, k, points):totalError = 0for i in range(0, len(points)):x = points[i, 0]y = points[i, 1]totalError += (y - (k * x + b)) ** 2return totalError
def gradient_descent_runner(points, starting_b, starting_k, learning_rate, num_iterations):b = starting_bk = starting_kerror=[]for i in range(num_iterations):b, k = step_gradient(b, k, np.array(points), learning_rate)error.append(compute_error_for_line_given_points(b, k, points))return [b, k,error]#画原始数据图
points = np.array(pd.read_table("data.csv", delimiter=",",names=['x1','y']))
plt.figure(figsize=(8,4),dpi=100)
plt.scatter(points[:,0],points[:,1],color='black')
plt.title("original_data")
#梯度下降法
plt.figure(figsize=(8,4),dpi=100)
plt.scatter(points[:,0],points[:,1],color='black')
learning_rate = 0.000001 #学习率
initial_b = 0 # initial y-intercept guess
initial_k = 0 # initial slope guess
num_iterations = 10000
[b, k, error] = gradient_descent_runner(points, initial_b, initial_k, learning_rate, num_iterations)
y_predict=points[:,0]*k+b
plt.plot(points[:,0],y_predict,color='blue')
plt.title("GradientDescent")
#迭代次数与误差的关系
plt.figure(figsize=(8,4),dpi=100)
plt.plot(range(num_iterations)[0:100],np.array(error)[0:100],color='red')
plt.xlabel("iteration number") #添加x轴的名称
plt.ylabel("error")

梯度下降法分类：
批量梯度下降算法，每一次迭代使用全部的样本，能达到全局最优解（凸函数情况下）当样本数目很多时，训练过程缓。
随机梯度下降算法，每一次更新参数只使用一个样本，进行多次更新，迭代速度快，准确度下降，每次不一定朝着收敛的方向，容易陷入局部最优，非凸函数情况下有可能跳出局部最优。
小批量梯度下降算法，更新每一参数时都使用一部分样本来进行更新，相对于随机梯度下降，mini-batch梯度下降降低了收敛波动性，即降低了参数更新的方差，使得更新更加稳定。相对于批量梯度下降，其提高了每次学习的速度。

线性回归的岭回归：

岭回归是在平方误差的基础上增加正则项（L2正则）：

其中L2正则项系数>0,通过确定L2正则项系数的值可以在模型方差（回归系数的方差）和偏差之间达到平衡，随着的增大，模型方差减小而偏差增大。
与线性回归一样，利用最小二乘法求解岭回归模型参数时，可得：

岭回归的模型选择，特征与L2正则项系数的关系图；

#岭回归系数选择
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
from math import exp
def load_data():data=pd.read_table('abalone.txt',sep=' ',names=['x1','x2','x3','x4','x5','x6','x7','x8','y'])feature=data.loc[:,'x1':'x8']label=data['y']return np.mat(feature), np.mat(label).T
def standarize(X):''' 标准化数据 (零均值, 单位标准差)'''std_deviation = np.std(X)mean = np.mean(X)return (X - mean)/std_deviation
def ridge_regression(X, y, lambd=0.2):''' 获取岭回归系数'''m,n = X.shapeI = np.matrix(np.eye(n))w = (X.T*X + lambd*I).I*X.T*yreturn w
def ridge_traj(X, y, k=20):''' 获取岭轨迹矩阵'''m, n = X.shapews = np.zeros((k, n))for i in range(k):w = ridge_regression(X, y, lambd=exp(i-10))ws[i, :] = w.Treturn ws
if '__main__' == __name__:# 加载数据X, y = load_data()# 标准化X, y = standarize(X), standarize(y)# 绘制岭轨迹k = 30ws = ridge_traj(X, y, k)lambdas = [exp(i-10) for i in range(k)]plt.semilogx(lambdas, ws)plt.show()

线性回归的Lasso回归：

Lasso回归模型求解方法：坐标轴下降法是沿着坐标轴的方向去下降，是在平方误差的基础上增加L1正则：

其中L1正则项系数>0 。通过确定L1正则项系数值可以使得在方差和偏差之间达到平衡，随着L1正则项系数值的增大，模型方差减小而偏差增大。

lasso回归的模型选择，特征与L1正则项系数的关系图；（lasso相当于降维的功能）

#lasso回归
import itertools
from math import exp
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import os
os.chdir('C:\\Users\\chen\\File')
def load_data():''' 加载数据'''data=pd.read_table('abalone.txt',sep=' ',names=['x1','x2','x3','x4','x5','x6','x7','x8','y'])feature=data.loc[:,'x1':'x8']label=data['y']return np.mat(feature), np.mat(label).T
def standarize(X):''' 标准化数据 (零均值, 单位标准差)'''std_deviation = np.std(X)mean = np.mean(X)return (X - mean)/std_deviation
def lasso_regression(X, y, lambd=0.2, threshold=0.1):''' 通过坐标下降(coordinate descent)法获取LASSO回归系数'''# 计算残差平方和rss = lambda X, y, w: (y - X*w).T*(y - X*w)# 初始化回归系数w.m, n = X.shapew = np.matrix(np.zeros((n, 1)))r = rss(X, y, w)# 使用坐标下降法优化回归系数wniter = itertools.count(1)for it in niter:for k in range(n):# 计算常量值z_k和p_kz_k = (X[:, k].T*X[:, k])[0, 0]p_k = 0for i in range(m):p_k += X[i, k]*(y[i, 0] - sum([X[i, j]*w[j, 0] for j in range(n) if j != k]))if p_k < -lambd/2:w_k = (p_k + lambd/2)/z_kelif p_k > lambd/2:w_k = (p_k - lambd/2)/z_kelse:w_k = 0w[k, 0] = w_kr_prime = rss(X, y, w)delta = abs(r_prime - r)[0, 0]r = r_prime#print('Iteration: {}, delta = {}'.format(it, delta))if delta < threshold:breakreturn w
def lasso_traj(X, y, k=20):''' 获取回归系数轨迹矩阵'''m, n = X.shapews = np.zeros((k, n))for i in range(k):w = lasso_regression(X, y, lambd=exp(i-10))ws[i, :] = w.T#print('lambda = e^({}), w = {}'.format(i-10, w.T[0, :]))return ws
if '__main__' == __name__:X, y = load_data()X, y = standarize(X), standarize(y)w = lasso_regression(X, y, lambd=10)# 绘制轨迹k = 20ws = lasso_traj(X, y, k)fig = plt.figure()lambdas = [exp(i-10) for i in range(k)]plt.semilogx(lambdas, ws)plt.show()

坐标轴下降法和梯度下降比较：
1.坐标轴下降法在每次迭代中在当前点处沿一个坐标方向进行一维搜索，固定其他的坐标方向，找到一个函数的局部极小值。而梯度下降总是沿着梯度的负方向求函数的局部最小值。
2.坐标轴下降优化方法是一种非梯度优化算法。在整个过程中依次循环使用不同的坐标方向进行迭代，一个周期的一维搜索迭代过程相当于一个梯度下降的迭代。
3.梯度下降是利用目标函数的导数来确定搜索方向的，该梯度方向可能不与任何坐标轴平行。而坐标轴下降法法是利用当前坐标方向进行搜索，不需要求目标函数的导数，只按照某一坐标方向进行搜索最小值。。
4.两者都是迭代方法，且每一轮迭代，都需要O(mn)的计算量(m为样本数，n为系数向量的维度)。

线性回归的弹性网：

弹性网是一种使用L1和L2作为正则化矩阵的线性回归模型，是岭回归和坐标轴下降回归的组合。弹性网结合了岭回归的计算精准的优点，同时又结合了坐标轴下降回归特征选择的优势。

三种方法的选择：

岭回归优缺点：计算精准（没有丢失特征），不具备特征选择功能。
Lasso回归优缺点：具有特征选择功能，但是将某些θ化为0，可能导致特征丢失，使得最终的模型的偏差比较大在进行正则化的过程中，通常要先使用岭回归。
当特征非常多时，尤其当多个特征和另一个特征相关的时弹性网络非常有用。

sklearn中的线性回归：

标准线性回归，
带交叉验证的线性回归，
利用网格搜索进行高效参数调优。

#标准线性回归
from sklearn import linear_model
reg = linear_model.LinearRegression()
reg.fit([[0, 0], [1, 1], [2, 2]], [0, 1, 2])#训练，拟合
reg.predict([[4,4],[3,4]])#预测#岭回归
from sklearn import linear_model
reg = linear_model.Ridge(alpha=.5)
reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1])
reg.predict([[4,4],[3,4]])#预测#带交叉验证的岭回归
from sklearn import linear_model
reg = linear_model.RidgeCV(alphas=(0.1,0.2,0.3,0.4,0.5),cv=3)
reg.fit([[0, 0], [0, 0], [1, 1]], [0.1, 1])
reg.alpha_#坐标轴下降回归
from sklearn import linear_model
reg = linear_model.Lasso(alpha=0.1)
reg.fit([[0, 0], [1, 1]], [0, 1])
reg.predict([[1, 1]])#带交叉验证的lasso回归
from sklearn.linear_model import LassoCV
from sklearn.datasets import make_regression
X, y = make_regression(noise=4)
reg = LassoCV(alphas=(0.1,0.2,0.3,0.4,0.5),cv=5).fit(X, y)
reg.alpha_#弹性网
from sklearn.linear_model import ElasticNet
from sklearn.datasets import make_regression
X, y = make_regression(n_features=2)
regr = ElasticNet()
regr.fit(X, y)#带交叉验证的弹性网
from sklearn.linear_model import ElasticNetCV
from sklearn.datasets import make_regression
X, y = make_regression(noise=4)
regr = ElasticNetCV(l1_ratio=[.1, .5, .7, .9, .95, .99, 1],alphas=[1,2],cv=5).fit(X, y)
regr.alpha_#模型评价
from sklearn import metrics
# 用scikit-learn计算MSE
print("MSE:",metrics.mean_squared_error(y_test, y_pred))
# 用scikit-learn计算RMSE
print("RMSE:",np.sqrt(metrics.mean_squared_error(y_test, y_pred)))

网格搜索
from sklearn import linear_model
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV
def load_data():''' 加载数据'''data=pd.read_table('abalone.txt',sep=' ',names=['x1','x2','x3','x4','x5','x6','x7','x8','y'])feature=data.loc[:,'x1':'x8']label=data['y']return np.mat(feature), np.mat(label).T
X, y = load_data()
X, y = StandardScaler().fit_transform(X), StandardScaler().fit_transform(y)
params={'alpha':[0.001,0.01,0.1,1,10,100]}#粗调节
ridge = linear_model.Ridge()#
grid_search=GridSearchCV(ridge,params,cv=10)
grid_search.fit(X, y)
grid_search.best_params_

模型评价MSE或者RMSE，如果我们用其他方法得到了不同的系数：
#均方差（Mean Squared Error, MSE）或者均方根差(Root Mean Squared Error, RMSE)。