bad fit with lsqnonlin LM

I am trying to create a solver that will fit the following. Y data has been normalised to 1.

param0 = 0;
opts = optimset('lsqnonlin'); % method
opts.Algorithm = 'levenberg-marquardt';
opts.Display = 'iter';

TR=0.0055

Xk = [0.0785398163397448 0.165806278939461 0.235619449019234 0.270526034059121 0.523598775598299];
Yk = [0.451237263464338 0.810771470160116 0.965065502183406 1 0.991266375545852]

modelfun = @(a) ((1-exp(-TR/a))*sin(Xk)) ./ (1-cos(Xk)*exp(-TR/a)) - Yk;
[estcoefs,resnorm,residual,exitflag,output] = lsqnonlin(modelfun,param0,[],[],opts);

I am getting unrealistic values from using lsqnonlin. For this data, I know the result should be 0.08. This can be verified if I put 0.08 into the model and compare with the data.

plot(Xk,Yk,'b')
plot(Xk,((modelfun(estcoefs))+Yk)/max(modelfun(estcoefs)+Yk),'r');
plot(Xk,(modelfun(0.08)+Yk)/max(modelfun(0.08)+Yk),'g')

Can anyone help me with this issue?

Answer 1

[Matt J]

Your criterion for a good fit is quite strange. Regardless, lsqnonlin uses a least squares fitting criterion and, with respect to that criterion, lsqnonlin has succeeded in giving you a superior solution:

>> norm(modelfun(estcoefs))

ans =

      1.3673

>> norm(modelfun(.08))

ans =

      1.6047

[Matt J]

I'm guessing that what you were really after the following,

modelfun = @(a) a(1)*((1-exp(-TR/a(2)))*sin(Xk)) ./ (1-cos(Xk)*exp(-TR/a(2))) - Yk;

param0 = [5,0.08];
[estcoefs,resnorm,residual,exitflag,output] = lsqnonlin(modelfun,param0,[0,0],[],opts);

estcoefs,
hold on
plot(Xk,Yk/max(Yk),'b--o')
plot(Xk,modelfun(estcoefs)+Yk,'r-*');
plot(Xk,modelfun(param0)+Yk,'g-+')
hold off

[Erin Wang]

Thanks Mark! Yes indeed I did have the a(1) in there initially, but got fed up and removed it, thinking it wouldn't make a difference.

Now I put it back, and realised that it is the initial parameters that was giving the problem. If I gave a(2) an initial of 0, it would just give me 0 back. Anything that is not 0 would give me the same results ~0.083 which is what I expect. Checking the norm(modelfun(estcoefs)) found that with a(2)=0 it is around 0.7, and a(2)~0.83 it is 0.02... Excuse my naivety with the fitting algorithms and optimisations, I didn't know that it would matter this much! Is this something that comes up a lot? I would expect regardless the initial parameters, it would try give a superior fit. But in this case it didn't..

[Matt J]

I would expect regardless the initial parameters, it would try give a superior fit.

It always tries, but can fail, depending on the quality of your initial guess.

In any case, a(2)=0 is not a legal guess, because it will cause divide-by-zero operations TR/a in your modelfun().