How to fit piecewise linear

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

Answer 1

[John]

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