MATLAB: Non Linear Fit help

non linear fit

I have been trying to fit some data with a model for some time now and am getting fairly close. In origin pro the fit seems to work ok and I get reasonable parameters back albeit the residuals are a little funny. Unfortunately in Matlab the fit is not so good. Could I send some one my code my function and some data to take a look at?

Oh and also there are imaginary numbers!

    clear all;
    close all
    close all hidden
    A = uiimport;
     clear xdata;
     clear ydata;
     xdata = A.data(:,1);
     ydata = A.data(:,2);
     Vg = xdata;
     Id = abs(ydata);
    x0 = [8.6E-10;1.7;0.8;5E6];
    options = optimset('Display','iter',...
             'TolFun',1E-100,'TolX',1E-100,'MaxIter',1000);
    [beta,resid,J,COVB,mse] = nlinfit(Vg,Id,@RcFun,[x0],options);
    [ci se] = nlparci(real(beta),resid,'covar',COVB);
    Id_new = ((((real(beta(1))*((Vg-real(beta(2))).^(real(beta(3))+1))).^(-1))+real(beta(4))).^(-1));
    plot(Vg,Id_new,'r',Vg,Id,'o');
    hold on;
    plot(Vg,resid,'+');
 function F = Rcfun(a,Vg)
  K = real(a(1));
  V = real(a(2));
  b = real(a(3));
  R = real(a(4));
F = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1);
 end

Data

A  B
0  3.03E-12
1  1.5E-13
2  1.58E-12
3  2.81E-12
4  2.55E-12
5  2.31E-12
6  4.13E-12
7  2.89E-12
8  4.99E-12
9  6.38E-12
10  1.068E-11
11  1.96E-11
12  5.343E-11
13  5.405E-11
14  5.347E-11
15  5.142E-11
16  2.4139E-10
17  7.4428E-10
18  1.5752E-9
19  2.7328E-9
20  4.3347E-9
21  6.5506E-9
22  9.5258E-9
23  1.31356E-8
24  1.72672E-8
25  2.17876E-8
26  2.66302E-8
27  3.18252E-8
28  3.7101E-8
29  4.23594E-8
30  4.78078E-8
31  5.32604E-8
32  5.86136E-8
33  6.39262E-8
34  6.93234E-8
35  7.47466E-8
36  8.01152E-8
37  8.54398E-8
38  9.08214E-8
39  9.62598E-8
40  1.0184E-7
41  1.074E-7
42  1.1322E-7
43  1.1876E-7
44  1.2432E-7
45  1.299E-7
46  1.3534E-7
47  1.4062E-7
48  1.4596E-7
49  1.5096E-7
50  1.558E-7
51  1.6118E-7
52  1.6616E-7
53  1.7064E-7
54  1.7546E-7
55  1.7946E-7
56  1.8402E-7
57  1.8776E-7
58  1.9138E-7
59  1.9584E-7
60  1.9992E-7

Best Answer

Thank you for clarifying ‘RcFun’. After working with your code for a bit, the problem definitely seems to be the ‘(Vg-V)’ term. The only way I am able to get a decent fit without complex parameters is to use ‘lsqcurvefit’ and constraining ‘V’ to not be greater than zero, but then ‘nlparci’ wouldn't calculate confidence intervals on the parameter estimates. [I was able to get a good fit with ‘nlinfit’ by redefining that term to ‘abs(Vg-V)’, but that changed your function.] With ‘lsqcurvefit’, your ‘options’ statement and structure remains the same, although I increased ‘MaxIter’ and ‘MaxFunEvals’ in mine for diagnostic purposes. I also made ‘RcFun’ a single-line anonymous function for convenience:

RcFun = @(a,Vg) 1./(1./(a(1).*(Vg-a(2)).^(a(3)+1)) + a(4));
Lobnd = [];
Upbnd = ones(4,1)*realmax;
Upbnd(2) = 0;        % Constrain ‘V’ to not be greater than zero
[beta,resnrm,resid,xitflg,outpt,lmda,J] = lsqcurvefit(RcFun,x0,Vg,Id,Lobnd,Upbnd,options);  % NOTE: ‘Answers’ wrapped this line

