Fit parameter not changed

Hi

I am trying to fit a function to some data points. My code is the following:

clc
clear all
close all

format longE

v = load('01.txt');
M=1600;
N=3000;
v(M:N, 2) = v(M:N, 2)*4.5e6; 

model = @(a, t)(...
    ((-a(1)+sqrt(a(1)*a(1)+4*a(2)*a(3)))/(2*a(2)))*...
    (1-exp(-(a(1)+2*a(2)*((-a(1)+sqrt(a(1)*a(1)+4*a(2)*a(3)))/(2*a(2))))*(t-a(4))))./...
    (1 +  ((a(2)*((-a(1)+sqrt(a(1)*a(1)+4*a(2)*a(3)))/(2*a(2))))/...
    (a(2)*((-a(1)+sqrt(a(1)*a(1)+4*a(2)*a(3)))/(2*a(2)))+a(1)))*...
    exp(-(a(1)+2*a(2)*((-a(1)+sqrt(a(1)*a(1)+4*a(2)*a(3)))/(2*a(2))))*(t-a(4)))));

figure(2)
temp = [100 1e-12 1e9 -0.012];
t=(-0.05:0.01:1)';
plot(t, model(temp, t), 'r', v(M:N,1), v(M:N,2), 'b')
xlim([-0.1 1])
ylim([0 2e7])

b1=temp;

figure(2)
x=-0.00:0.0001:0.3;
plot(x, model(b1, x), 'b--', v(M:N, 1), v(M:N, 2), '--')
ylim([-1 5])

b = lsqcurvefit(model, b1, v(M:N, 1),v(M:N, 2), [], []);

figure(3)
plot(x, model(b, x), '--',v(M:N, 1), v(M:N, 2), '--')

My data is here: http://uploading.com/files/get/449141c1/01.txt. First of all, the first degree of freedom a(1) must be positive (it is a lifetime). My problem is that the third degree of freedom is not varied much by the fitting routine, and I am worried about that as my guess is not optimal.

Am I using the routine as it is supposed to be used?

Best, Niles.

[Star Strider]

I couldn't download you data. That site wanted me to sign up for a ‘plan’ and pay for the privilege. I refused.

I wish TMW had a way to do this, at least temporarily.

If you have a small data set or can post representative (perhaps decimated) data here, that would be preferable. There is no formal protocol for that, so I suggest you used ‘Edit’ to add it at the end of your original post, label it as DATA, label the columns, and format it as ‘Code’.

How well do your data match the fit your function estimates? Do the fitted parameters themselves make sense? If a(1) must be positive, you can constrain it with lsqcurvefit, setting its limit to ‘eps’ or whatever value makes sense. If you don't want to constrain the others, you can set their lower limits at -Inf in your ‘lb’ vector.

The third parameter may not move because you made a good initial guess for it, or it may not be important to the model. If you want more information on the parameters, ask for more output arguments from lsqcurvefit, then use the ‘nlparci’ function. If the 95% confidence limits for the third parameter include zero, then it may not be necessary in the model.

[Niles Martinsen]

Hi

I am sorry about that, I thought a user was able to download it without any fuss. The data is too large to post it here, so I tried another upload-service: http://peecee.dk/upload/view/377162

This should work. If not, please let me know and I'll paste a representative sample of the data points manually. Sorry about the other service.

In my above script I plot obtained function alongside the data, and there is visually a good match. The parameters also make good sense. But the thing is that regardless of what initial value I choose for a(3), it stays the same. That is why I suspect something is wrong.

Answer 1

[Star Strider]

This time the ‘01.txt’ downloaded. You seem to be doing everything correctly, at least in a MATLAB sense.

I was able to run your code with your data. The reason a(3) doesn't change much is that it is likely not uniquely defined in the region of fit, but given other starting values (for example 1E+11) it will converge in the region of 1E+9. Using this change in your code:

[b,resnorm,residual,exitflag,output,lambda,J] = lsqcurvefit(model, b1, v(M:N, 1),v(M:N, 2), [], []);
ci = nlparci(b,residual,'jacobian',J)

The 95% confidence intervals only exist for a(4) and do not include zero. They are NaN for all the others. So your system may be over-determined in that you may need fewer than four parameters to fit it.

I got a good fit with these starting parameters:

100 1e-12 2.5E+9 -0.012

estimating these resulting parameters:

410.0000e+000   -16.0906e-006     2.5000e+009   -12.4858e-003

and this starting set:

100 1e-12 0.9E+9 -0.012

resulted in these estimates:

22.4418e+000     6.6366e-006   900.0000e+006   -14.2420e-003

with the same visual accuracy.

The function fit is sensitive to a(3), it is just that a guess on the order of 10^9 for it is close enough to the best fit for it that it doesn't have to change much. (Vary it to experiment with it.) I suspect you will get a range of parameter values that will fit your function, as I did.

[Niles Martinsen]

Hi Star Strider

Thanks for that, I see what you mean. Maybe it gets easier if I can constaint some of the parameters during the fit. Do you know how I can do that?

Best, Niles.

[Star Strider]

The lsqcurvefit function will allow you to set and lower and bounds on the parameters. The documentation for it refers to them as ‘lb’ and ‘ub’ respectively. Those are the two sets of empty brackets in the code you quoted above:

b = lsqcurvefit(model, b1, v(M:N, 1),v(M:N, 2), [], []);

All you have to do now is to put appropriate values for your parameter limits in them (in the same order and row vector or column vector orientation as your parameter vector). You can of course set these as variables and refer to them as such in the call to lsqnonlin.