How can I fit data to a piecewise function, where the breakpoint of the function is also a parameter to be optimised?

I have data with x and y values. This data should conform to a function: an assymmetric parabola. Here, the parameters that define the shape of the parabola should be different on either side of the maximum point of the parabola i.e. the breakpoint is where the maximum value of y occurs.

I was hoping to use 'fit' and to define an anonymous function for my data. But I'm not able to work out how to define an anonymous, piecewise function, especially where the breakpoint is one of the parameters to be determined by the fitting procedure, as it is not immediately clear from the data itself where the maximum value of y should occur.

Any help would be appreciated.

Answer 1

[Matt J]

Once you've chosen the coefficients of the first parabola [a1,b1,c1], the breakpoint is determined from,

d=-b1/(2*a1)

Only the leading coefficient of the second parabola is a free parameter:

F=@(x) asymParabola(-2,1,0,-0.6,x);
fplot(F,[-10,10]);axis padded  %example plot

ft = fittype(@(a1,b1,c1,a2, x) asymParabola(a1,b1,c1,a2, x) )

ft = 
     General model:
     ft(a1,b1,c1,a2,x) = asymParabola(a1,b1,c1,a2,x)

function y=asymParabola(a1,b1,c1,a2, x)

   d=-b1/(2*a1);
   b2=-d*2*a2;

   c2=polyval([a1,b1,c1],d)-polyval([a2,b2,0],d);
   
   left=(x<=d);
   
   y=x;
   y(left)=polyval([a1,b1,c1],x(left));
   y(~left)=polyval([a2,b2,c2],x(~left));

end

[Torsten]

From what the OP wrote, the breakpoint where the piecewise function attains its maximum doesn't need to be the maximum of one of the parabolas - it only has to be the maximum y-value. But maybe I didn't interpret it correctly.

[Matt J]

That would be very strange. If the given x,y values have no noise/errors, then the maximum y-value and the maximum of the parabola are one and the same. To say that you now want the fit to change when y is corrupted by noise is usually contrary to the goals of a fit.

[Rahul]

Thanks a lot Matt, I think this is right; I think there should only be 4 free parameters to be constrained overall. This simplifies things greatly.

[Torsten]

@Matt J

If the given x,y values have no noise/errors, then the maximum y-value and the maximum of the parabola are one and the same.

Why ? Both parabola can intersect below their respective maxima, and nonetheless the point of intersection can be the maximum y-value of the piecewise function.

But it seems you interpreted the question correctly.

Answer 2

[Walter Roberson]

(a1*x.^2 + b1*x + c1) .* (x <= d) + (a2*x.^2 + b2*x + c2) .* (x > d)

Note that for this to work, the coefficients must be constrained to be finite

[Paul]

Sounds like both sides of the function should have the same value at x = d, at least that's how interpret the question. If so, then I think the function would look something like

(a1*(x-d).^2 + b1*(x-d) + c) .* (x <= d) + (a2*(x-d).^2 + b2*(x-d) + c) .* (x > d)

[Rahul]

Thanks Paul, I think this would also work, but as Matt pointed out I think some of extra parameters are not needed (as they are actually constrained by the others). I am not sure how Matlab deals with this in the curve fitting/optimisation algorithms.

Answer 3

[Catalytic]

You can also parametrize the model function directly in terms of the break point coordinates (xbreak, ybreak) and two curvature parameters -

F= @(a1,a2,xbreak,ybreak, x) modelFun(a1,a2,xbreak,ybreak, x);

xbreak=3; ybreak=5;
fplot( @(x) F(-2,-0.6,xbreak,ybreak,x), [1,5]); 
xline(xbreak,'--')

fType = fittype(F);

function y=modelFun(a1,a2,xbreak,ybreak, x)

   X=x-xbreak;
   LHS=(X<=0);
   RHS=~LHS;
   
   y=X.^2;
   y(LHS)=a1.*y(LHS) + ybreak;
   y(RHS)=a2.*y(RHS) + ybreak;

end

Answer 4

[Matt J]

Why ? Both parabola can intersect below their respective maxima, and nonetheless the point of intersection can be the maximum y-value of the piecewise function.

If I understand @Torsten, that would be a 6-parameter function,

F=@(x) asymParabola(-2,1,0,-6,5,-20 ,x);
fplot(F,[-10,-1]);axis padded  %example plot

function y=asymParabola(a1,b1,c1, a2, s, rightSlope,  x)
%Requirements: a1<0, a2<0, s>=0, m<=0
   

   d=-b1/(2*a1)-s;

   c2=polyval([a1,b1,c1],d);
   
   left=(x<=d); 
   right=~left;
   xright=x(right);
   
   y=x;
   y(left)=polyval([a1,b1,c1],x(left));
   y(right)=a2*(xright-d).^2  + rightSlope*(xright-d) +c2;
   
end