lsqnonlin and stretched exponential function

Hi all,

I am currrently trying to fit my data with lsqnonlin instead lsqcurvefit for a comparative reason and I am facing some troubles.

Starting with the easiest example in https://se.mathworks.com/help/optim/ug/lsqnonlin.html, I modify the function as:

fun = @(r) r(1).*exp(-d.*r(2)).^r(3)-y

This function yields x = [1.0376 3.1962 0.4104].

From here, I calculate my errors in the way:

[x,resnorm,residual,exitflag,output,lambda,J] = lsqnonlin(fun,x0);

N = length(y(:,1));

[Q,R] = qr(J,0);

mse = sum(abs(residual).^2)/(size(J,1)-size(J,2));

Rinv = inv(R);

Sigma_var = Rinv*Rinv'*mse;

x_er = full(sqrt(diag(Sigma_var)));

Here, the values I get are not making any sense to me x_er = [0.02701 6458034.8422 829226.8633]; especially for x(2) and x(3).

Fit and errors are totally fine if I fit similar data in the form: (i) r(1).*exp(-d.*r(2)) or (ii) r(1).*exp(-d.*r(2))+r(3).*exp(-d.*r(4));

Could you please help me? Am I missing something?

Thanks a lot.

Best wishes

Alessandro

[Ameer Hamza]

Can you share a sample dataset? I suspect the error is related to the initial guess x0.

[Alessandro]

I used this:

rng default 
d = linspace(0,3);
y = exp(-1.3*d) + 0.05*randn(size(d));

I get the same issues with my experimental data.

I also attach everything else on my test script:

fun = @(r) r(1).*exp(-d.*r(2)) .^r(3)-y;
x0 = [1 4 0.2];
 
[x,resnorm,residual,exitflag,output,lambda,J]  = lsqnonlin(fun,x0)
 
N = length(y(:,1));
 
[Q,R] = qr(J,0);
mse = sum(abs(residual).^2)/(size(J,1)-size(J,2));
Rinv = inv(R);
Sigma_var = Rinv*Rinv'*mse;
x_er = full(sqrt(diag(Sigma_var)))';

figure(1)
hold off
plot(d,y,'ko',d,x(1).*exp(-x(2).*d).^x(3),'b-')

Thanks

A

[Ameer Hamza]

You code seems to give correct output

[Alessandro]

Yes, indeed but I need the errors as well.

When I insert the exponent in the function, i.e., r(3), something goes extremely wrong.

Any idea?

Answer 1

[Torsten]

error_bounds = nlparci(x,residual,'jacobian',J)

after the call to lsqnonlin.

If you don't have a toolbox with nlparci, search the web for nlparci.m.

[Alessandro]

Thanks. I am getting this output:

error_bounds = [ 0.985780584540947 1.08943145408670

-13252964.8453212 13252978.5087471

-372444.576861838 372444.960841471];

But the fitted params are still ~ [1.03760601931382 6.83171295901538 0.191989816441892]

Why is it generating these huge errors? Once again, this appears only when I add the strecth exponent to the function.

[Torsten]

exp(-d*r2)^r3 = exp(-d*r2*r3)

Thus r2 and r3 are dependent fitting parameters which explains the behaviour.

[Alessandro]

same result, nothing changes unfortunately.

fun = @(r) r(1).*exp(-d.*r(2).*r(3)) -y;

[Torsten]

What I meant is that - as written - your fitting problem is ill defined.

Say you get r2=3 and r3 = 5 as result. Then r2 = 1 and r3 = 15 give the same goodness of fit. r2 and r3 are dependent parameters and can be replaced by a single parameter r23.

Maybe you mean exp((-d*r2).^r3) in your model function ?

[Alessandro]

No, I use the fucntion as it is:

y = r(1)*exp(-d*r2)^r3

I cannot replace with one parameter (r23) as you had suggest because I need all three to characterize my samples etc. Why the addition of the r(3) is messing up the errors? This not happens in lsqcurvefit

[Torsten]

Your function to be fitted is

y= r1*exp(-d*r2*r3).

with three fitting parameters.

Say r1 = 1, r2 = 3 and r3 = 4 is what the solver returns.

But why should this result be better than

r1 = 1, r2 = 2 and r3 = 6 because both give the same fitted data.

So only the product r2*r3 is relevant, and the result returned by the solver was only arbitrary.

Mathematically speaking, parameters r1 and r2 do not contain individual information -only their product does.

Test it by choosing different initial values for r2. My guess is that the final r2 and r3 in each run will be different, but their product will remain constant.

If you search for "streched exponential function" in Wikipedia, it's of the form

exp(-r2*d^r3), not

exp(-r2*d)^r3.

[Alessandro]

Thanks. Right but I still don't get why the errors I got are making no sense to me

[Torsten]

There is no such thing as confidence intervals for dependent parameters.