How to fit piecewise linear

7 views (last 30 days)
John
John on 10 Oct 2017
Edited: John on 11 Oct 2017
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

Accepted Answer

John
John on 11 Oct 2017
Edited: John on 11 Oct 2017
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

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!