and got an acceptable fit with these parameters:

 K =     35.6744e-012
 V =    -32.7518e-003
 b =      1.1302e+000
 R =    145.7789e-003

The norm of the residuals was 5.2916E-015. I'm not certain what you're doing or what data you're modeling, so I can't interpret these or their validity.

In spite of the decent fit, the residuals had a definite oscillating trend.

Related Solutions

MATLAB: Lsqcurvefit problem. Error using snls (line 48) Objective function is returning undefined values at initial point. lsqcurvefit cannot continue.

There are multiple issues here to deal with.

First of all, why is the objective returning undefined values? Probably because you have a NaN in your data.

find(isnan(ydata))
ans =
   201
find(isnan(xdata))
ans =
  1×0 empty double row vector

So ydata(201) is a NaN.

Next I'll check your start point, and the numerics in general.

x0=[1.17e-4 0.00272 114.31287];
plot(xdata,ydata,'o')
grid on

You are using this model

y = a(1) + a(2)*exp(-x*/a(3))

which at a glance, seems fine.

So lets see how good your start values are, just by use of some pencil and paper estimates.

For that model, the asymptote at x==inf is a(1). Looking at the data in that plot, I'd guess 1.5e-4 for the asymptote. So 1.17e-4 is not bad, within 50% of my guess. Since the curve has essentially bottomed out over the last 25 points or so, this should give me a better guess:

nanmean(ydata(175:end))
ans =
   0.00015741

That seems to agree with my graphical estimate.

As far as the exponential rate parameter is concerned, a quick graphical estimate would have me guess a(3) should be roughly 100. So, it looks like you have chosen reasonably well.

Next, you have one NaN value in your data. The last point was a NaN in ydata. That won't really hurt anything as long as the tool you use for the fit excludes it. The problem is that lsqcurvefit does not seem to exclude NaNs in the data.

Next, I'll use my own fminspleas tool (found on the file exchange) to do the fit.

mdlterms = {1,@(C,x) exp(-x./C)};
[nlt,lint] = fminspleas(mdl,10,xdata(1:end-1),ydata(1:end-1))
nlt =
       119.06
lint =
   0.00019825
    0.0088797

So it estimates the model as:

1.9825e-4 + 8.8797e-3*exp(-x/119.06)

That seems in line with the estimates I found myself, and within shooting range of what you used as starting values.

plot(xdata,ydata,'o')
grid on
hold on
yhat = 1.9825e-4 + 8.8797e-3*exp(-xdata/119.06);
plot(xdata,yhat,'r-')

which, as curve fits go, is not too bad. I'd guess there is a little lack of fit in the tail. Along with that, it seems like there is a larger amount of noise at the bottom end of the curve. It seems like non-uniformity in the error, so heteroscadisticity, but that could be just the high slope of the curve on the left end is quite high, so we miss seeing any noise one that end.

A regression spline fit will resolve that question. So my SLM toolbox:

mdl = slmengine(xdata,ydata,'plot','on','decreasing','on','concaveup','on');

Then the residuals:

So here we see the residuals for a good fit are roughly the same size along the entire curve. The high slope just hides the noise on that end.

Otherwise, what you did should work acceptably.

fun=@(x,I) x(1)+x(2).*exp(-I/x(3));
ub=[1 100 1000];
lb=[1e-5 1e-5 10];
x0=[1.17e-4 0.00272 114.31287];
options=optimset('MaxIter',2000,'TolFun',1.e-5,'TolX',1.e-5,'Algorithm', ...
'trust-region-reflective',... 
'MaxFunEvals',2000,'Display','iter');
k = ~isnan(ydata) & ~isnan(xdata);
x = lsqcurvefit(fun,x0,xdata(k),ydata(k),lb,ub,options);

It now works fine.

                                           Norm of      First-order 
   Iteration  Func-count     f(x)          step          optimality
       0          4     0.000582026                          8.52
       1          8     1.55855e-06       0.079033       1.63e-05      
       2         12     1.29161e-06        0.10269       3.47e-06      
Local minimum found.
Optimization completed because the size of the gradient is less than
the selected value of the optimality tolerance.

What result came out?

x
x =
   0.00020904    0.0089523       117.73

