How to fit data on a specific region, using lsqcurvefit or any other?

I am trying to use lsqcurvefit to fit 4 parameters, I am usually able to do it but in this problem, the fitting function is only valid close to the plateau values. So I selected the values over which I apply the fitting. for example (14:37) in line 12.

The function seem to be a bit sensitive to initial guesses and that is not allowing me to achieve a good fit.

I wonder if the problem is that I am using lsqcurvefit and I should be using a different tool or if I need to create an iteration of initial guesses so that I can finally get the correct combination?

It is acceptable to increase or decrease a bit the range of the fitting function helps to get a better fit. Right now it is set for 14:37 for data set 1 and 12:37 for data set 2.

function:

code:

clear; clc; clf; close all;

x1 = [0.0005 0.0007 0.0008 0.001 0.0013895 0.0019307 0.0026827 0.0037276 0.0051795 0.0071969 0.01 0.013895 0.019307 0.026827 0.037276 0.051795 0.071969 0.1 0.13895 0.19307 0.26827 0.37276 0.51795 0.71969 1 1.3895 1.9307 2.6827 3.7276 5.1795 7.1969 10 13.895 19.307 26.827 37.276 51.795 71.969 100 138.95 193.07 268.27 372.76 517.95 719.69 1000];
y1 = [0.011171 0.015461 0.017605 0.021901 0.029975 0.041604 0.057316 0.078645 0.10726 0.14398 0.18833 0.24482 0.31533 0.38492 0.45433 0.51514 0.53269 0.54151 0.56182 0.56775 0.56635 0.56416 0.55786 0.55369 0.53971 0.52712 0.51333 0.50512 0.49684 0.48869 0.50622 0.51273 0.53057 0.55463 0.61141 0.78976 0.96969 1.1629 1.4075 1.6804 2.1478 2.6287 3.5048 4.4863 5.6525 7.3741];

x1_fit = x1(12:37);
y1_fit = y1(12:37);

x2 = [0.0005 0.0007 0.0008 0.001 0.0013895 0.0019307 0.0026827 0.0037276 0.0051795 0.0071969 0.01 0.013895 0.019307 0.026827 0.037276 0.051795 0.071969 0.1 0.13895 0.19307 0.26827 0.37276 0.51795 0.71969 1 1.3895 1.9307 2.6827 3.7276 5.1795 7.1969 10 13.895 19.307 26.827 37.276 51.795 71.969 100 138.95 193.07 268.27 372.76 517.95 719.69];
y2 = [0.0098762 0.013754 0.015562 0.019716 0.027332 0.038068 0.052772 0.073219 0.10152 0.13964 0.19046 0.25666 0.43342 0.53389 0.62881 0.7047 0.74682 0.73801 0.68675 0.64029 0.60505 0.59568 0.59975 0.61275 0.6274 0.64335 0.6483 0.64144 0.6314 0.64113 0.66358 0.70689 0.74901 0.79417 0.86603 0.93868 1.1244 1.3195 2.0893 2.6337 3.0766 3.6119 4.7272 6.0726 7.5609];

x2_fit = x2(14:37);
y2_fit = y2(14:37);

scale1 = 38.50;
scale2 = 48.27;

C1_guess = [1 1 0.001 500]; 
C2_guess = [1 1 0.001 500]; 


figure
subplot(1,2,1)
fun = @(C,x) C(1)*sinh(C(2)+(log10(x) - C(3)))+ C(4);
C_1 = lsqcurvefit(fun, C1_guess, x1_fit, y1_fit);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

a = loglog(x1*scale1, y1,'ks', 'MarkerFaceColor', 'c');
hold on
a1 = loglog(x1_fit*scale1, fun(C_1,x1_fit), 'r-', 'LineWidth', 1);
xlabel('x1')
ylabel('y1')

subplot(1,2,2)
C_2 = lsqcurvefit(fun, C2_guess, x2_fit, y2_fit);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

b = loglog(x2*scale2, y2,'ks', 'MarkerFaceColor', 'g');
hold on
b1 = loglog(x2_fit*scale2, fun(C_2,x2_fit), 'r-', 'LineWidth', 1);
xlabel('x2')
ylabel('y2')

[Matt J]

What makes you suppose the fit isn't good?

[Alfredo Scigliani]

The fitting function is supposed to pass through the data points almkst perfectly, lets say a deviation about 0.001

[Alfredo Scigliani]

