| 查看: 841 | 回复: 2 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
LM算法碰到分段函数该怎么处理 已有1人参与
|
||
|
各位虫友大家好，本虫在目前在用LM算法进行函数的曲线拟合，函数分两段，一部分类似于余弦，另一部分类似于幂函数。在循环中，因为分段点也是一个未知的变量，所以不能进行比较大小，特地求助广大虫友，感激不尽。这个问题卡了我好久。谢谢。下面附上一些代码。 代码1：单一函数的拟合，为了虫友能更好理解LM算法。 % 计算函数f的雅克比矩阵，是解析式 syms a b c d y x real; f=a+b*cos((pi/c)*(x-d)); Jsym=jacobian(f,[a b c d]) % 拟合用数据。 data_1=[0.25 0.5 1 1.5 2 3 4 6 8]; obs_1=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; % 2. LM算法 % 初始猜测s a0=4; b0=20;c0=14;d0=13; y_init = a0+b0*cos((pi/c0)*(x-d0)); % 数据个数 Ndata=length(obs_1); % 参数维数 Nparams=4; % 迭代最大次数 n_iters=50; % LM算法的阻尼系数初值 lamda=0.01; % step1: 变量赋值 updateJ=1; a_est=a0; b_est=b0; c_est=c0; d_est=d0; % step2: 迭代 for it=1:n_iters if updateJ==1 % 根据当前估计值，计算雅克比矩阵 J=zeros(Ndata,Nparams); J for i=1:length(data_1) i J(i,  =[ 1 cos((pi*(d_est - data_1))/c_est) (b_est*pi*sin((pi*(d_est - data_1))/c_est).*(d_est - data_1))/c_est^2 -(b_est*pi*sin((pi*(d_est - data_1))/c_est))/c_est];end % 根据当前参数，得到函数值 y_est =a_est+b_est*cos((pi/c_est)*(data_1-d_est)) ; % 计算误差 d=obs_1-y_est; % 计算（拟）海塞矩阵 H=J'*J; % 若是第一次迭代，计算误差 if it==1 e=dot(d,d); end end % 根据阻尼系数lamda混合得到H矩阵 H_lm=H+(lamda*eye(Nparams,Nparams)); % 计算步长dp，并根据步长计算新的可能的\参数估计值 dp=inv(H_lm)*(J'*d( );g = J'*d( ;a_lm=a_est+dp(1); b_lm=b_est+dp(2); c_lm=c_est+dp(3); d_lm=d_est+dp(4); % 计算新的可能估计值对应的y和计算残差e y_est_lm =a_lm+b_lm*cos((pi/c_lm)*(data_1-d_lm)) ; d_lm=obs_1-y_est_lm; e_lm=dot(d_lm,d_lm); % 根据误差，决定如何更新参数和阻尼系数 if e_lm<e lamda=lamda/10; a_est=a_lm; b_est=b_lm; c_est=c_lm; d_est=d_lm; e=e_lm; disp(e); updateJ=1; else updateJ=0; lamda=lamda*10; end end %显示优化的结果 a_est b_est c_est d_est 代码2（部分）：对代码1进行改进，想求分段函数的偏导。 % 计算函数f的雅克比矩阵，是解析式 syms a b c d f g x k y1 y2 real; if x<f y1=a+b*cos((pi/c)*(x-d)); y1x=jacobian(y1,[a b c d]); else y2=(a+g)+[b*cos((pi/c)*(f-d))-g]*exp(-(x-f)/k); y2x=jacobian(y2,[a b c d g ]) end %Jsym=jacobian(f,[a b c d f g h]); % 拟合用数据。 data_1=[0.25 0.5 1 1.5 2 3 4 6 8]; obs_1=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; % 2. LM算法 % 初始猜测s a0=4; b0=20;c0=14;d0=13;f0=17;g0=0.9; y_init = a0+b0*cos((pi/c0)*(x-d0)); % 数据个数 Ndata=length(obs_1); % 参数维数 Nparams=4; % 迭代最大次数 n_iters=50; % LM算法的阻尼系数初值 lamda=0.01; % step1: 变量赋值 updateJ=1; a_est=a0; b_est=b0; c_est=c0; d_est=d0; f_est=f0; g_est=g0; % step2: 迭代 for it=1:n_iters if updateJ==1 % 根据当前估计值，计算雅克比矩阵 J=zeros(Ndata,Nparams); a=ones(9,1); if x<h J=[ a cos((pi*(d_est - data_1))/c_est)' ((b_est*pi*sin((pi*(d_est - data_1))/c_est).*(d_est - data_1))/c_est^2)' (-(b_est*pi*sin((pi*(d_est - data_1))/c_est))/c_est)'] else J=[ a cos((pi*(d - f))/c)*exp((f - x)/k), (b*pi*exp((f - x)/k)*sin((pi*(d - f))/c)*(d - f))/c^2, -(b*pi*exp((f - x)/k)*sin((pi*(d - f))/c))/c, 1 - exp((f - x)/k), (b*pi*exp((f - x)/k)*sin((pi*(d - f))/c))/c] end 再次感谢虫友，谢谢！ |
回复此楼» 猜你喜欢
2026年循环经济功能材料国际会议（ICFMCE 2026）
已经有0人回复
-
2026年第五届电气、电子与信息工程国际会议（ISEEIE 2026）
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有137人回复
-
哈尔滨理工大学物理系招收物理学考研调剂
已经有5人回复
-
0702一志愿吉大B区求调剂
已经有5人回复
-
求调剂
已经有0人回复
-
0702一志愿吉大B区求调剂有论文
已经有0人回复
-
请问还有没有用Latex写文章的小伙伴们？
已经有0人回复
-
光学工程学硕调剂信息
已经有27人回复
-
欢迎加入课题组
已经有0人回复
-
散金币，求好运，祝面上顺利！
已经有43人回复
3楼2017-10-16 10:31:21
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
2楼2017-10-11 09:19:32
