Using curvefit for conditional functions

I have some data and I want to fit it with a function. Here is my attepmt

clear all
clc
X1=[-23.24 -16.79 -14.4 -10.39 -7.494 -4.984 -2.75 0 0 0 23.24 16.79 14.4 10.39 7.494 4.984 2.75];
Y1=[-4 -3 -2.6 -2 -1.6 -1.3 -1.1 -1 0 1 4 3 2.6 2 1.6 1.3 1.1];

plot(X1,Y1,'rO')
grid on
xlabel('X-Data')
ylabel('Y-Data')

Now I want to customize this plot with the function:

I have the x-data but I have to find a and b.

I used matlab lsqcurvefit, but failed! First I tried to make the function in a separate file:

function Vx=V2(x,u)
d = abs(x)<=1;
Vx(d) = 0;
d = abs(x)>1;
Vx(d) = u(1)*u(2)*sqrt((x(d).^2)-1);

end

And finally I added this code to my first code:

U=[2.3 3.5];
func = V2(X1,U);
Xf=lsqcurvefit(func,U,X1,Y1)

But I failed! How can I fit these data to the function?

Answer 1

[John D'Errico]

You CANNOT fit that model. Why not? a and b are not independent parameters. They look as if they are.

Ignoring the part where x is small, your model looks like:

a*b*sqrt(x^2 - 1)

First, I'll point out that it does not even fit your data. When x is near 1, your model predicts ZERO. Yet your data passes through 1 or -1 when x is near zero. So your model does not even have the correct behavior for x near 1. Look at the curve you want to use. a and b are irrelevant, as they just scale the thing up and down.

fun = @(x) sqrt(x.^2-1)

fun = function_handle with value:
    @(x)sqrt(x.^2-1)

fplot(fun,[1,5])
grid on
xlabel X
ylabel Y

As you should see, your data is flipped around. While both curves seems to be roughly asymptotically linear for large x, your data has a roughly zero slope near x == 1, whereas the model you want to use has an infinite slope near x==1.

But disregarding the fact that your model is complete crap for this data, could you even solve for both a and b?

NO!

Suppose you did find an optimal value for a and b. Now, consider what would happen if you tried instead a*2 and b/2 in that model? You would then notice that:

a*b*sqrt(x^2-1) == (a*2)*(b/2)*sqrt(x^2-1)

so if you think you have found an "optimal" set of values for a and b, note that they are not optimal. There are infinitely many possible sets of numbers you could have chosen, all equally optimal.

And that means you can arbitrarily set b==1. Now, in theory, you could hope to estimate a.

Except of course, that this model will never look like that data. lsqcurvefit will not help to fit a square peg into a round hole.

(Sigh. Do I need to tell you how to fix this? Probably.)

You seem to have the x and y relationships swapped around. That is, suppose you treat this as a problem of fitting x as a function of y? First, plot your data.

X1=[-23.24 -16.79 -14.4 -10.39 -7.494 -4.984 -2.75 0 0 0 23.24 16.79 14.4 10.39 7.494 4.984 2.75];
Y1=[-4 -3 -2.6 -2 -1.6 -1.3 -1.1 -1 0 1 4 3 2.6 2 1.6 1.3 1.1];
plot(Y1,X1,'rO')
xlabel Y
ylabel X
hold on
fplot(@(Y) 6*Y,'b')
hold off
grid on

On top of that, I can see the curve should be visually roughly asymptotic to the line X = 6*Y. I've plotted that line in blue. And that should tell you the rough value of a, as 6.

Now, consider the model as:

if abs(y) < 1, x == 0

if abs(y) >= 1, x = a*sqrt(y^2 - 1)

this has the desired behavior. When y == +/-1, we can see that x will be zero. And it is still asymptotially linear for large y. Finally, see that I dropped b completely. Assume b==1.

Can I fit that model? Yes. I'll use the curvefitting toolbox here, the I'll do it with lsqcurvefit, since you want to use that tool.

As well, I'll show you how to write that piecewise thing using only one line, using max.

mdl = fittype('a*sign(Y1).*sqrt(max(abs(Y1),1).^2 - 1)','indep','Y1')

mdl = 
     General model:
     mdl(a,Y1) = a*sign(Y1).*sqrt(max(abs(Y1),1).^2 - 1)

You should see that when abs(Y1) is 1 or larger, this function will have the shape a*sqrt(Y1^2 - 1). When abs(Y1) is 1 or less, the max function maps that to 1. But then it subtracts 1, to get 0. So this function now has the desired shape.

The sign function is NOT the trigonometric sine function.

estmdl = fit(Y1',X1',mdl,'start',1)

estmdl = 
     General model:
     estmdl(Y1) = a*sign(Y1).*sqrt(max(abs(Y1),1).^2 - 1)
     Coefficients (with 95% confidence bounds):
       a =       5.985  (5.971, 6)

So your model is now

if abs(y) < 1, x == 0

if abs(y) >= 1, x = sign(y)*5.985*sqrt(y^2 - 1)

plot(estmdl,'b')
hold on
plot(Y1,X1,'ro')
xlabel Y1
ylabel X1
grid on
hold off

As you should see, it does a remarkably good job now. lsqcurvefit will do as easily. (Do I need to show you how to do it with lsqcurvefit too?) But don't bother trying to estimate b. You cannot do so. b == 1.

[Pouyan Msgn]

OK Thanks for the advice!

But how can I define this function in Matlab ? I mean this function with condition: