24小时热门版块排行榜

返回列表

manmanbobo

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 25
红花: 1
帖子: 37
在线: 12.5小时
虫号: 2965461
注册: 2014-02-09
专业: 地理信息系统

[求助] 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
再次感谢虫友，谢谢！

回复此楼