Nonlinear fit comparison - Matlab vs OriginLab

Good afternoon, I'm posting this question because I'm not able to find a satisfactory result of a non-linear fit with Matlab. On the contrary, OriginLab returns a good fit together with errors for fit parameters. Here my x and y data: x = [177600,961200, 2504000, 4997000, 8884000]; y = [6.754, 24.416, 58.622, 107.980, 154.507]; yErr = sqrt(y);

The fitting function is the following: y = V + (A-V)*S*x*EXP(-S*x) where V, A and S are the parameters which I need to calculate from the fit. OriginLab gives me the following:

• V = 1.978 ± 0.889
• A = 585 ± 64
• S = 4.392E-8 ± 6.771E-9

On Matlab I tried all the possible form of non-linear regression: fitnlm, lsqcurvefit, nlinfit etc... every time the same warning: The Jacobian is ill conditioned. Here is an example:

    Nonlinear regression model:
    y ~ a1 + (a2 - a1)*a3*x*exp( - a3*x)

Estimated Coefficients:
         Estimate         SE           tStat         pValue  
        __________    __________    ___________    __________
  a1       -0.6789    2.9104e-19    -2.3327e+18    2.0265e-73
  a2        530.78    8.0894e-20     6.5614e+21    3.2371e-87
  a3    5.2511e-08    5.1865e-10         101.25    5.7063e-08

Of course such small standard errors are not reliable, but still I'm not able to reproduce the Origin results. Any suggestions?

[Star Strider]

Using the genetic algorithm ga function (Global Optimization Toolbox), the best fit I got was with:

V  =   7.550390625000001E-01
A  =   5.476210479736328E+02
S  =   4.920000007402336E-08
RN =   3.197638945910025E+00

where ‘RN’ is the norm of the residuals. (I constrained it to have only positive parameter estimates.)

[Stephan]

I had a result which was near to this from Star Strider using fit using this code:

x = [177600,961200, 2504000, 4997000, 8884000]';
y = [6.754, 24.416, 58.622, 107.980, 154.507]';
x0 = [585 0 2]; 
fitfun = fittype('V + (A-V) *S .*x .* exp(-S .* x)');
[fitted_curve,gof,output] = fit(x,y,fitfun,'Start',x0)
scatter(x, y, 'r+')
hold on
plot(x,fitted_curve(x))
hold off

Results are:

fitted_curve = 

       General model:
       fitted_curve(x) = V + (A-V) *S .*x .* exp(-S .* x)
       Coefficients (with 95% confidence bounds):
         A =         555  (423.4, 686.6)
         S =    4.82e-08  (2.599e-08, 7.04e-08)
         V =      0.8133  (-8.284, 9.91)

gof = 

    struct with fields:

             sse: 10.5406
         rsquare: 0.9993
             dfe: 2
      adjrsquare: 0.9986
            rmse: 2.2957

output = 

    struct with fields:

             numobs: 5
           numparam: 3
          residuals: [5×1 double]
           Jacobian: [5×3 double]
           exitflag: 2
      firstorderopt: 1.3449e+08
         iterations: 4
          funcCount: 20
       cgiterations: 0
          algorithm: 'trust-region-reflective'
           stepsize: 1.8486e-10
            message: 'Success, but fitting stopped because change in coefficients less than tolerance (TolX).'

[MZ123]

@Star Strider what about the standard errors associated with the parameters?

[MZ123]

