Not quite fitting the data using lsqcurvefit

Hello everyone. I'm trying to use lsqcurvefit to get optimal parameter value for K (see code). The fitting looks fine except for 2 things: fitted curves are always shifted down by some amount to respect to the data (see picture), and the values I get for K are complex when it should be a real number (but I guess that's a subproduct of the main problem here). K values for different initial guesses are similar, which is good. I really don't know what is going on, I hope someone can give me a hint.

Here's the code:

R=[0.1:0.1:0.4 0.6:0.2:1.8 2.1 2.4:0.4:3.2 3.8 4.8 6 8.4 12.8 20 36];
pks_locs1 = [114 93.5 144 167 199.5 225.5 275 271.5 282.5 301.5 315.5 319.5 348 356.5 362.5 390.5 408.5 420.5 464.5 449.5 457.5 477 494.5]';
CIS_H2 = 352.6;
H=0.0002529;
G=H./R';

fun_H2=@(K,R) pks_locs1(1)+CIS_H2*(K*G*(1+R)+1-sqrt((K*G*(1+R)+1)^2-R'*(2*K*G').^2))/(2*K*G');
K0_H2=100;
K=lsqcurvefit(fun_H2,K0_H2,R,pks_locs1(2:end))

plot([0 R],pks_locs1,'ko',R,fun_H2(K,R),'b-')
legend('Data','Fitted exponential')
title('Data and Fitted Curve')

The value I get for K is 9.3155e+03 - 1.1463e-01i. I have tried using lsqnonlin with similar results: fitted curve down-shifted and complex K values.

Thanks in advance!

[Walter Roberson]

R needs to be the same size as pks_locs1(2:end) but at the moment pks_locs1(2:end) is empty.

[Javier Agustin Romero]

Thanks for replying. I put pks_locs1 at the end. Here it is:

pks_locs1 = [114 93.5 144 167 199.5 225.5 275 271.5 282.5 301.5 315.5 319.5 348 356.5 362.5 390.5 408.5 420.5 464.5 449.5 457.5 477 494.5]';

Also, edited the question so pks_locs1 is there, sorry for that.

[Walter Roberson]

Your fun_H2 is not vectorized. You have ^ and you have / and you have some * that appear to be acting as algebraic matrix multiplication.

[Javier Agustin Romero]

What do you mean? give K any scalar value and fun_H2 is a 22x1 vector.

[Walter Roberson]

You have

fun_H2=@(K,R) pks_locs1(1)+CIS_H2*(K*G*(1+R)+1-sqrt((K*G*(1+R)+1)^2-R'*(2*K*G').^2))/(2*K*G');

Notice this includes sqrt((K*G*(1+R)+1)^2) . But your R is a vector, so the ^2 is being applied to a vector, unless the algebraic matrix multiplication by G collapses the vector (1+R) into a scalar. ^ is matrix exponential, not element-by-element exponential. * is algebraic matrix multiplication, not element-by-element multiplication. And further down in the expression you have /(2*K*G') where G is a vector, so the / is matrix division, not element-by-element division.

You should be using .* and .^ and ./ everywhere unless you deliberately want the values at one location to influence the calculation of values at all locations. The / operator is essentially a fitting operation rather than a division: if you want fitting to be taking place there then you should comment heavily .

[Matt J]

@Javier,

Incidentally, also, your code labels the fit as a "Fitted Exponential", but in fact there are no exponentials in your model function, so one wonders whether the model function as you have it now is missing some exponential terms.

Answer 1

[Matt J]

When your model function is fully vectorized, as suggested by Walter, the results are better, but only you can know for sure what your model function is supposed to be.

R=[0.1:0.1:0.4 0.6:0.2:1.8 2.1 2.4:0.4:3.2 3.8 4.8 6 8.4 12.8 20 36].';
pks_locs1 = [114 93.5 144 167 199.5 225.5 275 271.5 282.5 301.5 315.5 319.5 348 356.5 362.5 390.5 408.5 420.5 464.5 449.5 457.5 477 494.5]';
CIS_H2 = 352.6;
H=0.0002529;
G=H./R;
fun_H2= @(K,R)pks_locs1(1)+CIS_H2.*(K.*G.*(1+R)+1-sqrt((K.*G.*(1+R)+1).^2-R.*(2.*K.*G).^2))./(2.*K.*G);


K0_H2=5.0758e+04;
[K,fval]=lsqcurvefit(fun_H2,K0_H2,R,pks_locs1(2:end))
plot([0 R.'],pks_locs1,'ko',R,fun_H2(K,R),'b-')
legend('Data','Fitted exponential','location','southeast')
title('Data and Fitted Curve')

[Javier Agustin Romero]

Well yes, you are both right. I do not want values at one location to influence the calculation of values at all locations, that was my mistake and should use '.' with every operator. lsqcurvefit seems to be working fine but I am not entirely sure about the mathematical model now, it comes from a publication and it should fit the data. I'll keep working on it. Thank you very much for your comments.

[Javier Agustin Romero]

It turns out G was supposed to be a constant, not = H/R... I am working with chemists and my knowledge in chimestry is null. Making G = 0.0002529 (the value I had mistakenly given to H) the fit works fine. Thanks again.