一、一元线性回归

1.一元线性回归的梯度下降法

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  4. # 载入数据
  5. data = np.genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\data.csv", delimiter=",") #注意这是csv文件的读取，第二个参数是指定分隔符符号
  6. x_data = data[:,0] #选取csv的第一列数据
  7. y_data = data[:,1] #选取csv的第二列数据
  8. plt.scatter(x_data,y_data) #散点图显示，做回归之前最好还是先看看数据点的分布情况
  9. plt.show()
  11. # 学习率learning rate
  12. lr = 0.0001 #定义学习率，开始时设定较小一下
  13. # 截距
  14. b = 1 #指定初值
  15. # 斜率
  16. k = 1 #指定初值
  17. # 最大迭代次数
  18. epochs = 50 #指定迭代次数，一般梯度下降法都要设置一个迭代的次数保证算法的收敛性
  20. # 最小二乘法的应用
  21. def compute_error(b, k, x_data, y_data): #这里是构造损失函数
  22.     totalError = 0 #先设置代价函数的初值为0
  23.     for i in range(0, len(x_data)):
  24.         totalError += (y_data[i] - (k * x_data[i] + b)) ** 2 #相当于sum求和求残差平方和，y_data是真实值，k * x_data[i] + b是预测值
  25.     return totalError / float(len(x_data)) / 2.0 #这里返回的是之前代价函数的公式情况
  27. #梯度下降的具体应用
  28. def gradient_descent_runner(x_data, y_data, b, k, lr, epochs):
  29.     # 计算总数据量
  30.     m = float(len(x_data))
  31.     # 循环epochs次
  32.     for i in range(epochs):
  33.         b_grad = 0 #定义初值
  34.         k_grad = 0 #定义初值
  35.         # 计算梯度的总和再求平均
  36.         for j in range(0, len(x_data)):
  37.             b_grad += (1/m) * (((k * x_data[j]) + b) - y_data[j]) #这里的+=还是一个最小二乘中求和的形式
  38.             k_grad += (1/m) * x_data[j] * (((k * x_data[j]) + b) - y_data[j])
  39.         # 更新b和k
  40.         b = b - (lr * b_grad)
  41.         k = k - (lr * k_grad)
  42.         # 每迭代5次，输出一次图像
  43.         if i % 5==0:
  44.             print("epochs:",i)
  45.             plt.plot(x_data, y_data, 'b.')
  46.             plt.plot(x_data, k*x_data + b, 'r')
  47.             plt.show()
  48.     return b, k
  50. print("Starting b = {0}, k = {1}, error = {2}".format(b, k, compute_error(b, k, x_data, y_data))) #依次打印b的值，k的值，代价函数的值，大括号中的0，1，2表示传入数据的位次
  51. print("Running...")
  52. b, k = gradient_descent_runner(x_data, y_data, b, k, lr, epochs)
  53. print("After {0} iterations b = {1}, k = {2}, error = {3}".format(epochs, b, k, compute_error(b, k, x_data, y_data))) #输出迭代次数到50次时的值
  56. # 画图→迭代第50次的值
  57. plt.plot(x_data, y_data, 'b.')
  58. plt.plot(x_data, k*x_data + b, 'r')
  59. plt.show()

2.一元线性回归的sklearn库的应用

  1. from sklearn.linear_model import LinearRegression
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  5. # 载入数据
  6. data = np.genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\data.csv", delimiter=",")
  7. x_data = data[:,0]
  8. y_data = data[:,1]
  9. plt.scatter(x_data,y_data)
  10. plt.show()
  11. print(x_data.shape) #这里显示x_data的形状，会发现这是一个向量形式，后面要对它转变形式，sklearn库才会识别并拟合模型
  13. #x_data = data[:,0,np.newaxis] #这个代码表示100行1列的数据，这里相当于给数据加上一个维度，否则sklearn库模型拟合会出现错误
  14. #x_data.shape #重新检查x_data的形状
  17. x_data = data[:,0,np.newaxis] 
  18. y_data = data[:,1,np.newaxis]
  19. # 创建并拟合模型
  20. model = LinearRegression() #先创建一个模型
  21. model.fit(x_data, y_data) #进行建模
  22. display(model.intercept_)  #截距
  23. display(model.coef_)  #线性模型的系数
  24. model.score(x_data,y_data) #计算R方
  26. # 画图
  27. plt.plot(x_data, y_data, 'b.')
  28. plt.plot(x_data, model.predict(x_data), 'r')
  29. plt.show()

