How to fit piecewise linear

Question

John 2017-10-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360656-how-to-fit-piecewise-linear

编辑： John 2017-10-11

The code below generates data similar to experimental data. However, the curve fit does not appear to locate the knot point. Any suggestions?

function bilinearfit
close all; clc;
% Random coeffs
al = -rand*20;
ah = rand/10;
bl = (rand-0.5)*20;
isxnrand = 1;
if isxnrand
    xn = rand*8+2;
    bh = xn * (al - ah) + bl;
else
    bh = (rand-0.5)*20;
    xn = (bh - bl) / (al - ah);
end
% Random data...
xdata = linspace(2,10,101);
blfit = @(coeff, x) ( x < coeff(1) ) .* (coeff(2) * x + coeff(3)) ...
    + (~( x < coeff(1) )) .* (coeff(4) * x + coeff(5));
coeffNoise = ((rand(1,5)-0.5)*2/100+1);
coeff0 = [xn al bl ah bh] .* coeffNoise;
ydata = blfit(coeff0, xdata);
dataNoise = (rand(size(ydata))-0.5)*2/100+1;
ydata = ydata .* dataNoise;
plot(xdata,ydata,'ko')
% F = @(B,xdata) min(B(1),B(2)+B(3)*xdata);   %Form of the equation
F = @(coeff)(blfit(coeff,xdata));
IC = [1 1 1 1 1]; %Initial guess
LB = [min(xdata) -inf -inf -inf -inf];
UB = [max(xdata) inf inf inf inf];
B = lsqcurvefit( blfit,IC,xdata,ydata,LB,UB);
hold all;
plot(xdata,blfit(B,xdata),'r--');
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John 2017-10-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/360656-how-to-fit-piecewise-linear#answer_285228

编辑：John 2017-10-11

在 MATLAB Online 中打开

The problem appears to be by using five unknowns the problem is ill posed. By using 4 unknowns, the calculated line fits the data. For completeness, the bootstrap method is used to estimate the 90% confidence range of the unknown parameters.

function bilinearfit
% given xn, yn
% x < xn
%  y = al x + bl
%  yn = al xn + bl
%  bl = yn - al xn
%  y = al x + yn - al xn
%  y = al (x - xn) + yn
% x >= xn
%  y = au x + bu
%  yn = au xn + bu
%  bu = yn - au xn
%  y = au x + yn - au xn
%  y = au (x - xn) + yn
close all; clc;
ntrys = 10; npoints = 30;
blfit = @(coeff, x) ( x < coeff(1) ) .* (coeff(3) * (x-coeff(1)) + coeff(2)) ...
    + (~( x < coeff(1) )) .* (coeff(4) * (x-coeff(1)) + coeff(2));
lsqOpts = optimset('Display','off');
for itry=1:ntrys
    % Random coeffs
    al = -rand*9-1;  % Between -1 and -10
    ah = (rand-0.5)*2/10; % Between -0.1 and 0.1
    xn = rand*6+2;   % Between 2 and 8
    yn = rand+.2; % Between 0.2 and 1.2
      % Random data...
      xMeasValues = linspace(2,10,6);
      indxMeasValues = randi([1 length(xMeasValues)],1,npoints);
      xMeas = xMeasValues(indxMeasValues);
      coeffNoise = ((rand(1,4)-0.5)*2/100+1);
      coeffAct = [xn yn al ah];
      coeffWNoise =  coeffAct .* coeffNoise;
      yTrue = blfit(coeffWNoise, xMeas);
      dataNoise = ((rand(size(yTrue))-0.5)*2)*(30/100*max(yTrue));
      yMeas = yTrue + dataNoise;
      nBootStrap = 100;
      percentBootStrap = 90;
      nPointsBootStrap = floor(percentBootStrap/100*npoints);
      coeffBootStrap = nan(nBootStrap,4);
      for iBootStrap = 1:nBootStrap
          indxMeasValues = randi([1 length(xMeas)],1,nPointsBootStrap);
          selXMeas = xMeas(indxMeasValues);
          selYMeas = yMeas(indxMeasValues);
          IC = [(min(selXMeas)+max(selXMeas))/2 1 1 1]; %Initial guess
          LB = [min(selXMeas) -inf -inf -inf];
          UB = [max(selXMeas) inf inf inf];
          coeffBootStrap(iBootStrap,:) = lsqcurvefit( blfit,IC,selXMeas,selYMeas,LB,UB,lsqOpts);
      end
      IC = [(min(xMeas)+max(xMeas))/2 1 1 1]; %Initial guess
      LB = [min(xMeas) -inf -inf -inf];
      UB = [max(xMeas) inf inf inf];
      coeffEst = lsqcurvefit( blfit,IC,xMeas,yMeas,LB,UB,lsqOpts);
      figure;
      set(gcf,'position',[140.2000   66.6000  908.0000  695.2000]);
      subplot(3,1,1);    
      plot(xMeas,yMeas,'k.')
      hold all;
      xfit = [min(xMeas) coeffEst(1) max(xMeas)];
      yfit = blfit(coeffEst,xfit);
      plot(xfit, yfit, 'r-', 'LineWidth', 1.5);
      plot(coeffEst(1), coeffEst(2), 'ro');
      xlabel('x'); ylabel('y');
      title(sprintf('Try %d: xn %.2f yn %.2f al %.3f au %.3f', ...
          itry, coeffWNoise));
      for ivar=1:4
          subplot(3,2,2+ivar);
          hist(coeffBootStrap(:,ivar),20);
          switch(ivar)
              case 1
                  varname = 'xn';
              case 2
                  varname = 'yn';
              case 3
                  varname = 'al';
              case 4
                  varname = 'au';
          end
          prct = .9;
          prctRng = [(1-prct)/2 1-(1-prct)/2;];
          indxRng = [ceil(prctRng(1)*nBootStrap) floor(prctRng(2)*nBootStrap)];
          sortCoeff = sort(coeffBootStrap(:,ivar));
          mnBSCoeff = mean(coeffBootStrap(:,ivar));
          rngBSCoeff = sortCoeff(indxRng);
          % deltaBSCoeff = rngBSCoeff - mnBSCoeff;
          titstr = sprintf('Actual %.3f mn %.3f %d%%: [%.3f, %.3f]', ...
              coeffWNoise(ivar), mnBSCoeff, round(prct*100), rngBSCoeff);
          title(titstr);
          xlabel(varname);
          ylabel('# Occur');
      end
      drawnow;
  end
  end