So quite reasonable. As I said, maybe a little lack of fit in the model at worst. I'd still use fminspleas for the fit here, because it will be more accurate, converge faster, and be more robust. For example, fminspleas found a better solution than did lsqcurvefit.

nanstd(ydata - (lint(1) + lint(2)*exp(-xdata/nlt)))
ans =
   8.0291e-05
nanstd(ydata - (x(1) + x(2)*exp(-xdata/x(3))))
ans =
   8.0552e-05

but not a huge difference.

MATLAB: Curve fitting of two equations

I'll expand on this answer in the morning. But, yes, you can use lsqcurvefit here, or lsqnonlin.

The y equation in particular suggests logging it however. Anytime that you see large differences in scale across the range of your problem, it makes sense that taking a log will be beneficial. Essentially this presumes proportional error in the fit is appropriate.

If you don't log the y expression, then the problem is, lsqnonlin presumes additive, homoscedastic error. So it assumes that the point with y == 0.01 will have the same absolute importance as the point where y=5.445.

As well, the exponentials in there also suggest taking logs.

x = [36.78 ,37.53 ,38.28 ,39.06 ,39.85 ,40.65 ,41.47 ,42.3162082317748,43.17 ,44.04 ,44.93 ,45.84 ,46.76 ,47.71 ,48.67 ,49.65 ]
y = [0.01 ,0.0152 ,0.023 ,0.035,0.0536 ,0.0816 ,0.12 ,0.1891 ,0.287 ,0.438 ,0.6668 ,1.01494 ,1.544 ,2.35097424365239,3.578 ,5.445 ]
f=0:0.2:3;

So, you have TWO equations that will allow you to predict x and y, as a function of f, and the parameters teta and a.

y=teta*exp((a+1)*f);
x=exp((a-1)*f-1)/(a-1)^2;

Logging the y equation gives us:

log(y) = log(teta) + (a+1)*f;

If we had no more than tis equation, then we could compute teta and a directly using simple linear regression.

The x equation need not be logged, because all the values of x are essentially the same size. So additive versus proportional error is not a real factor. I'll leave it alone, although there is no problem either way.

So, now lets formulate an objective function. I'll use lsqnonlin here,as I prefer it for this specific case. It is not a strong preference though.

I'll implicitly put teta into params(1), a as params(2).

fun = @(params) [x - exp((params(2)-1)*f-1)/(params(2)-1)^2, log(y) - log(params(1)) - (params(2)+1)*f];

Now, just call lsqnonlin on fun.

You probably want to put lower bounds on teta as 0, as you don't want it trying to try out teta less than zero. In fact, I'd suggest setting the lower bound for teta as something like 1e-10, to prevent it from even trying to go negative. If it does, lsqnonlin will fail, due to the complex numbers thereby created.

Similarly, set the lower bound on a as something like -1+1e-10.

The essential idea here is to pack the x and log(y) equations together into one set, then using lsqnonlin to solve the entire set for the combined paramters teta and a.

Good starting values, just be looking at the equations MIGHT be easily gained, from the log(y) equation. Simply fit a linear polynomial to that curve.

p0 = polyfit(f,log(y),1)
p0 =
      2.102      -4.6117

That implies

a+1 = 2.102

and

log(teta) = -4.6117

solving for a and teta, we get

teta0 = exp(p0(2))
teta0 =
    0.0099353
a0 = p0(1) - 1
a0 =
        1.102

These are logical numbers. we know that a MUST be slightly smaller or larger than 1. As a test, I'll solve the problem twice, starting a from slightly above or below 1.

params = lsqnonlin(fun,[teta0,a0],[1e-10,-1+1e-10])
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the function tolerance.
<stopping criteria details>
params =
    0.0099644          1.1
>> 
>> params = lsqnonlin(fun,[teta0,0.99],[1e-10,-1+1e-10])
Local minimum possible.
lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.
<stopping criteria details>
params =
     0.013194      0.91282

Best Answer

Related Solutions

MATLAB: Lsqcurvefit problem. Error using snls (line 48) Objective function is returning undefined values at initial point. lsqcurvefit cannot continue.

MATLAB: Curve fitting of two equations

Related Question