二、多元线性回归

1.多元线性回归的梯度下降法

  1. import numpy as np
  2. from numpy import genfromtxt
  3. import matplotlib.pyplot as plt  
  4. from mpl_toolkits.mplot3d import Axes3D  #这个库是画3D图的
  6. # 读入数据 
  7. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\Delivery.csv",delimiter=',')
  8. print(data)
  10. # 切分数据
  11. x_data = data[:,:-1] #取前两列数据，作为x向量的值,属性是numpy的一个二维数组
  12. y_data = data[:,-1] #取y值
  13. print(x_data)
  14. print(y_data)
  16. # 学习率learning rate
  17. lr = 0.0001
  18. # 参数
  19. theta0 = 0 
  20. theta1 = 0
  21. theta2 = 0 #设置三个参数
  22. # 最大迭代次数
  23. epochs = 1000
  25. # 最小二乘法
  26. def compute_error(theta0, theta1, theta2, x_data, y_data):
  27.     totalError = 0
  28.     for i in range(0, len(x_data)):
  29.         totalError += (y_data[i] - (theta1 * x_data[i,0] + theta2 * x_data[i,1] + theta0)) ** 2 #最小二乘法，这里相当于求和
  30.     return totalError / float(len(x_data))
  32. def gradient_descent_runner(x_data, y_data, theta0, theta1, theta2, lr, epochs):
  33.     # 计算总数据量
  34.     m = float(len(x_data))
  35.     # 循环epochs次
  36.     for i in range(epochs):
  37.         theta0_grad = 0
  38.         theta1_grad = 0
  39.         theta2_grad = 0 #先进行初始化
  40.         # 计算梯度的总和再求平均
  41.         for j in range(0, len(x_data)):
  42.             theta0_grad += (1/m) * ((theta1 * x_data[j,0] + theta2*x_data[j,1] + theta0) - y_data[j]) #利用for循环进行求和过程！预测值减真实值！
  43.             theta1_grad += (1/m) * x_data[j,0] * ((theta1 * x_data[j,0] + theta2*x_data[j,1] + theta0) - y_data[j])
  44.             theta2_grad += (1/m) * x_data[j,1] * ((theta1 * x_data[j,0] + theta2*x_data[j,1] + theta0) - y_data[j])
  45.         # 更新b和k
  46.         theta0 = theta0 - (lr*theta0_grad)
  47.         theta1 = theta1 - (lr*theta1_grad)
  48.         theta2 = theta2 - (lr*theta2_grad)
  49.     return theta0, theta1, theta2
  51. print("Starting theta0 = {0}, theta1 = {1}, theta2 = {2}, error = {3}".
  52.       format(theta0, theta1, theta2, compute_error(theta0, theta1, theta2, x_data, y_data)))
  53. print("Running...")
  54. theta0, theta1, theta2 = gradient_descent_runner(x_data, y_data, theta0, theta1, theta2, lr, epochs)
  55. print("After {0} iterations theta0 = {1}, theta1 = {2}, theta2 = {3}, error = {4}".
  56.       format(epochs, theta0, theta1, theta2, compute_error(theta0, theta1, theta2, x_data, y_data)))
  58. #一下代码为画图
  59. ax = plt.figure().add_subplot(111, projection = '3d') 
  60. ax.scatter(x_data[:,0], x_data[:,1], y_data, c = 'r', marker = 'o', s = 100) #点为红色三角形  
  61. x0 = x_data[:,0]
  62. x1 = x_data[:,1]
  63. # 生成网格矩阵
  64. x0, x1 = np.meshgrid(x0, x1) #注意！这是生成一个较为密集的网格矩阵，便于观察点的位置
  65. z = theta0 + x0*theta1 + x1*theta2
  66. # 画3D图
  67. ax.plot_surface(x0, x1, z)
  68. #设置坐标轴  
  69. ax.set_xlabel('Miles')  
  70. ax.set_ylabel('Num of Deliveries')  
  71. ax.set_zlabel('Time')  
  73. #显示图像  
  74. plt.show()

