一杯清酒邀明月
天下本无事,庸人扰之而烦耳。

计算步骤如下:

下面使用书中的练习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'

迭代结果,其中散点为带噪声数据,红线为原始模型,绿线为解算模型

posted on 2020-09-10 14:42  一杯清酒邀明月  阅读(1170)  评论(0编辑  收藏  举报