And if you extract the second output of lsqcurvefit (which is the rsb) it is about 0.20

[Matt J]

The fitting function is supposed to pass through the data points almkst perfectly, lets say a deviation about 0.001.

That cannot happen with the sinh model that you've shown. The sinh function is unimodal, whereas the plateau in your data has oscillations considerably bigger than 0.001.

Answer 1

[Mathieu NOE]

hello

fyi, some results obtained with the poor's man solution (based on fminsearch)

NB : I changed your log10 to log in my code

for the first set of data , this is what I could obtain

sol = 0.0019 2.0232 0.0002 0.5240

plot in log log scale

x1 = [0.0005 0.0007 0.0008 0.001 0.0013895 0.0019307 0.0026827 0.0037276 0.0051795 0.0071969 0.01 0.013895 0.019307 0.026827 0.037276 0.051795 0.071969 0.1 0.13895 0.19307 0.26827 0.37276 0.51795 0.71969 1 1.3895 1.9307 2.6827 3.7276 5.1795 7.1969 10 13.895 19.307 26.827 37.276 51.795 71.969 100 138.95 193.07 268.27 372.76 517.95 719.69 1000];
y1 = [0.011171 0.015461 0.017605 0.021901 0.029975 0.041604 0.057316 0.078645 0.10726 0.14398 0.18833 0.24482 0.31533 0.38492 0.45433 0.51514 0.53269 0.54151 0.56182 0.56775 0.56635 0.56416 0.55786 0.55369 0.53971 0.52712 0.51333 0.50512 0.49684 0.48869 0.50622 0.51273 0.53057 0.55463 0.61141 0.78976 0.96969 1.1629 1.4075 1.6804 2.1478 2.6287 3.5048 4.4863 5.6525 7.3741];

x = x1(12:46);
y = y1(12:46);

% curve fit using fminsearch
% We would like to fit the function
f = @(a,b,c,d,x) a*sinh(b+(log(x) - c))+ d;
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y);
C1_guess = [0.001 1 0.001 1];
sol = fminsearch(obj_fun, C1_guess); % 
a_sol = sol(1);
b_sol = sol(2);
c_sol = abs(sol(3));
d_sol = sol(4);
y_fit = f(a_sol, b_sol,c_sol, d_sol, x);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

loglog(x, y_fit, '-',x,y, 'r .-', 'MarkerSize', 25)

title(['Fit equation to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Mathieu NOE]

And for the second data set

sol = 0.0142 0.3511 0.0112 0.6948

x2 = [0.0005 0.0007 0.0008 0.001 0.0013895 0.0019307 0.0026827 0.0037276 0.0051795 0.0071969 0.01 0.013895 0.019307 0.026827 0.037276 0.051795 0.071969 0.1 0.13895 0.19307 0.26827 0.37276 0.51795 0.71969 1 1.3895 1.9307 2.6827 3.7276 5.1795 7.1969 10 13.895 19.307 26.827 37.276 51.795 71.969 100 138.95 193.07 268.27 372.76 517.95 719.69];
y2 = [0.0098762 0.013754 0.015562 0.019716 0.027332 0.038068 0.052772 0.073219 0.10152 0.13964 0.19046 0.25666 0.43342 0.53389 0.62881 0.7047 0.74682 0.73801 0.68675 0.64029 0.60505 0.59568 0.59975 0.61275 0.6274 0.64335 0.6483 0.64144 0.6314 0.64113 0.66358 0.70689 0.74901 0.79417 0.86603 0.93868 1.1244 1.3195 2.0893 2.6337 3.0766 3.6119 4.7272 6.0726 7.5609];

x = x2(12:45);
y = y2(12:45);


% curve fit using fminsearch
% We would like to fit the function
f = @(a,b,c,d,x) a*sinh(b+(log(x) - c))+ d;
obj_fun = @(params) norm(f(params(1), params(2), params(3), params(4),x)-y);
C1_guess = [0.001 1 0.001 1];
sol = fminsearch(obj_fun, C1_guess) % 
a_sol = sol(1);
b_sol = sol(2);
c_sol = abs(sol(3));
d_sol = sol(4);
y_fit = f(a_sol, b_sol,c_sol, d_sol, x);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

loglog(x, y_fit, '-',x,y, 'r .-', 'MarkerSize', 25)

title(['Fit equation to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Mathieu NOE]

hello

problem solved ?