注：网格矩阵np.meshgrid()的用法

  1. x0,x1 = np.meshgrid([1,2,3],[4,5,6])
  2. x0 #代表横坐标
  3. array([[1, 2, 3],
  4.        [1, 2, 3],
  5.        [1, 2, 3]]) #x0的output
  6. x1 #代表纵坐标
  7. array([[4, 4, 4],
  8.        [5, 5, 5],
  9.        [6, 6, 6]]) #x1的output
  10. #上述的x0，x1可以表示成
  11. (1,4)(2,4)(3,4)
  12. (1,5)(2,5)(3,5)
  13. (1,6)(2,6)(3,6)
  15. plt.scatter(x0,x1)
  16. plt.show() #生成一个网格矩阵的模型样式

image.png 此为图的output情况

2.多元线性回归的sklearn库方法

（1）基本形式

  1. import numpy as np
  2. from numpy import genfromtxt
  3. from sklearn import linear_model
  4. import matplotlib.pyplot as plt  
  5. from mpl_toolkits.mplot3d import Axes3D  
  7. # 读入数据 
  8. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\Delivery.csv",delimiter=',')
  9. print(data)
  11. # 切分数据
  12. x_data = data[:,:-1]
  13. y_data = data[:,-1]
  14. print(x_data)
  15. print(y_data)
  17. # 创建模型
  18. model = linear_model.LinearRegression()
  19. model.fit(x_data, y_data) #利用fit方法来进行模型拟合
  21. # 系数
  22. print("coefficients:",model.coef_)
  23. # 截距
  24. print("intercept:",model.intercept_)
  25. # 测试
  26. x_test = [[102,4]]
  27. predict = model.predict(x_test)
  28. print("predict:",predict)
  30. ax = plt.figure().add_subplot(111, projection = '3d') 
  31. ax.scatter(x_data[:,0], x_data[:,1], y_data, c = 'r', marker = 'o', s = 100) #点为红色三角形  
  32. x0 = x_data[:,0]
  33. x1 = x_data[:,1]
  34. # 生成网格矩阵
  35. x0, x1 = np.meshgrid(x0, x1)
  36. z = model.intercept_ + x0*model.coef_[0] + x1*model.coef_[1]
  37. # 画3D图
  38. ax.plot_surface(x0, x1, z)
  39. #设置坐标轴  
  40. ax.set_xlabel('Miles')  
  41. ax.set_ylabel('Num of Deliveries')  
  42. ax.set_zlabel('Time')  
  44. #显示图像  
  45. plt.show()

（2）其他用法

详见https://blog.csdn.net/weixin_40014576/article/details/79918819查询具体用法

三、多项式回归

