作业2-线性回归的Matlab代码实现

一、前言

相关配置：Matlab 2020a（版本的影响应该不大，.m代码基本都能运行，个人感觉就是Simulink对版本的要求高一些）

二、任务描述

基于近两节课的理论推导，用代码实现线性回归，并对预测结果进行分析。同时体会学习率、迭代次数对系统的影响（收敛速度、代码运行速度等）。

三、代码实现

我们将任务假设为房价的预测，研究单变量的线性回归，即输入特征只有一个（房屋面积）。使用一次函数来生成相关样本（这里的一次函数为笔者随意设置，并不与现实的情况相符合）。

% ----------------------原始数据集---------------------- %
x = linspace(0, 1000, 200); % 生成横坐标（房屋面积）数据
y = 1.5 * x + 100; % 生成房价
data = [x; y]; % 组合成二维数据

这样生成的样本（以50个样本为例）连起来就是一条直线。然而现实情况是，样本数据或多或少会存在“噪声”，房价与房屋面积也不会呈现完美的一次函数关系。所以我们需要对样本增加“噪声”。代码示例如下：

noise_level = 30;  % 定义噪声参数noise = noise_level * randn(size(data)); % 生成与数据相同大小的高斯噪声
noisy_data = data + noise; % 将噪声添加到数据中noisy_x = noisy_data(1,:); % 添加噪声后的特征（房屋面积）
noisy_y = noisy_data(2,:); % 添加噪声后的结果（房价）

其中randn 函数用于生成符合标准正态分布（均值为 0，标准差为 1）的随机数。增加“噪声”后的样本示例如下图所示。我们可以发现，样本点出现了一定的变化，呈现一种“随机”的状态。

在进行梯度下降迭代之前，我们还需要将样本数据进行归一化。对自变量和因变量进行归一化或标准化，使其值更为均衡，有利于帮助模型更快收敛。（笔者在之前的测试中发现，若不进行归一化，输出的截距θ0会为0或者其他奇奇怪怪的问题。）

【特别注意】最终我们预测结果的那条线的相关参数需要进行反归一化才行！！！

std_x = (noisy_x - mean(noisy_x)) / std(noisy_x); % 标准化x
std_y = (noisy_y - mean(noisy_y)) / std(noisy_y); % 标准化y

接下来就是进行梯度下降的迭代了。代码主要分为以下三步

①、依据输入特征，获得预测值（预测的房价）

②、计算梯度，获得当前的θ0和θ1（实际情况是需要对θ0和θ1求偏导，代码中直接给出了计算结果，详见《机器学习-周志华》P54的式3.5和式3.6）

③、更新参数，获取新的θ0和tθ1（这里需要注意学习率α不宜设置的过大。过大的α可能会导致损失函数J无法收敛）

④、计算损失值

主循环代码如下所示：

% ------------------------------主循环代码-------------------------------
% GPT的例程是固定一个迭代次数，例如1000次，次数到了就退出循环，输出结果
% 笔者尝试固定一个损失值，当损失小于固定值时才退出循环，但运行效果不佳，
% 故放弃了这种方案
% ----------------------------------------------------------------------
for i = 1 : iteration_numy_pred = theta1 * std_x + theta0; % 预测值% 计算梯度d_theta1 = (1/m) * (theta1 * sum(std_x.^2) - sum((std_y - theta0) .* std_x));d_theta0 = theta0 - (1/m) * sum(std_y - theta1 * std_x);% 更新参数theta1 = theta1 - alpha * d_theta1;theta0 = theta0 - alpha * d_theta0;loss = (1/m) * sum((y_pred - std_y).^2); % 计算损失J_history(i) = loss;
end

主循环结束运行的方式笔者认为可以有两种方法。第一种是如上图所示的，通过设定迭代次数，当迭代次数到了后，退出循环。第二种是人为设定一个损失值，当计算的损失值小于设定的损失值时，退出循环。在测试中，第二种方法效果并不好，有时代码会长时间运行，无法跳出循环，遂放弃。

到这一步，我们就完成“模型的训练”啦，在画图之前，别忘了我们还需要进行反归一化。

% 反归一化后，才是真正的theta0和theta1
theta1 = theta1 * (std(noisy_y) / std(noisy_x)); % 恢复原始斜率
theta0 = mean(noisy_y) - theta1 * mean(noisy_x); % 恢复原始截距

四、输出结果与分析

下图为迭代500次的预测结果，其中绿色的回归线与数据点基本吻合，说明模型预测效果较好。

此外迭代500次的损失函数的收敛图如下图所示，我们可以发现，损失值随着迭代次数的增加逐渐减小，呈现一种收敛的状态。

下图是迭代1000次的损失函数的收敛图，我们可以发现，迭代次数超过500后，损失值几乎不再发生变化（本身“损失”就已经很小了），所以无脑的增加迭代次数并不一定是好的，适可而止，在大模型训练中也能节约计算时间。

以下Matlab“命令窗口”打印的数据结果，我们可以发现，最终输出的线性回归方程与前文设置的一次方程 $y=1.5x+100$ 相近似。

迭代次数：500次
loss损失：0.011024
最终线性回归方程：y = 1.49x + 104.89
theta1 = 1.49, theta0 = 104.89
>>

通过绘制 $J(\theta _{0},\theta _{1})$ 的图像，我们也可以发现，最终的损失函数J会有一个全局最小的点。(这只是一个简单的测试样例）

五、程序代码

clc;
clear;
close all;% ----------------------原始数据集---------------------- %
x = linspace(0, 1000, 200); % 生成横坐标（房屋面积）数据
y = 1.5 * x + 100; % 生成房价
data = [x; y]; % 组合成二维数据% -------------------------变量------------------------ %
iteration_num = 500; % 迭代次数
noise_level = 30;    % 定义噪声参数
theta1 = 0;          % 初始化斜率
theta0 = 0;          % 初始化截距
alpha = 0.01;        % 学习率
m = length(y);       % 样本数量
J_history = zeros(iteration_num, 1); % 记录每次迭代的损失值noise = noise_level * randn(size(data)); % 生成与数据相同大小的高斯噪声
noisy_data = data + noise; % 将噪声添加到数据中noisy_x = noisy_data(1,:); % 添加噪声后的特征（房屋面积）
noisy_y = noisy_data(2,:); % 添加噪声后的结果（房价）std_x = (noisy_x - mean(noisy_x)) / std(noisy_x); % 标准化x
std_y = (noisy_y - mean(noisy_y)) / std(noisy_y); % 标准化y% ------------------------------主循环代码-------------------------------
% GPT的例程是固定一个迭代次数，例如1000次，次数到了就退出循环，输出结果
% 笔者尝试固定一个损失值，当损失小于固定值时才退出循环，但运行效果不佳，
% 故放弃了这种方案
% ----------------------------------------------------------------------
for i = 1 : iteration_numy_pred = theta1 * std_x + theta0; % 预测值% 计算梯度d_theta1 = (1/m) * (theta1 * sum(std_x.^2) - sum((std_y - theta0) .* std_x));d_theta0 = theta0 - (1/m) * sum(std_y - theta1 * std_x);% 更新参数theta1 = theta1 - alpha * d_theta1;theta0 = theta0 - alpha * d_theta0;loss = (1/m) * sum((y_pred - std_y).^2); % 计算损失J_history(i) = loss;
end% ------------------------------绘制图像1------------------------------- %
subplot(2, 1, 1); % 创建两个子图
plot(noisy_data(1,:), noisy_data(2,:), '.b');
grid on; % 添加网格
xlim([0, 1000]); % 设置x轴范围
ylim([0, 1800]); % 设置y轴范围
hold on; % 保持图形，防止被后续图形覆盖% 反归一化后，才是真正的theta0和theta1
theta1 = theta1 * (std(noisy_y) / std(noisy_x)); % 恢复原始斜率
theta0 = mean(noisy_y) - theta1 * mean(noisy_x); % 恢复原始截距y_fit = theta1 * noisy_x + theta0;  % 计算回归线的 y 值
plot(noisy_x, y_fit, '-g', 'LineWidth', 2); % 绘制回归线title('运行结果');          % 标题
xlabel('房屋面积');         % 横坐标标签
ylabel('房屋价格');         % 纵坐标标签
legend('数据点', '回归线');  % 图例hold off;% ------------------------------绘制图像2------------------------------- %
subplot(2, 1, 2); % 创建两个子图
plot(1 : iteration_num, J_history, '-r', 'LineWidth', 2);
grid on; % 添加网格
title('损失函数收敛图'); % 标题
xlabel('迭代次数');      % 横坐标标签
ylabel('损失值');        % 纵坐标标签% ------------------------------输出结果------------------------------- %
fprintf('迭代次数：%d次\n', iteration_num); 
fprintf('loss损失：%.6f\n', loss);
fprintf('最终线性回归方程：y = %.2fx + %.2f\n', theta1, theta0);
fprintf('theta1 = %.2f, theta0 = %.2f\n', theta1, theta0);% -----------------------------END OF FILE---------------------------- %