【MATLAB技巧】【免费的小例程】线性拟合的多种方法：polyfit、矩阵、带权重、稳健回归

最新推荐文章于 2025-07-30 20:07:32 发布

MATLAB卡尔曼

最新推荐文章于 2025-07-30 20:07:32 发布

阅读量571

点赞数 3

CC 4.0 BY-SA版权

分类专栏： MATLAB技巧免费的小例程文章标签： matlab 矩阵回归线性拟合

本文链接：https://blog.csdn.net/callmeup/article/details/149666051

MATLAB技巧同时被 2 个专栏收录

64 篇文章

订阅专栏

免费的小例程

10 篇文章

订阅专栏

在这里插入图片描述

本文给出MATLAB进行线性拟合的详细示例及多种方法说明，附例程代码和运行结果

基础线性拟合（polyfit函数）

适用场景：单变量线性回归
代码实现：

% 生成示例数据
x = [0, 0.5, 1, 2, 4, 6]';
y = [0, 0.21723, 0.43445, 0.99924, 2.43292, 4.77895]';

% 一阶多项式拟合（线性）
P = polyfit(x, y, 1);% P返回系数 [β1, β0]
y_fit = polyval(P, x);

% 绘图
plot(x, y, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(x, y_fit, 'r--', 'LineWidth', 2);
xlabel('x'); ylabel('y');
legend('原始数据', '拟合直线');

输出说明：

P(1)为斜率（β1），P(2)为截距（β0）。
结果误差可通过残差平方和评估：RSS = sum((y - y_fit).^2)。

在这里插入图片描述

矩阵法求解正规方程

适用场景：多元线性回归（含截距项）
数学原理：
通过最小二乘法求解参数向量：
$\hat{\beta} = (X^T X)^{-1} X^T Y$
其中设计矩阵 $X$ 需包含全1列以表示截距。

代码实现：

clc; clear; close all;

% 数据矩阵构造（参考正规方程理论）
x = [0, 0.5, 1, 2, 4, 6]';
y = [0, 0.21723, 0.43445, 0.99924, 2.43292, 4.77895]';
% 正规方程求解（矩阵法）
beta_hat = (x'*x) \ (x'*y); % 核心公式：β = (X^T X)^{-1} X^T y

% 残差分析与误差计算
y_fit = x * beta_hat;
% 绘图
plot(x, y, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(x, y_fit, 'r--', 'LineWidth', 2);
xlabel('x'); ylabel('y');
legend('原始数据', '拟合直线');

优势：适用于多变量场景，可手动控制误差分析（如协方差矩阵计算）。
运行结果：
在这里插入图片描述

带权重的线性拟合

适用场景：数据存在异方差性（误差非均匀）
代码扩展：

clc; clear; close all;
x = [0, 0.5, 1, 2, 4, 6]';
y = [0, 0.21723, 0.43445, 0.99924, 2.43292, 4.77895]';
y_error = [1,1,1,1,10,1]';
weights = 1 ./ y_error.^2;% 假设y_error为测量误差向量
% 加权最小二乘求解
W = diag(weights);  % 权重矩阵
beta_hat = (x' * W * x) \ (x' * W * y);
y_fit = x * beta_hat;

% 绘图
plot(x, y, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(x, y_fit, 'r--', 'LineWidth', 2);
xlabel('x'); ylabel('y');
legend('原始数据', '拟合直线');

说明：

加权最小二乘法通过调整weights参数实现，常用于仪器精度不同的数据。
运行结果：

在这里插入图片描述

稳健回归（LIBRA工具箱）

适用场景：数据含异常值
实现步骤：
下列代码含内嵌的函数，不能直接在命令行窗口运行，需要先粘贴到空脚本中，然后运行脚本

clc; clear; close all;
x = [0, 0.5, 1, 2, 4, 6]';
y = [0, 0.21723, 0.43445, 0.99924, 2.43292, 4.77895]';
beta_huber = huber_wls(x, y, 20);
y_fit = x * beta_huber;

% 绘图
plot(x, y, 'ko', 'MarkerFaceColor', 'k'); hold on;
plot(x, y_fit, 'r--', 'LineWidth', 2);
xlabel('x'); ylabel('y');
legend('原始数据', '拟合直线');

function beta = huber_wls(X, y, max_iter)
    % 初始化参数
    beta = X \ y;  % OLS初步估计
    c = 1.345;     % Huber参数
    
    for iter = 1:max_iter
        residuals = y - X*beta;
        sigma = 1.4826 * median(abs(residuals - median(residuals)));  % MAD估计
        r_std = residuals / sigma;
        
        % 计算Huber权重
        w = ones(size(r_std));
        outlier_idx = abs(r_std) > c;
        w(outlier_idx) = c ./ abs(r_std(outlier_idx));
        
        % 加权最小二乘更新
        W = diag(w);
        beta = (X'*W*X) \ (X'*W*y);
    end
end