最小二乘法拟合多项式（Matlab）

贡献者： addis

预备知识　最小二乘法

　　事实上，Matlab 的 polyfit 函数已经可以实现多项式拟合，但为了教学我们重新实现一次。

图 1：15 阶多项式拟合正弦函数（蓝点为数据点，红线为多项式）

图 2：20 阶多项式拟合锯齿波（蓝线为数据点，红线为多项式）

　　以下代码使用式 6 解出多项式拟合的最小二乘法系数 $c_i$。

\begin{equation} p(x) = c_0 + c_1 x + c_2 x^2 + \dots~ \end{equation}

代码 1：ls_poly.m

% 多项式最小二乘法拟合
% p(x) = c(1)x^order + c_2 x^(order-1) + ... c(order+1)
% numel(x) = numel(y)
function c = ls_poly(x, y, order)
Sx = zeros(1, 2*order+1);
Sx(1) = numel(x);
for n=1:2*order
    Sx(n+1) = sum(x.^n);
end

Sxy = zeros(order+1, 1);
for n=0:order
    Sxy(n+1) = sum(y .* x.^n);
end
Sxy = flip(Sxy);

A = zeros(order+1, order+1);
for i = 0:order
    for j = 0:order
        A(order-i+1,order-j+1) = Sx(i+j+1);
    end
end
c = A \ Sxy; % 解线性方程组
end

　　使用例子：拟合正弦函数。

代码 2：ls_poly_demo.m

% ls_poly_demo
x = linspace(-10, 10, 100);
y = sin(x);
figure; plot(x, y);
c = ls_poly(x, y, 15);
y1 = polyval(c, x); % 求值
hold on; plot(x, y1);
axis([-10, 10, -1, 1]);

运行结果如图 1 。

代码 3：ls_poly_demo2.m

% ls_poly_demo2
x = linspace(-1, 1, 50);
y = zeros(size(x));
y(x > 0) = -x(x > 0);
y(x > 0.5) = x(x > 0.5)-1;
y(x < 0) = x(x < 0);
y(x < -0.5) = -x(x < -0.5)-1;
figure; plot(x, y, '.');
c = ls_poly(x, y, 20);
y1 = polyval(c, x);
hold on; plot(x, y1);

致读者：小时百科一直以来坚持所有内容免费无广告，这导致我们处于严重的亏损状态。长此以往很可能会最终导致我们不得不选择大量广告以及内容付费等。因此，我们请求广大读者热心打赏 ，使网站得以健康发展。如果看到这条信息的每位读者能慷慨打赏 20 元，我们一周就能脱离亏损，并在接下来的一年里向所有读者继续免费提供优质内容。但遗憾的是只有不到 1% 的读者愿意捐款，他们的付出帮助了 99% 的读者免费获取知识，我们在此表示感谢。