lsqcurvefit with 2 points trouble

Hello Matlab Community,

I need help with a little code ...

I want to do this : , obtain by this code :

xdata2 = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata2 = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
fun = @(x,xdata2)x(1)-x(2)*xdata2(1:2).^2;
x0 = [100,-1];
x = lsqcurvefit(fun,x0,xdata2,ydata2(1:2))

times = linspace(xdata2(1),xdata2(end));
figure
hold on
plot(xdata2,ydata2,'ko')
plot(times,x(1)-x(2)*times.^2,'b--')
ylim([-10 500])
xlim([0 100])
legend('Data','Fitted curve')
title('Data and Fitted Curve')

with my own data but I obtain this : with the following code:

xdata2 = 2e-5:2e-5:18e-5;
ydata2 = [0.997 0.991 0.983 0.974 0.964 0.954 0.944 0.935 0.926];
fun = @(x,xdata2)x(1)-x(2)*xdata2(1:2).^2;
x0 = [100,-1];
x = lsqcurvefit(fun,x0,xdata2,ydata2(1:2))

times = linspace(xdata2(1),xdata2(end));
figure
hold on
plot(xdata2,ydata2,'ko')
plot(times,x(1)-x(2)*times.^2,'b--')
legend('Data','Fitted curve')
title('Data and Fitted Curve')

. I expect a line which passes by at least the 2 first points... Do you know why please?

Thanks in advance.

Regards,

[Mathieu NOE]

hello

I am not sure to understandwhat shape / model you want to fit

from your function description, looks like you want to fit an inverted parabola ? in the first picture the legend speaks about an exponential fit , so what do you want ?

[Matt J]

It is strange that you don't just use polyfit. There are no bounds being applied in your fitting process.

Answer 1

[Mathieu NOE]

maybe you want a parabola , you can use polyfit to obtain the best fit :

x = 2e-5:2e-5:18e-5;
y = [0.997 0.991 0.983 0.974 0.964 0.954 0.944 0.935 0.926];
% y = a - b*x²
% Fit a polynomial p of degree "degree" to the (x,y) data:
degree = 2;
p = polyfit(x,y,degree);

% Evaluate the fitted polynomial p and plot:
f = polyval(p,x);
eqn = poly_equation(flip(p)); % polynomial equation (string)
Rsquared = my_Rsquared_coeff(y,f); % correlation coefficient

figure(3);plot(x,y,'o',x,f,'-')
legend('data',eqn)
title(['Data fit - R squared = ' num2str(Rsquared)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function eqn = poly_equation(a_hat)
eqn = " y = "+a_hat(1);
for i = 2:(length(a_hat))
    if sign(a_hat(i))>0
        str = " + ";
    else
        str = " ";
    end
    if i == 2
        eqn = eqn+str+a_hat(i)+"*x";
    else
        eqn = eqn+str+a_hat(i)+"*x^"+(i-1)+" ";
    end
end
eqn = eqn+" ";
end
        

[Mathieu NOE]

FYI, maybe this decaying exponential model (for the first data set)

x = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
y = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];

f = @(a,b,c,x) a.*exp(-(x./b).^c);
obj_fun = @(params) norm(f(params(1), params(2), params(3),x)-y);
sol = fminsearch(obj_fun, [y(1),1,1]);
a = sol(1)

a = 495.4161

b = sol(2)

b = 10.0348

c = sol(3)

c = 1.0200

yfit = f(a, b,c, x);
Rsquared = my_Rsquared_coeff(y,yfit); % correlation coefficient

figure;
plot(x, y, '+', 'MarkerSize', 10, 'LineWidth', 2)
hold on
plot(x, yfit, '-');grid on
xlabel('time t');
ylabel('P(t)');

title(['Data fit - R squared = ' num2str(Rsquared)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

[Mindy]

Thank you so much for your help.

Sorry I wasn't enough clear but yes, I want to use an osculating parabola on my first two points from my time series. I forget to update the legend ("fitted exp" -> "fitted curve"). Thanks for sharing the one with the exp. fitting also. It will be usefull for my next step!

[Mathieu NOE]

My pleasure !

Answer 2

[Matt J]

Change the units of your xdata to something more natural:

xdata2 = (2e-5:2e-5:18e-5)*1e5;

and then,

ydata2 = [0.997 0.991 0.983 0.974 0.964 0.954 0.944 0.935 0.926];
fun = @(x,xdata2)x(1)-x(2)*xdata2(1:2).^2;
x0 = [100,-1];
x = lsqcurvefit(fun,x0,xdata2,ydata2(1:2))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 1×2

0.9990 0.0005

times = linspace(xdata2(1),xdata2(end));
figure
hold on
plot(xdata2,ydata2,'ko')
plot(times,x(1)-x(2)*times.^2,'b--')
legend('Data','Fitted curve')
title('Data and Fitted Curve')

[Matt J]

You can, of course, change xdata2 back to their original units by multiplying x(2) by 1e10. Don't know why you would prefer the original units, however.

Alternatively, you can introduce a scaling coefficient in your model equation:

xdata2 = (2e-5:2e-5:18e-5);
ydata2 = [0.997 0.991 0.983 0.974 0.964 0.954 0.944 0.935 0.926];
fun = @(x,xdata2)x(1)-1e10*x(2)*xdata2.^2;
x0 = [100,-1];
x = lsqcurvefit(fun,x0,xdata2(1:2),ydata2(1:2))

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 1×2

0.9990 0.0005

times = linspace(xdata2(1),xdata2(end));
figure
hold on
plot(xdata2,ydata2,'ko', times,fun(x,times),'b--')
legend('Data','Fitted curve')
title('Data and Fitted Curve')
hold off

[Mindy]

Thank you so much for your help. xdata is expressed in seconds. These are my 10 first values ​​of my dataset (xdata = times, ydata = velocity from cross-correlation function). In my case, acquisitions of this frequency level are useful for calculating the size of turbulence.