1.sklearn库的多项式回归应用

  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. from sklearn.preprocessing import PolynomialFeatures
  4. from sklearn.linear_model import LinearRegression
  6. # 载入数据
  7. data = np.genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\job.csv", delimiter=",")
  8. x_data = data[1:,1]
  9. y_data = data[1:,2]
  10. plt.scatter(x_data,y_data) #先观察散点情况，根据分布情况拟合曲线
  11. plt.show()
  13. x_data = x_data[:,np.newaxis] #原始的x_data形式是不合适的，sklearn库封装的函数使用这种形式会报错(需要二维格式),于是需要修改!
  14. y_data = y_data[:,np.newaxis]
  16. # 创建并拟合模型(一元线性回归形式)
  17. model = LinearRegression()
  18. model.fit(x_data, y_data)
  19. # 画图
  20. plt.plot(x_data, y_data, 'b.')
  21. plt.plot(x_data, model.predict(x_data), 'r')
  22. plt.show()
  23. model.score(x_data,y_data) #其实可以发现R方比较小，拟合效果比较差
  25. # 定义多项式回归,degree的值可以调节多项式的特征
  26. poly_reg  = PolynomialFeatures(degree=5) #degree表示拟合的最高次数
  27. # 特征处理
  28. x_poly = poly_reg.fit_transform(x_data) #此步最前面会产生1这个偏置值，这是用来拟合β0的
  29. # 定义回归模型
  30. lin_reg = LinearRegression()
  31. # 训练模型
  32. lin_reg.fit(x_poly, y_data)
  34. #注：x_poly的值如下表现
  35. array([[1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00, 1.0000e+00,
  36.         1.0000e+00],
  37.        [1.0000e+00, 2.0000e+00, 4.0000e+00, 8.0000e+00, 1.6000e+01,
  38.         3.2000e+01],
  39.        [1.0000e+00, 3.0000e+00, 9.0000e+00, 2.7000e+01, 8.1000e+01,
  40.         2.4300e+02],
  41.        [1.0000e+00, 4.0000e+00, 1.6000e+01, 6.4000e+01, 2.5600e+02,
  42.         1.0240e+03],
  43.        [1.0000e+00, 5.0000e+00, 2.5000e+01, 1.2500e+02, 6.2500e+02,
  44.         3.1250e+03],
  45.        [1.0000e+00, 6.0000e+00, 3.6000e+01, 2.1600e+02, 1.2960e+03,
  46.         7.7760e+03],
  47.        [1.0000e+00, 7.0000e+00, 4.9000e+01, 3.4300e+02, 2.4010e+03,
  48.         1.6807e+04],
  49.        [1.0000e+00, 8.0000e+00, 6.4000e+01, 5.1200e+02, 4.0960e+03,
  50.         3.2768e+04],
  51.        [1.0000e+00, 9.0000e+00, 8.1000e+01, 7.2900e+02, 6.5610e+03,
  52.         5.9049e+04],
  53.        [1.0000e+00, 1.0000e+01, 1.0000e+02, 1.0000e+03, 1.0000e+04,
  54.         1.0000e+05]])
  56. # 画图
  57. plt.plot(x_data, y_data, 'b.')
  58. plt.plot(x_data, lin_reg.predict(poly_reg.fit_transform(x_data)), c='r') #这里要传入经过特征处理的数据才对
  59. plt.title('Truth or Bluff (Polynomial Regression)')
  60. plt.xlabel('Position level')
  61. plt.ylabel('Salary')
  62. plt.show()
  63. lin_reg.score(poly_reg.fit_transform(x_data),y_data) #计算R方
  64. lin_reg.coef_ #计算回归系数
  65. lin_reg.intercept_ #计算截距

四、线性回归的标准方程法

