最小二乘法matlab代码

最小二乘法是一种常见的数值计算方法，常用于拟合数据和解决参数估计问题。MATLAB中可用polyfit函数实现最小二乘拟合。

假设有n个数据点(x_i, y_i)， 拟合函数为y = f(x)，其中f(x)为n次多项式，即f(x) = a_n * x^n + a_{n-1} * x^(n-1) + ... ... + a_1 * x + a_0，需要找到a_0, a_1, ..., a_n，使得拟合函数y = f(x)与真实数据点的离散程度最小。

实现步骤如下：

1. 计算n个数据点的平均值，即 x_bar = (x_1 + x_2 + ... ... + x_n) / n 和 y_bar = (y_1 + y_2 + ... ... + y_n) / n。将数据点相对于(x_bar, y_bar)平移，变为(x_i - x_bar, y_i - y_bar)。

2. 计算n个数据点的最小二乘矩阵 A 和 b（参考PPT最小二乘法）：


其中A为n * (n+1)/2阶Vandermonde矩阵，b为n * 1的向量。

3. 解Ax = b，求出拟合多项式的系数x = [a_0, a_1, ..., a_n]。

4. 绘制原始数据点和拟合曲线。

完整的MATLAB代码如下：

% 最小二乘拟合多项式函数
% 输入：x, y - 数据点，n - 多项式次数
% 输出：p - 拟合多项式系数，f - 拟合曲线
function [p, f] = polyfit_least_square(x, y, n)
    % 计算n个数据点的平均值
    x_bar = mean(x);
    y_bar = mean(y);

    % 将数据点相对于(x_bar, y_bar)平移
    x_shift = x - x_bar;
    y_shift = y - y_bar;

    % 构造最小二乘矩阵
    A = zeros(n+1, n+1);
    b = zeros(n+1, 1);
    for i = 0:n
        for j = i:n
            A(i+1, j+1) = sum(x_shift.^(i+j)); % Vandermonde矩阵
            if i == 0
                b(j+1) = sum(y_shift .* x_shift.^j);
            end
        end
    end

    % 解Ax = b，求出拟合多项式系数
    p = A \ b;

    % 构造拟合曲线
    f = polyval(p, x);
    
    % 绘制原始数据点和拟合曲线
    plot(x, y, 'o')
    hold on
    plot(x, f)
    xlabel('x')
    ylabel('y')
    title(['Least-squares polynomial fit of degree ', num2str(n)])
    legend({'Data points', 'Fitted polynomial'})
end

使用示例：

% 生成样例数据
x = linspace(0, 2*pi, 10);
y = sin(x) + randn(size(x)) * 0.1; % 加入噪声

% 最小二乘拟合
[p, f] = polyfit_least_square(x, y, 3); % 拟合3次多项式

% 输出拟合多项式系数
disp('Polynomial coefficients:')
disp(p')

参考文献：

[1] https://en.wikipedia.org/wiki/Least_squares

[2] https://www.mathworks.com/help/matlab/ref/polyfit.html

[3] https://www.mathworks.com/help/matlab/ref/vander.html