@Stephan almost close to Origin, but V has an enormous confidence interval.. :( thank you

Answer 1

[Matt J]

Well, I guess I would credit OriginLab with knowing how to pre-scale the x,y data appropriately. If you do this manually, then lsqcurvefit works fine, and gives a slightly better fit (according to resnorm) than OriginLab,

c=10*max(x); %pre-scaling factor

a0=[1.9780 ,585, 4.3920e-08*c];
F=@(a,x) a(1)+(a(2)-a(1)).*a(3).*x/c.*exp(-a(3).*x/c);

[a,resnorm,~,~,~,~,J]=lsqcurvefit(F,[1,100,1],x,y,[],[],...
                       optimoptions('lsqcurvefit','Display','iter'));

solution=[a(1), a(2), a(3)/c]

resnorm %10.1985
condJ = cond(full(J))  %510.8393

plot(x,y,'x',x,F(a,x),'--'); xlabel 'x', xlabel 'y';

[MZ123]

I didn't get the point. I noticed you just applied a scaling factor, but how to get the Standard Error with the parameters? Also, how do I re-map the scaled coefficient to those similar to Origin? Thank you for your help

[Matt J]

I didn't get the point.

The point is that we reached a better fit than OriginLab - the residuals are lower. Also, there are no ill-conditioned Jacobian problems once you scale the data properly. With c=1 (no scaling) you will find that cond(J)=Inf whereas with scaling, it drops nicely down to 60 or so.

Also, how do I re-map the scaled coefficient to those similar to Origin?

I already remapped them for you in the line

solution=[a(1), a(2), a(3)/c]

how to get the Standard Error with the parameters?

You have the Jacobian (a sensible one now) and the residuals, so you can use them to compute whatever gof metrics that you like.

[Matt J]

You have the Jacobian (a sensible one now) and the residuals, so you can use them to compute whatever gof metrics that you like.

which leads me to

• V = 0.6652 ± 1.7676
• A = 545 ± 33
• S = 4.9629e-08 ± 5.1273e-09

Answer 2

[MZ123]

Another frustrating scenario.. simple polynomial fit y = a+bx+cx^2 I get another error on Matlab "Polynomial is badly conditioned". On the other hand Origin does it better.. why? To me this is really embarrassing.

[Stephan]

Could you provide your code? Yesterday in the evening i tried to understand the problem you have. I did not get the exactly same result you provided from OriginLab, but it was a good fit i thinkand there was no warning. So it would be interesting to see what you are doing.

One more question to your fitting equation:

In your question it is:

y = V + (AV) * S * x * EXP (-S * x)

which should mean

A * V

In the Second Part of your question it is:

y ~ a1 + (a2 - a1) * a3 * x * exp (- a3 * x)

which is

A - V

[MZ123]

Thank you everybody for your answers. For Stephan: in my function I have (A - V) and also the exponential part. In my second part I just put a very simple parabola-like polynomial fit to emphasize why it is (to me) so difficult to obtain in Matlab the same results which I got with Origin.

For what concerns the code I'm running... I would say that I'm applying straightforwardly what I learned in the documentation. For example in the parabola case what I run is this:

x = [177600,961200, 2504000, 4997000, 8884000];
y = [6.754, 24.416, 58.622, 107.980, 154.507];
beta0 = [1.5,1E-5,-1E-12];
myfun = 'y~b1+b2*x+b3*x^2';
mdl1 = fitnlm(x,y,myfun,beta0)

Estimated Coefficients:
         Estimate          SE          tStat         pValue  
        ___________    __________    __________    __________
  b1         1.0237    2.3928e-13    4.2781e+12    2.8166e-38
  b2     2.6025e-05    4.9385e-07        52.699    1.5049e-05
  b3    -9.8098e-13    6.2704e-14       -15.645    0.00056757

Warning: The Jacobian at the solution is ill-conditioned, and some model parameters may not be estimated well (they are not identifiable). Use caution in making predictions.

Of course parameter b1 has a acceptable estimate value, but the SE is not accurate at all. Origin retrieves a SE = 0.5 more or less..

Answer 3

[Alex Sha]

The global solution looks like below:

Root of Mean Square Error (RMSE): 1.42818042380826

Sum of Squared Residual: 10.1984966147457

Correlation Coef. (R): 0.999655170571204

R-Square: 0.999310460049744

Adjusted R-Square: 0.998620920099487

Determination Coef. (DC): 0.999310460049744

Chi-Square: 0.212210508410738

F-Statistic: 1449.24229503583

Parameter Best Estimate

---------- -------------

v 0.665211685653936

a 544.888484613983

s 4.9628897712431E-8