多项式线性回归实现波斯顿房价预测任务

1）采用纯python语言实现不能调用pytorch库函数

2）绘制损失函数变化曲线

3）分析模型超参数对算法性能的影响。

以下为初步学习，没有增加偏置以及划分数据集，单纯理解学习的过程以及特征多项式特征矩阵形成的过程，还存在一些数值溢出问题。但是能看出损失值在不断学习中变小了

# 导入必要的库
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
import itertools  # 用于生成特征组合

# 定义标准化函数，用于对数据进行标准化处理
# 波士顿房价数据集包含13个特征，例如犯罪率、住宅平均房间数等，
# 这些特征的数值范围可能非常不同。比如犯罪率可能在0到1之间变化，而住宅平均房间数则可能在几个到十几个之间变化。
def standardize(X):
    mean = X.mean(axis=0)  # 计算每列特征的平均值
    std = X.std(axis=0)  # 计算每列特征的标准差
    # 对数据进行标准化：(原始值 - 平均值) / 标准差
    return (X - mean) / std, mean, std
# 标准化之后的数据将具有相同的尺度，通常情况下，它们会被缩放至均值为0，标准差为1的标准正态分布形式。
# 例如，犯罪率特征的原始值可能从0.00632到0.87100不等，标准化后它的值将在大约-1到+1之间。
# 同样地，住宅平均房间数特征的原始值可能从3到8不等，标准化后也将位于-1到+1之间

# 定义多项式特征生成函数
def polynomial_features(X, degree):
    n_samples, n_features = X.shape  # 获取样本数量和特征数量
    # 定义一个内部函数，用于生成多项式特征的索引组合,会返回所有 1到degree阶数 的特征组合
    # 如特征有012，那么若degree=2，会返回0，1，2，01，02，12，00，11，22 这个就是输出的特征数量
    def index_combinations():
        for i in range(1, degree + 1):  # 遍历多项式的阶数
            for combo in itertools.combinations_with_replacement(range(n_features), i):  # 生成特征的组合
                yield combo  # 返回一个组合
    # 计算输出特征的数量
    n_output_features = len(list(index_combinations()))
    X_poly = np.empty((n_samples, n_output_features))  # 初始化多项式特征矩阵
    # 构建多项式特征矩阵
    for i, combo in enumerate(index_combinations()):  # 遍历所有组合
        X_poly[:, i] = np.prod(X[:, combo], axis=1)  # 计算组合特征的乘积
    return X_poly  # 返回多项式特征矩阵
# 上述函数：大概为比如原始特征有犯罪率，房间个数等等，规范化数据之后根据特征式阶数得出一堆输出特征：
# （犯罪率×房间数，房间数²等等）然后 n_output_features是输出特征的数量，作为初始化多项式特征矩阵的列数
# 行数就是n_samples，也就是样本数量（脑补一下这个矩阵）通过计算得出该矩阵列数为104，也就是新特征的数量

# # 定义均方误差计算函数
# def mean_squared_error(y_true, y_pred):
#     return np.mean((y_true - y_pred) ** 2)  # 均方误差公式
# # 定义均方误差计算函数
def mean_squared_error(y_true, y_pred):
    diff = y_true - y_pred
    # 使用 np.nanmean 以防有 NaN 值
    return np.nanmean(diff ** 2)

# 定义梯度下降优化算法
def gradient_descent(X, y, learning_rate, epochs, poly_degree):
    m, n = X.shape  # 获取样本数量和特征数量
    theta = np.random.randn(n)  # 初始化权重向量
    costs = []  # 存储每次迭代的损失值
    for _ in range(epochs):  # 迭代次数
        predictions = np.dot(X, theta)  # 预测值，用新特征矩阵和每个新特征的权重点乘得出预测该样本（楼房）的价格
        error = predictions - y  # 预测误差，预测价格和实际价格（标签）相比较
        gradient = 2 * np.dot(X.T, error) / m  # 计算梯度
        theta -= learning_rate * gradient  # 更新权重
        cost = mean_squared_error(y, predictions)  # 计算损失值
        costs.append(cost)  # 存储损失值
    return theta, costs  # 返回最终权重和所有损失值

# 加载波士顿房价数据集
boston = load_boston()
X = boston.data  # 特征数据
y = boston.target  # 目标数据（房价）

# 标准化特征
X, _, _ = standardize(X)

# 定义超参数
learning_rate = 0.01  # 学习率
epochs = 1000  # 迭代次数
poly_degree = 2  # 多项式的阶数

# 生成多项式特征矩阵，这个矩阵的列数有104，也就是新特征的数量是104个
X_poly = polynomial_features(X, poly_degree)

# 训练模型，将生成的新的特征式矩阵输入进去，会返回出每个新特征值的权重
theta, costs = gradient_descent(X_poly, y, learning_rate, epochs, poly_degree)


# 绘制损失曲线
plt.figure(figsize=(10, 5))
plt.plot(range(1, epochs+1), costs)
plt.xlabel("Epoch")  # x轴标签
plt.ylabel("Cost")  # y轴标签
plt.title("Cost Function Over Time")  # 图表标题
plt.show()  # 显示图表

# 打印最终系数
print("Coefficients:", theta)

# 分析不同超参数的影响
degrees = [1, 2, 3]  # 多项式阶数
learning_rates = [0.001, 0.01, 0.1]  # 学习率
# 对不同的多项式阶数和学习率进行测试
for degree in degrees:
    X_poly = polynomial_features(X, degree)
    for lr in learning_rates:
        _, costs = gradient_descent(X_poly, y, lr, epochs, degree)
        print(f"Degree={degree}, Learning Rate={lr}: Final Cost={costs[-1]}")