我有一些实验数据(信号与时间):
我想把这些曲线拟合成线性分段函数。我需要它,因为:
1) 我对“跳跃”的时间感兴趣(中间部分)
2) 我想整理一下右上角的曲线。
要进行拟合,我尝试使用解决方案,如中所述
stackoverflow
,但这个方法(或至少我的实现)对初始节点位置非常敏感。我的代码和结果:
%% Generate dummy data
x=linspace(-2,2,100); y=sinh(x)./cosh(x).^2;
% add noise
y=y+(rand(1,numel(y))-0.5)*0.1;
%% Estimate knots
d=(max(x)-min(x));
X0_badEstimation=[min(x)+d/5, min(x)+4*d/5];
X0_goodEstimation=[min(x)+d/3, min(x)+2*d/3];
%% Estimate piecewise fit
idx=1:min(10,numel(x));
p1=polyfit(x(idx),y(idx),1);
idx=max(0,round(numel(x)/2-10)):min(round(numel(x)/2+10),numel(x));
p2=polyfit(x(idx),y(idx),1);
idx=max(0,(numel(x)-10)):numel(x);
p3=polyfit(x(idx),y(idx),1);
%% estimate slopes
s1=p1(1); s2=p2(1); s3=p3(1);
%% estimate offsets
o1=p1(2); o2=p2(2); o3=p3(2);
%% model in form
% y=(o1+s1*x)*((x0-x)>=0)+(o2+s2*x)*((x0-x)<0)*((x0+delta)-x>=0)+(o3+s3*x)*(((x0+delta)-x)<0)
model=@(P,x) (P(6)+P(3).*x).*((P(1)-x)>=0)+...
(P(7)+P(4).*x).*((P(1)-x)<0).*((P(1)+P(2))-x>=0)+...
(P(8)+P(5).*x).*(((P(1)+P(2))-x)<0);
%% Initial parameters:
x0=X0_goodEstimation(1);
delta=diff(X0_goodEstimation);
P0_goodEstimation=[x0,delta,s1,s2,s3,o1,o2,o3];
x0=X0_badEstimation(1);
delta=diff(X0_badEstimation);
P0_badEstimation=[x0,delta,s1,s2,s3,o1,o2,o3];
%% fit it!
Pfit_goodEstimation = lsqcurvefit(model,P0_goodEstimation,x,y);
Pfit_badEstimation = lsqcurvefit(model,P0_badEstimation,x,y);
%% plot results
hold all
plot(x,y,'LineWidth',3)
plot(x,model(Pfit_badEstimation,x))
plot(x,model(Pfit_goodEstimation,x))
legend({'data','bad knots estimation','good knots estimation'},'FontSize',24)
关于代码的几点评论:
1) 我用来生成虚拟数据的函数:
y=sinh(x)./cosh(x).^2
对我来说没有意义,我是凭经验发现的,只用于提供运行的代码。
2) 我尝试了不同的方法来更好地估计结的位置,但在实际数据上我没有成功,所以这里我使用一些简单的方法JU进行演示。
问题:
你能帮我试穿一下吗?我在实现中是否有错误,或者在我的情况下应该使用不同的方法?
