图

最小二乘法的数值计算

预备知识 最小二乘法, Nelder-Mead 算法

   我们在 “最小二乘法” 见到的三个例子中, 方差函数都是待定系数的线性组合, 这种情况下我们令偏导为零后得到的是线性方程组, 便于求解. 然而当方差不是待定系数的线性组合时, 得到的方程组往往非常复杂, 这时就需要借助数值计算. 相比用数值计算解 $N$ 元的非线性方程组, 更简单的方法是直接用数值方法寻找方差函数的极小值(如 Nelder-Mead 算法) . 实践证明, 大多数情况下极小值点仅有一个, 那就是最小值点.

   为了验证结果的正确性, 我们先来用数值方法拟合 $A\cosRound {x + \varphi_0} + C$, 并与“最小二乘法” 中的方法比较结果.

图
图1:运行结果

  curveFit.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
function curveFit
close all;
% 随机生成简谐曲线
N = 20;
x0 = linspace(0, 2*pi, N);
y0 = 5*rand * sin(x0 + 2*pi * rand) + 10 * rand;
y0 = y0 + 2*rand(1,20); % 产生随机误差
scatter(x0, y0, '+'); % 画出散点
hold on;

% Nelder-Mead 求方差最小值点
f = @(x) norm( x(1)*sin(x0 + x(2)) + x(3) - y0 )^2;
c = NelderMead(f, [5, 1, 5], 1e-7);
% 画拟合结果
x = linspace(0, 2*pi, 50);
y1 = c(1) * sin(x + c(2)) + c(3);
plot(x, y1);

% 解线性方程组求方差最小值点
c = sinfit(x0, y0);
% 画拟合结果
y2 = c(1)*cos(x) + c(2)*sin(x) + c(3);
plot(x, y2, '.');
end

% 拟合简谐曲线
% y = C(1)*cos(x) + C(2)*sin(x) + C(3)
function C = sinfit(x, y)
N = numel(x);
cosx = cos(x); sinx = sin(x);
sc = sum(sinx.*cosx);
s = sum(sinx); c = sum(cosx);
% 系数矩阵
M = [sum(cosx.^2), sc          , c;
               sc, sum(sinx.^2), s;
                c,            s, N];
b = [sum(y.*cosx); sum(y.*sinx); sum(y)];
C = M\b; % 解线性方程组
end

   运行结果如图 1 所示, 可见两种方法拟合出的曲线一致(红色的曲线和黄色的点线). 注意第 13 行使用了“Nelder-Mead 算法” 中的函数 NelderMead 求函数句柄 f 的最小值.

致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利