线性回归(梯度下降)

首先说案例：

房子的价格和所占面积有着很大的关系，假如现在有一些关于房子面积和价格的数据，我要如何根据已经有的数据来判断未知的数据呢？

假如x(房屋面积)，y(房屋价格)

x=[ 56 72 69 88 102 86 76 79 94 74]

y=[92, 102, 86, 110, 130, 99, 96, 102, 105, 92]

单变量梯度下降

我们可以看到x和y有一定的关系。假如，我们令

$y=\omega x + b$

这里有两个未知数( $\omega$ ,b)，我们要如何求得这两个未知数?

我们现在知道的是什么？是所有的x、y的真实值。我们想干什么？想要用公式来预测。预测什么？根据x预测y。我们可以不断尝试改变w和b的值，来判断其预测的准确性。因此这里的w和b是变量。我们要如何衡量这个变量取值后预测的准确性呢？我们可以看下面的图：红线是预测直线，点是真实的面积和价格。可以看到其中不乏存在误差，我们的目标就是使误差达到最小值。那么我们如何来衡量误差呢？我们加下来引进以下术语：

模型函数： $y= \omega x + b$ ，也可写成： $h_{\omega,b}(x)=\omega x + b$

规定：

$x^{(i)}$	第i条记录(样本)
m	总样本数

损失函数：

我们要定义一个函数来衡量误差的大小，这个函数称为损失函数

那么如何定义损失函数才能够准确衡量呢？

首先要体现预测值与真实值的误差，那就需要有预测值与真实值的差。

$J(\omega,b)=\frac{\sum_{i=1}^{m}(h_{w,b}(x^{(i)})-y^{(i)})^{2}}{2m}$

这里用到平方使值都为正数，以便不会相加的时候正负抵消，这里除以m是为了使随着m的增大让损失函数不会递增，这里分母上的2是为了后续计算简便。

我们现在的任务就是求解损失函数J的最小值时对应的w和b。

这里我们引入梯度下降法。网上有很多讲解，这里就不展开。

这里展示对J求导的具体过程：

$\frac{\partial J(\omega,b)}{\partial \omega} =\frac{1}{2m}\frac{\partial (\sum_{i=1}^{m}(\omega x^{(i)}+b-y^{(i)})^{2})}{\partial \omega} =\frac{\sum_{i=1}^{m}2 (\omega x^{(i)}+b-y^{(i)})x^{(i)}}{2m} =\frac{\sum_{i=1}^{m}(h_{\omega,b}(x^{(i)})-y^{(i)})x^{(i)}}{m}$

梯度下降：

$\omega=w-\alpha \frac{\partial J(\omega,b)}{\partial w}$

$b=b-\alpha \frac{\partial J(\omega,b)}{\partial b}$

多变量梯度下降

下面看向量化的内容：

上面的参数只有一个，但是现实生活中影响因子可不止一个，假如有10000个我们肯定不能像上面那样写成一个乘一个。

波士顿房价数据集是机器学习中非常知名的数据集，它被用于多篇回归算法研究的学术论文中。该数据集共计 506 条，其中包含有 13 个与房价相关的特征以及 1 个目标值（房价）。

数据在这里：https://labfile.oss.aliyuncs.com/courses/1081/course-5-boston.csv

该数据集统计了波士顿地区各城镇的住房价格中位数，以及与之相关的特征。每列数据的列名解释如下：

CRIM: 城镇犯罪率。
ZN: 占地面积超过 2.5 万平方英尺的住宅用地比例。
INDUS: 城镇非零售业务地区的比例。
CHAS: 查尔斯河是否经过 (=1 经过，=0 不经过)。
NOX: 一氧化氮浓度（每 1000 万份）。
RM: 住宅平均房间数。
AGE: 所有者年龄。
DIS: 与就业中心的距离。
RAD: 公路可达性指数。
TAX: 物业税率。
PTRATIO: 城镇师生比例。
BLACK: 城镇的黑人指数。
LSTAT: 人口中地位较低人群的百分数。
MEDV: 城镇住房价格中位数。（需要预测的）

我们现在有n个参数 $\omega_{1}\omega_{2}\omega_{3} \cdots \omega_{n-1}\omega_{n}$ ,写成向量的形式就是 $\vec{\omega}=[\omega_{1}\cdots \omega_{n}]$

现在模型就是： $f_{\vec{\omega},b}=\vec{\omega}\cdot \vec{x}+b$

损失函数是： $J(\vec{\omega},b)$

梯度下降重复操作是：

repeat:

$\omega_{j}=w_{j}-\alpha \frac{\partial J(\vec{\omega},b)}{\partial \omega_{j}}$

$b=b-\alpha \frac{\partial J(\vec{\omega},b)}{\partial b}$

接下来我们具体对损失函数求导：

接下来介绍一些线性回归的一些API：

import sklearn.linear_model as lm
#创建模型
model = lm.LinearRegression()
#训练模型
#输入为一个二维数组表示的样本矩阵
#输出为每个样本最终结果
model.fit(输入,输出)#
#预测输出
#输入array是一个二维数组，每行是一个样本，每一列是一个特性
result = model.perdict(array)

假如下面是总体数据：

x1	x2	x3	y

那么这个就是输入：

这就是输出：

这里使用上面的数据进行测试：

我们的思路如下：

import pandas as pd
import numpy as np
import matplotlib

matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression

features = pd.read_csv("https://labfile.oss.aliyuncs.com/courses/1081/course-5-boston.csv")
target = features['medv']
split_num = int(features.shape[0] * 0.7)  # 得到 70% 位置  # 目标值数据
X_train = features[:split_num]  # 训练集特征
y_train = target[:split_num]  # 训练集目标

X_test = features[split_num:]  # 测试集特征
y_test = target[split_num:]  # 测试集目标

model = LinearRegression()  # 建立模型
model.fit(X_train, y_train)
print("系数为", model.coef_)
print("截距为", model.intercept_)

pre = model.predict(X_test)


# 计算平均绝对误差（MAE）：
def mae_value(y_true, y_pred):
    n = y_true.shape[0]
    mae = sum(np.abs(y_true - y_pred)) / n
    return mae


# 计算均方误差（MSE）
def mse_value(y_true, y_pred):
    n = y_true.shape[0]
    mse = sum(np.square(y_true - y_pred)) / n
    return mse


mae = mae_value(y_test.values, pre)
mse = mse_value(y_test.values, pre)

print("MAE: ", mae)
print("MSE: ", mse)
plt.figure(figsize=(10, 6))
plt.scatter(y_test.index,y_test.values, color='blue', label='真实值', marker='o')
plt.scatter(y_test.index,pre, color='red', label='预测值', marker='x')
plt.title("预测值与真实值比较")
plt.xlabel('样本编号')
plt.ylabel('房价')
plt.legend(loc='best')
plt.grid(True)
plt.show()

结果：
系数为 [-3.89487871e-15 -5.03069808e-15 3.26128013e-16 9.53456200e-15
1.23485884e-14 -1.95123865e-15 -4.26741975e-15 6.41414005e-16
2.91894330e-16 2.81025203e-16 1.00033913e-15 -5.89805982e-17
-1.12323345e-16 1.00000000e+00]
截距为 2.4513724383723456e-13
MAE: 6.197381789039163e-14
MSE: 7.683577875093926e-27

因为数据给的有点刻意，所以预测结果非常匹配。

注意：这里使用的学习率和最大迭代次数都是固定的，如果想要修改可以通过如下：

model = SGDRegressor(learning_rate='constant', eta0=0.01, max_iter=1000, tol=1e-3)

其中：

learning_rate：
- 这个参数决定了学习率的调度方法（schedule）。它可以取以下值之一：
  - 'constant'：使用恒定的学习率，即 eta0。
  - 'optimal'：使用一个自适应的学习率，这个学习率取决于初始学习率、当前步数和正则化项。
  - 'invscaling'：使用逆比例调度的学习率，学习率会随时间步数的增加而减少。
  - 'adaptive'：如果一次训练中的损失未减少，学习率会逐步减少。
eta0：
- 初始学习率。当 learning_rate 设为 'constant'、'invscaling' 或 'adaptive' 时，这个值就决定了模型的初始学习率。对于 SGDRegressor，eta0 是模型学习参数更新的步长。
max_iter：
- 最大迭代次数。它决定了训练模型时允许的最大迭代次数。对于大型数据集，可能需要更多的迭代次数才能使模型收敛。
tol：
- 容差（tolerance）。训练过程会在损失函数的变化小于这个容差值时提前停止，意味着模型已基本收敛。如果设置得过小，模型可能会运行更多的迭代次数才能收敛。

下面看一个纯手敲的版本：
添加一个截距项，方便矩阵相乘得出，偏移项b

# 增加一个截距项
X_train = np.c_[np.ones(X_train.shape[0]), X_train]
X_test = np.c_[np.ones(X_test.shape[0]), X_test]

这个代码：X_train.shape得到的是一个元组(行数，列数),再[0]取出行数。然后用np.ones创建一个这么多行的全为1的矩阵。np.c_是将两个矩阵连接（水平连接）

接下来是梯度下降的过程：

def gradient_descent(X, y, learning_rate=0.01, max_iter=1000, tol=1e-3):
    m, n = X.shape
    theta = np.zeros(n)
    for iteration in range(max_iter):
        gradient = X.T.dot(X.dot(theta) - y) / m
        new_theta = theta - learning_rate * gradient
        if np.linalg.norm(new_theta - theta, ord=1) < tol:
            break
        theta = new_theta
    return theta

输入参数:
- X: 特征矩阵，其中每行代表一个样本，每列代表一个特征。
- y: 目标变量的向量。
- learning_rate: 学习率，控制每次更新的步长，默认为0.01。
- max_iter: 最大迭代次数，用于控制梯度下降的迭代次数，默认为1000。
- tol: 容忍度，用于控制算法收敛的阈值，默认为1e-3。
初始化参数:
- m, n = X.shape: 获取特征矩阵 X 的行数（样本数）和列数（特征数）。
- theta = np.zeros(n): 初始化参数向量 theta，全部设为零，长度与特征数相同。
梯度下降迭代:
- for iteration in range(max_iter):：循环执行梯度下降迭代。
- gradient = X.T.dot(X.dot(theta) - y) / m: 计算梯度，这里使用了线性回归模型的损失函数的梯度公式。
- new_theta = theta - learning_rate * gradient: 更新参数向量 theta，向梯度的负方向移动一定步长。
- if np.linalg.norm(new_theta - theta, ord=1) < tol:：检查参数更新的变化是否小于阈值 tol，如果满足要求则停止迭代。
- theta = new_theta: 将更新后的参数向量赋值给 theta 继续下一次迭代。
返回结果:
- 返回最终收敛的参数向量 theta，该向量能够使得线性回归模型拟合样本数据最佳。