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