标准方程法即利用代数方法将系数以矩阵的形式求出，由多元线性回归知识可以得到：
βhat = （XX）XY
注：矩阵不可逆情况：
1.多重共线性
2.特征数据太多，样本数小于特征数量

  1. #这里是一元线性回归的方式，多元也可以类似
  2. import numpy as np
  3. from numpy import genfromtxt
  4. import matplotlib.pyplot as plt  
  6. # 载入数据
  7. data = np.genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\data.csv", delimiter=",")
  8. x_data = data[:,0,np.newaxis]
  9. y_data = data[:,1,np.newaxis]
  10. plt.scatter(x_data,y_data)
  11. plt.show()
  13. print(np.mat(x_data).shape) #转化成矩阵形式
  14. print(np.mat(y_data).shape) #转化成矩阵形式
  15. # 给样本添加偏置项
  16. X_data = np.concatenate((np.ones((100,1)),x_data),axis=1) #np.ones((100,1))这里是生成100行1列的数据为1的矩阵，添加偏置项；np.concatenate是用来进行矩阵的合并的，axis是指定合并的方向，1是列方向
  17. print(X_data.shape) #查看新数据的格式
  19. #注：观察X_data的形式
  20. # print(X_data[:3])
  21. #[[ 1.         32.50234527]
  22.  [ 1.         53.42680403]
  23.  [ 1.         61.53035803]] output形式
  27. # 标准方程法求解回归参数
  28. def weights(xArr, yArr): #定义回归系数公式形式
  29.     xMat = np.mat(xArr) #转化矩阵形式
  30.     yMat = np.mat(yArr) #转化矩阵形式
  31.     xTx = xMat.T*xMat # 矩阵乘法
  32.     # 计算矩阵的值,如果值为0，说明该矩阵没有逆矩阵
  33.     if np.linalg.det(xTx) == 0.0: #如果行列式等于0，则没有逆
  34.         print("This matrix cannot do inverse")
  35.         return
  36.     # xTx.I为xTx的逆矩阵
  37.     ws = xTx.I*xMat.T*yMat
  38.     return ws
  39. ws = weights(X_data,y_data)
  40. print(ws)
  42. # 画图
  43. x_test = np.array([[20],[80]])
  44. y_test = ws[0] + x_test*ws[1]
  45. plt.plot(x_data, y_data, 'b.')
  46. plt.plot(x_test, y_test, 'r')
  47. plt.show()

五、其他库的回归模型用法（statsmodel库的用法）

这个库可以跟R语言一样直接输出其他统计量的信息，详细用法请见https://blog.csdn.net/chongminglun/article/details/104242342

六、岭回归

为了解决(XX)无法计算的情况（不可逆），于是引入岭回归概念，在标准方程法中，令
w = （XX+λI)Xy
其中，λ是岭系数，I为单位矩阵
岭回归是一个有偏估计
λ的选取：选择λ值，使得：各回归系数的岭估计基本稳定；残差平方和（代价函数的第一个部分）增大不太多

1.sklearn库的岭回归

  1. import numpy as np
  2. from numpy import genfromtxt
  3. from sklearn import linear_model
  4. import matplotlib.pyplot as plt
  6. #读入数据 
  7. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\longley.csv",delimiter=',') #getfromtext只能获取数值数据，不能获取字符串等数据
  8. print(data)
  10. # 切分数据
  11. x_data = data[1:,2:]
  12. y_data = data[1:,1]
  13. print(x_data)
  14. print(y_data)
  16. # 创建模型
  17. # 生成50个值
  18. alphas_to_test = np.linspace(0.001, 1) #从0.001到1取50个λ值进行岭回归实验
  19. # 创建模型，保存误差值
  20. model = linear_model.RidgeCV(alphas=alphas_to_test, store_cv_values=True) #CV代表交叉验证，利用交叉验证法训练模型
  21. model.fit(x_data, y_data)
  22. model.coef_ #岭回归系数
  24. # 岭系数
  25. print(model.alpha_) #选取最合适的λ值
  26. # loss值
  27. print(model.cv_values_.shape) #model.cv_values_会产生16行50列的数据，16行是指交叉验证法在16个样本中随机选取一个样本当测试集，其余的样本当作训练集进行拟合，50列是指50个λ值对应的loss值 
  29. # 画图
  30. # 岭系数跟loss值的关系
  31. plt.plot(alphas_to_test, model.cv_values_.mean(axis=0)) #axis = 0表示对行16个值求平均值（每个值对应一个模型的loss）
  32. # 选取的岭系数值的位置
  33. plt.plot(model.alpha_, min(model.cv_values_.mean(axis=0)),'ro')
  34. plt.show()
  36. model.predict(x_data[2,np.newaxis]) #查看第二行数据对应预测的y值，这里要传入2维数据，否则模型会报错

2.标准方程法的岭回归

