计算步骤如下:
下面使用书中的练习y=exp(a*x^2+b*x+c)+w这个模型验证一下,其中w为噪声,a、b、c为待解算系数。
代码如下:
1 clear all; 2 close all; 3 clc; 4
5 a=1;b=2;c=1; %待求解的系数 6
7 x=(0:0.01:1)';
8 w=rand(length(x),1)*2-1; %生成噪声 9 y=exp(a*x.^2+b*x+c)+w; %带噪声的模型 10 plot(x,y,'.') 11
12 pre=rand(3,1); %步骤1 13 for i=1:1000
14
15 f = exp(pre(1)*x.^2+pre(2)*x+pre(3)); 16 g = y-f; %步骤2中的误差 17
18 p1 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x.^2; %对a求偏导 19 p2 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x; %对b求偏导 20 p3 = exp(pre(1)*x.^2+pre(2)*x+pre(3)); %对c求偏导 21 J = [p1 p2 p3]; %步骤2中的雅克比矩阵 22
23 delta = inv(J'*J)*J'* g; %步骤3,inv(J'*J)*J'为H的逆 24
25 pcur = pre+delta; %步骤4 26 if norm(delta) <1e-16
27 break; 28 end 29 pre = pcur; 30 end 31
32 hold on; 33 plot(x,exp(a*x.^2+b*x+c),'r'); 34 plot(x,exp(pre(1)*x.^2+pre(2)*x+pre(3)),'g'); 35
36 %比较一下 37 [a b c] 38 pre'
迭代结果,其中散点为带噪声数据,红线为原始模型,绿线为解算模型