标准方程法需要自己进行交叉验证，这里书写是比较麻烦的，所以调包可能会更好一点，在这里只写算法思想之类的代码，后续需要再补充（学习交叉验证法的手写命令！）。

  1. import numpy as np
  2. from numpy import genfromtxt
  3. import matplotlib.pyplot as plt  
  5. # 读入数据 
  6. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\longley.csv",delimiter=',')
  7. print(data)
  9. # 切分数据
  10. x_data = data[1:,2:]
  11. y_data = data[1:,1,np.newaxis]
  12. print(x_data)
  13. print(y_data)
  15. print(np.mat(x_data).shape) #将数据转化成矩阵形式，方便标准方程法
  16. print(np.mat(y_data).shape) #将数据转化成矩阵形式，方便标准方程法
  17. # 给样本添加偏置项
  18. X_data = np.concatenate((np.ones((16,1)),x_data),axis=1)
  19. print(X_data.shape)
  21. #可以查看新的数据集的内容
  22. #print(X_data[:3])
  24. # 岭回归标准方程法求解回归参数
  25. def weights(xArr, yArr, lam=0.2):
  26.     xMat = np.mat(xArr)
  27.     yMat = np.mat(yArr)
  28.     xTx = xMat.T*xMat # 矩阵乘法
  29.     rxTx = xTx + np.eye(xMat.shape[1])*lam #np.eye(xMat.shape[1]）传入单位矩阵，这里要传入和xTx相同维度的的单位矩阵，于是利用.shape函数
  30.     # 计算矩阵的值,如果值为0，说明该矩阵没有逆矩阵
  31.     if np.linalg.det(rxTx) == 0.0:
  32.         print("This matrix cannot do inverse")
  33.         return
  34.     # xTx.I为xTx的逆矩阵
  35.     ws = rxTx.I*xMat.T*yMat
  36.     return ws
  38. ws = weights(X_data,y_data)
  39. print(ws) #打印岭回归系数值
  41. # 计算预测值
  42. np.mat(X_data)*np.mat(ws)

七、lasso回归

由于损失函数和lasso系数范围的切点多在坐标轴上，因此采用lasso回归会出现系数稀疏的现象，即有很多系数的估计值都为零（在λ较小的情况下时），因此lasso回归有时比岭回归具有优势， 因为可以剔除模型中不显著的变量，使得简化模型

  1. import numpy as np
  2. from numpy import genfromtxt
  3. from sklearn import linear_model
  5. # 读入数据 
  6. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\longley.csv",delimiter=',')
  7. print(data)
  9. # 切分数据
  10. x_data = data[1:,2:]
  11. y_data = data[1:,1]
  12. print(x_data)
  13. print(y_data)
  15. # 创建模型
  16. model = linear_model.LassoCV() #这里是应用lasso的交叉验证法择λ值
  17. model.fit(x_data, y_data)
  18. # lasso系数
  19. print(model.alpha_)
  20. # 相关系数
  21. print(model.coef_)
  23. model.predict(x_data[-2,np.newaxis]) #传入倒数第二行进行预测y值，还是要传入二维数组形式才不会报错

八、弹性网（Elastic Net）

在弹性网中，正则项部分为线性回归（python实现） - 图2，效果会比lasso和岭回归要好一点

  1. import numpy as np
  2. from numpy import genfromtxt
  3. from sklearn import linear_model
  5. # 读入数据 
  6. data = genfromtxt(r"C:\Users\小新\Desktop\课程pdf\py\机器学习\线性回归及非线性回归\longley.csv",delimiter=',')
  7. print(data)
  9. # 切分数据
  10. x_data = data[1:,2:]
  11. y_data = data[1:,1]
  12. print(x_data)
  13. print(y_data)
  15. # 创建模型
  16. model = linear_model.ElasticNetCV() #利用弹性网交叉验证法求λ值
  17. model.fit(x_data, y_data)
  19. # 弹性网系数
  20. print(model.alpha_)
  21. # 相关系数
  22. print(model.coef_)
  24. model.predict(x_data[-2,np.newaxis]) #预测