Non Linear Fit help?

Thanks for all the help so far, its kinda making more sense to me but i cant seem to get the way you have done it to converge? What initial guess are you using? How many iterations are you using?

Also the R that has come out of your fit seems quite low despite all the other parameters fitting, do you mean e+3?

If you like to know the formula I am using for this fit is

Id = Vd/((1/(W/L*C*Mu*(Vg-V)))+R)

In the data I gave you Vd is 1. W, L, C are all constants (channel width, channel length and capacitance of a transistor)

Mu is the mobility of the holes in the device and i am assuming that it obeys a power law Mu = K(Vg-V)^b and is Vg(gate voltage dependant)

R is the contact resistance of the device.

V is known as the threshold voltage and is where the dcevice switches on. I chose to change change the values of gate voltage to the opposite sign, but I'm not sure why I did that. if you want the actual raw data change column A to negative numbers.

Thanks for the help

I meant to say Vd was -1, maybe thats why I put everything negative

My pleasure. The initial guess I used for that attempt was:

x0 = rand(4,1);

Since I don't have ‘x0’ display to the Command Window, I don't know what they actually were. I went back and checked to be sure about ‘R’; the parameters I quoted were copy-pasted from the Command Window.

Changing ‘Vg’ to ‘-Vg’ could make a significant difference in the success of the fit.

When I just now changed ‘Vg’ to ‘-Vg’, and changed the initial parameter guess to:

x0 = rand(4,1).*[1E-10;  1;  1;  1E+6];

the estimated parameters were:

K =    80.5824e-012
V =   326.4389e-003
b =   549.8791e-003
R =   388.7843e+003

(That fit did not produce a plot because it threw an error on complex arguments to ‘nlparci’ before it got to the plot.) The problem is that now the residuals are complex! The norm of the residuals is 564.3979E-015, and as might be expected from complex residuals, ‘J’ is complex. (It gave no errors other than its inability to calculate confidence intervals on complex parameters, which may mean the complex Jacobian.)

I also gave it a go with ‘nlinfit’ now that I no longer have to constrain the parameters, and got complex coefficients (with close-to-equal real and imaginary parts) and complex results for ‘Id_new’. I have no experience with transistor physics, so I cannot suggest any changes in your model.

I'll be glad to continue to help.

I think we'll manage to solve this eventually!

For the '-Vg' fit did you also use Vd = -1, as in made the function -ve?

I think the problem is that we are not allowin V to be positive. Origin for the same data gives a good fit with

k = 8.6E-10 V = 17.3 b = 0.618 R = 2.3E6

this is with the function still posititve as when it is negative it doesn't seem to work at all.

In fact using lsqcurvefit makes the fit woese than the original values

When I plugged those parameters into the function I've been using and removed the negative from ‘Vg’ and ‘RcFun’, it behaved similarly to when I used ‘abs(Vg-V)’. The values MATLAB gave using this ‘RcFun’:

RcFun = @(a,Vg) (1./(1./(a(1).*(Vg-a(2)).^(a(3)+1)) + a(4)));

for the range 0V to 20V was:

         Vg                            Id_new          
     0.0000e+000                    41.0816e-009 - 68.1918e-009i
     1.0000e+000                    36.7197e-009 - 63.0046e-009i
     2.0000e+000                    32.5968e-009 - 57.8053e-009i
     3.0000e+000                    28.7182e-009 - 52.6210e-009i
     4.0000e+000                    25.0870e-009 - 47.4801e-009i
     5.0000e+000                    21.7052e-009 - 42.4118e-009i
     6.0000e+000                    18.5731e-009 - 37.4463e-009i
     7.0000e+000                    15.6895e-009 - 32.6150e-009i
     8.0000e+000                    13.0521e-009 - 27.9500e-009i
     9.0000e+000                    10.6578e-009 - 23.4849e-009i
    10.0000e+000                     8.5028e-009 - 19.2546e-009i
    11.0000e+000                     6.5831e-009 - 15.2962e-009i
    12.0000e+000                     4.8952e-009 - 11.6496e-009i
    13.0000e+000                     3.4370e-009 -  8.3593e-009i
    14.0000e+000                     2.2090e-009 -  5.4771e-009i
    15.0000e+000                     1.2174e-009 -  3.0677e-009i
    16.0000e+000                   479.2402e-012 -  1.2228e-009i
    17.0000e+000                    44.4394e-012 -114.2460e-012i
    18.0000e+000                   482.3752e-012                
    19.0000e+000                     2.0200e-009                
    20.0000e+000                     4.2480e-009                

It's possible that the model isn't appropriate for ‘Vg’ < ‘V’.

There nevertheless must be some significant difference between the Origin and MATLAB nonlinear regression functions, because even when I restrict ‘Vg’ >= 18V I'm not getting as good a fit with either ‘nlinfit’ or ‘lsqcurvefit’ as you are with Origin. (I find this extremely disconcerting!) I'm still getting complex parameters for ‘K’, ‘V’, and ‘b’. Those latest (with ‘nlinfit’) are:

beta_est =

     3.3672e-009 - 11.3930e-012i
    19.2532e+000 -682.8867e-015i
   141.3304e-003 +987.7662e-006i
   554.7623e+003                

absbeta_est =

     3.3672e-009
    19.2532e+000
   141.3339e-003
   554.7623e+003

I also get different results from ‘nlinfit’ (better) and ‘lsqcurvefit’. (I'm going to look into Origin. I'd like to know why it comes up with better parameter estimates than MATLAB does.) I wonder if it would be worthwhile to bring this up with TMW Tech Support?

I have the Symbolic Math Toolbox, so I'm going to see if calculating the Jacobian with it and giving that analytically to the functions could solve the problem. That may take a day or so, assuming the Symbolic Math Toolbox can do that. (It too has been having problems of late in some situations.)

Sorry also for the delay. I ran a number of test fits while I was writing this comment.

It does seem odd. Origin seems to be able to get the fit in only aboput 18 iterations or so as well.

you are right about Vg < V. In origin it cuts off the fit there and in in matlab the curve with parameters from origin give an upward curve at values of Vg < V. Origin use Levenberg-Marquardt algorithm for there fits.

Is there someway of stopping the fit from trying to go below V before it has calculated it?

I have been looking into lm.m and LMFsolve.m to see if they could find a better fit but I am strugling a little in getting them to work.

I apprectiate all the help, I hope you are getting something out of this as well.

The vectorized Jacobian (Symbolic Math Toolbox) of:

F = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1) 
	w.r.t:             K,     V,     b,     R 
        respectively       a(1),  a(2),  a(3),  a(4)

is:

J = [1./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)) - (a(1).*a(4))./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2; 

 -(a(1).*(a(3) + 1))./((Vg - a(2)).^(a(3) + 2).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2); 

 (a(1).*log(Vg - a(2)))./((Vg - a(2)).^(a(3) + 1).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2); 

 -a(1).^2./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2]; 

MATLAB wants it as a row vector, I listed it as a column vector for readability. I'll work on it in a few minutes. Add to ‘optimset’ the argument pair 'Jacobian', 'on' to use it, add it to ‘RcFun’, and change ‘RcFun’ to

function [F, J] = RcFun(a,Vg)  

I'm posting it now so you can work on it in parallel.

I'm not at all happy that MATLAB can't produce as good a fit to your function as Origin obviously can, even when ‘Vg’ > ‘V’. If an analytic Jacobian doesn't work to fit your function with precision similar to what Origin is producing, I'm going to contact TMW Tech Support to see if they have a solution. At the very least they need to know about the problem.

I haven't tried lm.m and LMFsolve.m. I wasn't aware of them. I've been pursuing the Jacobian hypothesis, so our doing those in parallel seems an effective strategy.

I'm taking a short break, then will be back to see if the calculated Jacobian can improve the fit. That's the only possibility I can think of that we haven't collectively already tried.

It's my pleasure to help! For my part, I'm learning about organic FETS (new and interesting) and more about curve fitting and parameter estimation in MATLAB. I probably learn more than I contribute.

To be honest, I'm not really sure what you meant above with the jacobian but I have done the above and for some reason am getting a better fit for data with Vg > V but it is still different from Origin and R (a(4)) does not change from the initial guess???

I'm going to give it a break for a bit and head home.

Well ... unfortunately, the analytic Jacobian with ‘nlinfit’ didn't improve anything. I'm still getting complex parameter estimates and a complex fit, although the imaginary parts of the parameters are much less (by 1E-3) than the real parts. Those parameter estimates are:

beta_est =

       3.1982e-009 +  2.7072e-012i
      19.2149e+000 +  3.0356e-015i
     162.7345e-003 -229.5976e-006i
     685.9447e+003                

absbeta_est =

       3.1982e-009
      19.2149e+000
     162.7347e-003
     685.9447e+003

I think it's time we contacted TMW Tech Support with this. If Origin can fit the function, MATLAB should be able to as well.

As for stopping the fit with ‘V’ > ‘Vg’, the only way I can think of is to add a statement at the beginning of ‘RcFun’:

if a(2) > Vg
    F = 0;
    return
end

That I also incorporated the latest run. The ‘return’ breaks out of the objective function without calling the Jacobian, and returns control to the calling program. It might throw an error if it wants the Jacobian at that step, but we can deal with that later if necessary.

Add this instead of your current ‘options’ statement before you call ‘nlinfit’:

options = statset('Display','final', 'TolFun',1E-100,'TolX',1E-100,'MaxIter',7500, 'MaxFunEvals', 5000, 'Jacobian', 'on', 'RobustWgtFun','bisquare');

It seems to improve the fit and parameter estimates significantly, approaching those Origin is producing. It might be interesting to experiment with various types of 'RobustWgtFun' options and 'Tune' tuning constants. I need to learn more about these.

The parameter estimates with these options are:

	Parameter Estimates:
		K =   4.013845660779272E-09    2.971810901348672E-11j
		V =   1.973624752002914E+01    1.566502377201617E-02j
		b =   8.384772342971474E-02   -2.018805502904597E-03j
		R =   3.163204974971465E+05    0.000000000000000E+00j

	Parameter Estimates:
		|K| =   4.013955674214642E-09
		|V| =   1.973625373683734E+01
		|b| =   8.387202334512081E-02
		|R| =   3.163204974971465E+05

The Jacobian is the multidimensional directional derivative that the gradient-descent algorithms use to determine the direction the algorithm proceeds to find the multivariate least-squares minimum. The default is to use a finite-difference numerical approximation of the Jacobian, but the analytic functions are always superior, and especially so when the response surface slope is shallow.

does 'RobustWgtFun' only work for nlinfit?

Also I get a warning (with nlinfit) that says the jacobian is near zero sometimes and that means the model is insensitive to some of the parameters which is clear as R never changes from the initial guess. Why is this? Do you get the same problem? Have you tried lsqnonlin?

If I use R = 100 i get

 K =

   4.8770e-09

 V =

   20.2103

 b =

    0.0164

 R =

   100

Or R = 2.32E6

K =

   8.7809e-10

 V =

   17.3575

 b =

    0.6139

 R =

     2320000

The final result is actually very sensitive to the initial guess

using x0 = rand(4,1).*[8E-10; 17; 0.6; 2.3E+6];

and doing the fit a few times actually can give very different answers, especially considering the fit does not allow R to change each time....

Perhaps if we can overcome this problem then we may hav cracked it?

Here is an example of the warnings

Warning: Rank deficient, rank = 3, tol =  7.443910e-06. 
> In nlinfit>LMfit at 378
In nlinfit>nlrobustfit at 532
In nlinfit at 215
In nlinfitRc at 26 
Iterations terminated: relative change in SSE less than OPTIONS.TolFun
Warning: Some columns of the Jacobian are effectively zero at the solution, indicating
that the model is insensitive to some of its parameters.  That may be because those
parameters are not present in the model, or otherwise do not affect the predicted
values.  It may also be due to numerical underflow in the model function, which can
sometimes be avoided by choosing better initial parameter values, or by rescaling or
recentering.  Parameter estimates may be unreliable. 
> In nlinfit at 253
In nlinfitRc at 26 
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =
4.718304e-17. 
> In nlinfit at 274
In nlinfitRc at 26 
Warning: Imaginary parts of complex X and/or Y arguments ignored 
> In nlinfitRc at 32 

Catching up ...

I believe 'RobustWgtFun' only works for ‘nlinfit’. I can't find any mention of it in the ‘optimset’ function that ‘lsqcurvefit’ uses. Note that ‘nlinfit’ is in the Statistics Toolbox and ‘lsqcurvefit’ in the Optimization Toolbox. The two Toolboxes and their functions don't seem to have much in common.

I got the same errors with ‘nlinfit’ as you did. I believe it's simply the nature of the function.

I haven't tried ‘lsqnonlin’ specifically. It uses the same algorithm and option structure as as ‘lsqcurvefit’, so I doubt there would be much difference in performance.

Trying to get some insight into the Jacobian, I created these plots:

Jcbn = @(a) ([1./((Vg - a(2)).^(a(3) + 1).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2), -(a(1).*(a(3) + 1))./((Vg - a(2)).^(a(3) + 2).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2), (a(1).*log(Vg - a(2)))./((Vg - a(2)).^(a(3) + 1).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2), -a(1).^2./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2]); 

Fit_J = Jcbn(beta);

Fit_J_a = abs(Fit_J);

figure(4099)
semilogy(Vg,Fit_J_a,'-')
% hold on
% hold off
legend('K','V','B','R', 4)
title('Jacobian Absolute Value')
xlabel('Vg')
ylabel('log_{e} |J|')
grid

figure(5000)
mesh(Fit_J_a)
set(gca, 'XTickLabel', 'K|V|B|R');
xlabel('Parameters')
grid on

I have the feeling that we may be approaching the point at which we've exhausted our options. The one idea I want to consider is the possibility of using a weighting vector to set the weights below 18V to zero and the rest to Vg or something else. That way it will ignore those with zero weight and weigh appropriately those we want it to. What parts of the function are most important for you to fit most closely?

We will solve this! I don't want to be responsible for your not being selected for FRS.

:)

No problem.

What I did remember during my interrupted sleep last night (thanks to thinking about this) was that one of the boundary conditions for the model is that Vd (which we have set as 1 here) must satisfy Vd<Vg-V so for us this means that Vg-V >= 1 and coincides with what you have said above.

Is there some way of including this condition in the regression with some kind of if loop or anything?

Wow

Well thats ruined mine and star striders fun.

Thanks

Couple of questions

Why log10 K and R? Why do the residuals have an oscillating pattern?

Might be able to move back on to my GUI now... Woo!

From what I understand, the log transformations normalize the parameters. I honestly have problems with that, although linear normalizing and scaling of the data — and not the parameters — are routine in most statistical applications.

I haven't run or plotted the posted code, but oscillating residuals — a problem I encountered earlier — indicates unmodeled behaviour in the fitted objective function. Any residuals with a non-random trend indicates a suboptimal fit.

I had problems with ‘lsqcurvefit’ miscalculating the Jacobian. I submitted a bug report to TMW Tech Support about it, since the same function and Jacobian worked perfectly in ‘nlinfit’. That won't have a significant effect on your problem since curve fitting doesn't seem to be the solution, but is an issue I want to resolve with TMW.

I have a couple errands to run, but will continue working on this with GlobalSearch later.

You could probably do without the logs, I'm sure it will still work just fine. But it just makes the numbers "nicer" and keeps all of the variables more or less of the same magnitude, which sometimes makes things easier numerically.

Instead of having to deal with very different numbers such as 1e-9 and 17 and 5000000, now the fitting routine can work with -9 and 17 and 6.7.

The residuals have an oscillating pattern because those oscillations were present in your original data, it's just hard to see them.

You can see this by plotting the gradient of your data, and comparing that with the gradient of the fitted function:

figure;
plot(diff(Id));
hold on;
plot(rcfun(BETA,Vg))

Ok, so the fit is probably good, or the model needs to include the oscilating trend in the data?

Thansk guys for all the help so far.

Would you both agree that the fit is good enough to use? Or should I keep working on optimising it?

Also would it be usefull to run RobustWgtFun or not for this instance? Should I keep using the Jacobian? and ow do I get the errors back to normal for the parameters that have been logged?

I could probably figure out some of these questions myself (;}) but i figured as you're here...

It seems to have trouble plotting without the logs?????

Right, there now seem to be a lot of options here to choose from.

Teja? I seem unable to change the code to allow me to get a good fit without using the rescalling, do you have the same problem

Star? Could you expain to me exactly what the jacobian function does in terms of making a better fit and what does this bit do

if nargout > 1
  J = Jfcn(a,Vg);
  if any(size(J) == 1)
      J = Jfcn(a,Vg');
  elseif ~any(size(J) == 1)
      return
  end
end

As there is not a unique solution, could this be a problem and also the solutions do not match that of Origin, but is this just because the fit is better as the residuals are no longer oscillating.

Does the code above still give imaginary numbers and is that a problem? Can we change the boundary condition to include the drain voltage?

Sorry for all the questions, i just want to get my head round this and use the best method to fit the data..

Thanks

Oh and I can't get your function work for me star, I'm not sure why...

The Jacobian is the directional derivative matrix that tells the gradient descent algorithm the direction in which to proceed to minimise the objective function with respect to each parameter. The gradient-descent algorithms all use the Jacobian for every problem, but usuallly use a finite-difference approximation to the Jacobian. The problem with this is that the approximation is only accurate to about 1E-10 or so. (I've read MATLAB documentation on this but I can't find it just now.) In situations such as yours, an analytic Jacobian is more accurate and more precise and as a bonus is also faster. If you want, you can experiment with that idea with this snippet with various values of ‘deltaw’:

        w1 = [0:0.1:2*pi]';
        deltaw = 1E-13;
        cosapp = (sin(w1+deltaw)-sin(w1))/deltaw;
        errs = [cosapp  cos(w1)  cos(w1)-cosapp]

It has nothing directly to do with your current problem but illustrates the difference between the accuracy and precision of finite-difference calculations and analytic calculations of derivatives.

The code bit you quoted has to do with changes I had to make to the objective function work with both the Statistics and Optimization Toolboxes.

First, I inserted the anonymous function ‘Jcbn’ because it's easier than duplicating the code for the Jacobian matrix in the ‘if’ block. I then call it with ‘J = Jcbn(a,Vg)’ or ‘J = Jcbn(a,Vg')’ later. The code is more efficient that way.

Second, the Optimization Toolbox does two things differently from the Statistics Toolbox. The Optimization Toolbox does not always call the Jacobian when it calls the objective function, so according to the documentation it's necessary to test to see if the calling routine wants the Jacobian or only wants the function. The Statistics Toolbox passes the data (in this case, Vg) to the objective function as a column, while the Optimization Toolbox passes it as a row, but both want a Jacobian matrix returned. When the Optimization Toolbox passes row-vector data to the Jacobian set up for the Statistics Toolbox, it returns a vector. In order for the objective function to return the Jacobian matrix to the Optimization Toolbox, it's necessary to transpose the data vector and calculate the Jacobian with it. This is what I spent a few hours yesterday figuring out.

The code I sent you was exactly the code that worked for me for both ‘nlinfit’ and ‘lsqcurvefit’. I copy-pasted it from the objective function I'm using. I have no idea why it didn't work for you, but I'll troubleshoot it with you. What about it doesn't work? Does it throw an error?

There is probably a ‘best’ solution, a problem I intend to use GlobalSearch to find, but the response surface is shallow enough that numerical methods can only approximate it. That's the reason I suggested you run the fit several times to get an idea of the range of acceptable parameter estimates. I then suggest doing some basic statistics on a set of 10-20 good fits to get an idea of how the parameters behave and how they are related.

The code and fits I quoted give real numbers for the parameter estimates and the function approximation, so no imaginary numbers appear. The residuals with the parameter estimates I quoted are very small and have no discernable oscillation.

I'm not sure what you mean by ‘change the boundary condition to include the drain voltage’. We're regressing drain current against gate voltage as I understand it. If we're going to add Vd as a parameter to be estimated, this would probably require rewriting the objective function and the Jacobian. Is there any way we could get around that by calculating it from ‘Id’ and ‘R’? It would be indirect rather than a specific parameter estimate, but considering the nature of the function we're fitting, might be preferable.

No worries with your questions. The ‘best’ method to fit the data may be GlobalSearch. That's the reason I'm going to give it a shot as soon as I teach myself to use it. In your particular situation, gradient-descent algorithms are likely not the best. By the same token, there may be no one ‘best’ way to approximate your or other problems with difficult response surfaces. We simply have to live with the uncertainties in the estimates. I note that the 95% confidence intervals calculated for some of the parameter estimate sets do not include parameters estimated in other runs (that appear to be as valid) so the parameter estimates may not be normally distributed. This is the reason I suggest you do several (perhaps 20 or so) runs, save the parameter estimates that are real and give good fits, and then explore their statistical distributions. (See ‘kstest’ for details on how to do this.) I may do that as well, but not until after I give GlobalSearch a go.

In that case would this be all you need to use for just nlinfit?

function  [F, J] = RcFun(a,Vg)

 if a(2) > Vg
    F = 0;
    return
end

 F = a(1)./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1));

 J = [1./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)) - (a(1).*a(4))./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2;  -(a(1).*(a(3) + 1))./((Vg - a(2)).^(a(3) + 2).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2);  (a(1).*log(Vg - a(2)))./((Vg - a(2)).^(a(3) + 1).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2);  -a(1).^2./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2]';
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The rest of my code then looks like

   x0 = rand(4,1).*[1E-13;  1;  1;  1E+7];  

 options = statset('Display','iter', 'TolFun',1E-30,'TolX',1E-30...
    ,'MaxIter',7500, 'MaxFunEvals', 5000, 'Jacobian', 'on'...
    ,'Robust','on','WgtFun','bisquare');

 [beta,resid,J,COVB,mse] = nlinfit(Vg,Id,@RcFun,x0,options);
 [ci se] = nlparci(real(beta),resid,'covar',COVB);

 Id_new = RcFun(real(beta),Vg);

But for some reason that gives me a really bad fit

Am I doing something fundamentally wrong?

If however I use

function F = rcfun(X,Vg)
F = (((10.^(X(1)).*(Vg-X(2)).^(X(3)+1)).^(-1))+10.^(X(4))).^(-1);
F = F.*(~imag(F));

I get the same fit every time, which looks ok and gives approximately the same estimates as matlab, but also gives the same oscillation in the residuals. Which I don't really like...

I believe we may be close to a saddle point or some other discontinuity in the response surface. I had to do the fit several times each to get the parameter estimates I posted. It won't always converge to the optimal parameter set, but the optimal sets give quite good fits. There aren't any discernable or repeatable oscillations in those residuals, and the residuals are very small in both absolute and relative magnitude.

If you're only going to use ‘nlinfit’ or other Statistics Toolbox functions, then use the function you quoted in your last comment. It's only necessary to use mine (with the ‘if’ blocks) if you want to use it with Optimization Toolbox functions.

This is a sort of heuristic effort. Just now, it took me six attempts to get one set of decent parameter estimates:

ci =

      73.5905e-012   385.7190e-012
      13.6982e+000    16.2602e+000
     819.6940e-003     1.2000e+000
       2.7308e+006     3.2005e+006

	Parameter Estimates:
		K =   2.296547485965766E-10    0.000000000000000E+00j
		V =   1.497918790161365E+01    0.000000000000000E+00j
		b =   1.009832375543948E+00    0.000000000000000E+00j
		R =   2.965617775329871E+06    0.000000000000000E+00j

	Parameter Estimates:
		|K| =   2.296547485965766E-10
		|V| =   1.497918790161365E+01
		|b| =   1.009832375543948E+00
		|R| =   2.965617775329871E+06

This is the rule more than the exception in nonlinear parameter estimation, especially for difficult problems such as yours. That's the reason I'm going to give GlobalSearch a go. I suspect (and I sincerely hope) that it will be able to go better than I can in finding the optimal set of parameter estimates.

I've just been trying to use lsqcurveift, but the problem seems to be that is doesn't like using the Jacobian.

I get errors of reference to non existing field JacobMult if I have 'Jacobian','on' and if I turn it 'off' i get a reference to non existing field JacobPattern?

With nlinfit I get a lot of warnings that say rank deficient, if I use my initial guess, if I use yours I get an error of

??? Error using ==> nlinfit>checkFunVals at 356
MODELFUN has returned Inf or NaN values.

 Error in ==> nlinfit>LMfit at 303
    if funValCheck && ~isfinite(sse), checkFunVals(r); end

 Error in ==> nlinfit at 171
    [beta_ls,J] = LMfit(X,y, model,beta,options,0,maxiter);

 Error in ==> nlinfitRc at 22
[beta,resid,J,COVB,mse] = nlinfit(Vg,Id,@RdFun,x0,options);

I might start having a look at this global search thing...

I have GlobalSearch running in the background. After running it a few times over a small region to be sure my code for it worked, I started it using:

x0 = [1E-15; 1; 1; 1E+9];

about 13 hours ago. It's still running now and hasn't found all the minima! On a smaller scale, it found 113 local minima, some with complex parameters, but none of those with complex parameters was a preferred fit.

The ‘lsqcurvefit’ function wants its Jacobian to be the transpose of the one ‘nlinfit’ wants. (That was the reason I had to insert the second ‘if’ block inside the first that ‘lsqcurvefit’ required.) I didn't get the error of ‘Inf or NaN values’, in my parameter estimation attempts, at least that I recall. The ‘lsqcurvefit’ function has the option of using a ‘trust-region-reflective’ algorithm that has advantages over the Levenberg-Marquardt algorithm, but it has some peculiarities and specific requriements I haven't studied. (I opted to learn about GlobalSearch instead. MATLAB is still running your and my GlobalSearch problem, so until it finishes, I can't do anything else with it.) You can specify the algorithm ‘lsqcurvefit’ uses.

While I'm thinking about it, the code I'm using for GlobalSearch is:

GSRcFun = @(a) a(1)./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1));
GlblSrchRcFun = @(a) sum([GSRcFun(a) - Id].^2);

problem = createOptimProblem('fmincon','objective',GlblSrchRcFun, 'x0',x1, 'lb',ones(4,1)*eps);
gs = GlobalSearch;
[xg fg flg og sl] = run(gs,problem)

I'm including this for your information so that if you find a computer with the Global Optimization Toolbox and want to run it yourself, you can use it rather than having to start from scratch.

I'll definitely let you know what GlobalSearch comes up with. I have no idea how long it will take for GlobalSearch to explore the response surface, discover all the local minima, and ideally identify the global minimum.

In that case, what version of Matlab are you using?

Maybe that's my problem?

Because otherwise I think I have exactly the same code as you for lsqcuvefit an nlinfit

I'm running 2012a with all the latest Toolbox versions (of the ones I subscribe to). Note that 2012b is still in beta test.

‘GlobalSearch’ is still running 13½ hours after it started ...

Quite funny,

At some point, i'm going to have hundreds of this kind of data to plough through. Not sure if global search is the way forward :)

I think i'll get on an download the latest version. Pretty sure we have a licence for it....

If you're going to have to analyse hundreds of such data, it definitely seems to be worthwhile to consider a more efficient strategy. (I thought this was a one off.) Is there any way to measure any of your parameters directly, or do a load line analysis or something similar?

Just an idea. I haven't need to do a load line in a very long time.

T+14 hours and counting ...

Pretty much I need to take this data of each device that I make. And any changes I make to the material, i'll need to make more measurements. I'm not sure about load lines but I don't think it will give me the info I need.

Anyway, I think there is one other way to get what I want, I'll be getting some data for the other method soon. But I would like to be able to use a few methods. Would you like some other data to try and fit?

I have managed to solve one of the problems

I think lsqcurvefit didn't like using statset for options, so writing the options like this solves the problem

optionsold = statset('Display','iter', 'TolFun',1E-30,'TolX',1E-30...
  ,'MaxIter',7500, 'MaxFunEvals', 5000, 'Jacobian', 'on'...
  ,'Robust','on','WgtFun','bisquare');

options = optimset(optionsold);

but lsqcurvefit will still not converge or change R.

With regards the other method for determining the parameters. What I have found out is that I can make some more measurements and calculate R separately, effectively taking it out of the equation and leaving us with just three parameters which for my issue would be brilliant as the only problem i have is due to trying to estimate R.

There's no statistical magic about my code or what I'm doing. I'm getting the same errors you are most of the time. It's simply that you have a difficult problem and a difficult data set, and the usual gradient-descent (and I suspect many other different) algorithms are going to have trouble with it. (I came late to statistics, by the way. My undergraduate concentration was Chemistry, and one professor remarked that ‘with Avogadro's number of molecules, you don't need statistics’! I learned differently later when I started physiological system modeling and parameter estimation and other fun activities.)

The ‘lsqcurvefit’ function wants ‘optimset’. I suggest you change the code you quoted above to something like:

options = optimset('Display','final', 'TolFun',1E-100,'TolX',1E-100,'MaxIter',10000, 'MaxFunEvals', 10000, 'Jacobian','on');

I suggest something like this for the call to ‘statset’ for ‘nlinfit’:

options = statset('Display','final', 'TolFun',1E-100,'TolX',1E-100,'MaxIter',7500, 'MaxFunEvals', 5000, 'Jacobian', 'on', 'RobustWgtFun','bisquare');

I'm not sure that the options structure for one would translate to the other.

I estimated R from the linear region of the fit by doing the linear regression:

B = polyfit(Vg, Id, 1);

and taking advantage of Ohm's law. In that region (that I would consider best to be bounded by 20V to 60V to avoid the discontinuity), we can take advantage of the relation that:

Id = Vg*B(1) + B(2),    where:  R = 1/B(1) and V = -B(2)/B(1)

Setting R ( i.e. a(4) ) to a constant using something like a linear regression to calculate it, and estimating only the other three parameters, would definitely improve the fit. It might also be worthwhile to consider restricting:

Vg > ceil(-B(2)/B(1))

(That's not valid MATLAB code. I simply formatted it that way.)

Let me know what you want to do with R and if and how you want to incorporate it as a constant. If it's going to be the same for all data sets, then it can be incorporated as a constant in the model, or it could simply be edited into ‘RcFun’ for each data set. Otherwise, there may be a way to pass it as an extra parameter. I'm considering these possibilities with respect to ‘nlinfit’. It would be easy with ‘lsqcurvefit’ since that would require simply setting a narrow bound for it.

I'm simply thinking out loud here as to the possibilities I'd consider. This is definitely your call.

Well. This is annoying.

Why does it work when I set the initial guess and model to include R as log10®?

Anyway, a little reluctantly I may have to go away and take some new measurements(could take a while) and then come back and try with a calculated R.

I just don't understand why Matlab can't converge R with the same algorithm as Origin? apparently nlinfit and origin both use Levenberg-Marquardt. In fact i did just try and change lsqcurvefit to the Levenberg-Marquardt algorithm but then I realised that that made it essentially exactly the same as nlinfit.

Do you think it could still be a problem with the complex numbers, because even with

 for a(2) > Vg
  F = 0

the parameters still come out complex and that then seems to effect nlparci as the jacobian comes out complex as well.

If you have any bright ideas let me know.

I noticed that you are on another thread realting to a fit with four parameters, so I will keep an eye on that, could come in handy...

I just though I should mention this as it may be important

But the original data went from 0 to -ve 60 (rather than +ve 60) and the equation includes a constant which I set to +ve 1 which was actually -ve 1. I did this because it seemed to be what they did in the paper and made the fit work in origin. But if you haven't already you could change the data and the equation to include these negatives. It may improve the fitting ability?

Another question.

How is it that you found the Jacobian and Taylor equations, did you do the differentiation yourself?

First, the other thread was an entirely different (and much easier) problem. I've done a fair amount of nonlinear parameter estimation and have made about every mistake possible by now (I hope), so I at least look at many of them and answer the ones I believe I can help with. Yours is the most difficult nonlinear parameter estimation problem I've ever encountered.

In some regions of your response surface, and for reasons I don't understand (I haven't analyzed them), some local solutions to your problem produce complex parameters and complex function estimates even when Vg (and Id) are restricted to the region where Vg > V. The GlobalSearch function found several local minima (I believe at one run it discovered nearly 200) and many with complex coefficients but rejected all of the complex coefficients as suboptimal fits.

Please send me the equation you referred to where you fitted the data from 0V to -60V and the constant that was actually -1. (I'm lost.)

I used the MATLAB Symbolic Math Toolbox to calculate the Jacobian and Taylor equations. It's easier than doing it by hand, and most importantly it avoids the inevitable algebra errors.

I'll try and explain. Hopefully not being too patronisingly. Basically we have a transistor where we apply a votage between a source and drain electrode (Vd). In the case of the paper I sent you they ad this set at -0.5 V. For my case it is -1V. We then sweep the gate voltage (Vg). In the paper the graph shows it going positive but i am pretty sure, due to the material they are using that it should be going negative. Which is what I did with mine, so I swepth the gate voltage from 20 to -60 v, i decided to just take 0 to -60 as it didn't make any difference. So the equation that describes the drain voltage is

 Id = Vd/(R+1/(CW/L)u(Vg-V)) i think?

Where C is capacitance, W is channel width and L is channel length and they are all known quantities. u is the mobility and we are saying it obeys a power law with the gate voltage u = k(Vg-V)^(b+1), combining the othe constants with k to give the k we use in the equation.

Anyway, i noticed that in the paper they seemed to of switched the polarity of theor gate votage, so I did it too, but that needed a switch to the drain votage as well, I think, so I did that. Seems to have been working. I think it also means that V is switched???

So really the equation should have a -1 at the beginning and the data should be reversed.

What did that

 F = F.*(~imag(F));

do? Could that be utalised to lose the imaginary parts of the function? Can it also be use on the jacobian?

I believe we did that before as well, without much difference. I may give it a go again to see what the result is. If I do, I'll likely recalculate the Jacobian to be sure it reflects the changes.

If we switched the polarities of ‘Id’ and ‘Vg’, negated the equation, and changed (Vg - V) to (V - Vg), I'm not certain we'd be any closer to a definitive solution unless the parameter values we're now estimating don't make sense ( i.e. agree with literature values except perhaps for sign ). I'm not sure that negating everything would get us anywhere.

I've designed biomedical instrumentation with transistors, but I've never done what you are doing with them. (Truth be told, I prefer FET op-amps.)

The statement (that I disagree with) that sets F = 0 for imaginary results:

F = F.*(~imag(F)); 

will introduce a significant discontinuity in the estimation procedure because the function isn't actually zero with imaginary parameters, and neither the analytic nor the finite-difference Jacobians would reflect the function arbitrarily being set to zero if it produced complex parameters or a complex fit to the real-valued data. (I don't know how the gradient-descent algorithm handles those problems.) However, restricting the fit region to Vg > 20 avoids the nonlinearity for V > Vg that the model fails to account for.

I may experiment with estimating R uniquely and fitting only the remaining three parameters. The easiest way to do that is likely to use ‘lsqcurvefit’ and restrict the upper and lower bounds for R to a narrow range around the estimated value. That would be easy to do in code with the current scripts.

Do you not find that R does not vary at all during the fit. In fact even when I estimate it to a high degree of accuracy (like 5 or 6 decimal places) It still stays the same.

Yes. It seems R varies least, and in the best-fitting parameter estimates I've recorded R varies from a bit more than 2E+6 to a bit less than 3E+6 Ohms. The latest (and best fit yet) parameter estimates are:

	Parameter Estimates:
		|K| =   1.102954516033334E-09
		|V| =   1.776473601749216E+01
		|b| =   5.426049408039866E-01
		|R| =   2.160189159613943E+06 

95%CI =

     843.3331e-012     1.3626e-009
      17.2964e+000    18.2331e+000
     470.0186e-003   615.1913e-003
       1.9953e+006     2.3251e+006

MSE = 4.873345E-019

Those with nlinfit. Other low-MSE nlinfit values for R I've recorded are: 2.594274479774597E+06, 2.589905421589671E+06, and 2.914483559512664E+06. For some reason, lsqcurvefit doesn't seem to do nearly as well as nlinfit with your problem.

Have you tried changing lsqcurvefit options to

 options = optimset('Algorithm','Levenberg-Marquardt'...

I assume it will come out with the same answers then

I just did a fit with nlinfit

 x0 = [2.6278E-10 15.215 0.97130 2.91E6];

beta
1.72920230731699e-10 - 1.65034719949234e-11i	
14.0401059416831 - 0.248767477152502i	
1.06796781086032 + 0.0225819426437233i	
2910000.00000000 + 0.00000000000000i

ci
1.51792677273391e-10	1.94047784190007e-10
13.6332651919646	14.4469466914015
1.02822032207520	1.10771529964544
2773641.38125087	3046358.61874914

R just doesn't budge at all

I haven't but I will. Good point. I was initially keeping lsqcurvefit with the trust-region-reflective option to compare it, but the Levenberg-Marquardt option now seems best.

I suggest bounding R in that situation to the range I quoted earlier: 2E+6 to 3E+6.

You can't bound R with Levenberg-Marquardt :(

After reading your last comment and re-reading the documentation for lsqcurvefit, I decided to experiment. First, I haven't been able to independently estimate R with regression or any other techniques. Second, the Taylor Series attempt to estimate the best initial parameter set failed.

Actually, nlinfit does reasonably well at estimating the parameters. It simply takes it a few attempts to get there.

If you're certain R isn't going to change in value between your organic FETs, it might simply be easier to set it as a constant and estimate the other three parameters.

Unfortunately, I'm pretty sure R will actually change with Vg, but we wont worry about that for the time being... Although it sounds to me like you're ready to call this a day???

What I cant get my head round is why nlinfit doesn't even bother changing R at all. Is this because it reaches a local minima before it has to? And if so could we set R = a(1) instead, and would that help?

Any chance you could explain how to use Symbolic Math Toolbox to caluculate the jacbobian?

This is what I have

F = k./(k.*R + 1./(x - V).^(b + 1));

jacobian(F,x)
Error: Expecting vector function in form of a list or vector.     [linalg::jacobian]

Also I have been looking into it and I think that scaling the parameters might be the only way to get the fitting to work neatly. Is this a bad Idea? and if not what is the best way to go about it?

Errands this morning. Late start here.

First, I believe I have a perfectly legitimate solution that will work for you now and in the future: linearly scale ‘Id’. I have the feeling that with ‘Id’ on the order of 1E-7, we're encountering Jacobian matrix condition and other numeric problems, as the errors have described. The disparate magnitudes of ‘Id’ and ‘Vg’ force ‘K’ and ‘R’ to the limits of numeric precision. Scaling ‘Id’ to be closer to ‘Vg’ avoids these extremes. Then to recover the actual magnitudes of ‘K’ and ‘R’, re-scale them by the same factor after the fit.

I experimented using ‘lsqcurvefit’ with this small change:

IdS = Id/median(Id);    % ‘IdS’ = ‘Id Scaled’ -- median(Id) = 1E-7

then with Lobnd and Upbnd = [ ]:

LQFopts = optimset('Display','final', 'TolFun',1E-12, 'TolX',1E-12, 'MaxIter',10000, 'MaxFunEvals', 10000, 'Jacobian','on');

[beta,resnrm,resid,xitflg,outpt,lmda,J] = lsqcurvefit(@RcFun,x0,Vg,IdS,Lobnd,Upbnd,LQFopts);

yields very quickly and with an excellent fit to the data:

ElpTimeSec = 1.0250e+000 

MSE = 34.9104e-006    % NOT RE-SCALED

	Parameter Estimates:
		K =   9.334400477356170E-03    
		V =   1.744670226837812E+01    
		b =   5.977271970075951E-01    
		R =   2.262890264338273E-01    

To re-scale ‘K’ and ‘R’:

K = K * median(Id); 

R = R / median(Id); 

that here are:

K = 9.2457143384208E-010

R = 2.2845962129575E+006

Parameters ‘V’ and ‘b’ don't have to be re-scaled, because the scaling of ‘Id’ doesn't affect them. Only scaling ‘Vg’ would, and there's no need to scale it here.

NOTE that this is not ‘scaling the parameters’. I'm scaling the data with a very basic linear transformation — this is done commonly — then rescaling the parameters affected by the scaling of ‘Id’, using the same scaling factor I used to scale ‘Id’ initially.

Second, with respect to calculating the Jacobian with the Symbolic Math Toolbox, here's what I did:

F = (((K.*(V-Vg).^(b+1)).^(-1))+R).^(-1)

J_MOSFET = jacobian(F, [K V b R])

J_MOSFET = simplify(collect(J_MOSFET))

J_MOSFET_a = subs(J_MOSFET, {K V b R}, {'a(1)' 'a(2)' 'a(3)' 'a(4)'})

vectorize(J_MOSFET_a)

that yields:

ans =

matrix([[1./((a(2) - Vg).^(a(3) + 1).*(a(1).*a(4) + 1./(a(2) - Vg).^(a(3) + 1)).^2), (a(1).*(a(3) + 1))./((a(2) - Vg).^(a(3) + 2).*(a(1).*a(4) + 1./(a(2) - Vg).^(a(3) + 1)).^2), (a(1).*log(a(2) - Vg))./((a(2) - Vg).^(a(3) + 1).*(a(1).*a(4) + 1./(a(2) - Vg).^(a(3) + 1)).^2), -a(1).^2./(a(1).*a(4) + 1./(a(2) - Vg).^(a(3) + 1)).^2]]) 

From here it's easy to convert it to an anonymous function such as the one I sent along earlier for it.

After more than a week of exploring various options and experimenting, this may be the solution you're looking for. Please let me know your thoughts.

Hmmm

It looked like it was going to work great, but I can't get it to do what yours did. I can only get it to work well if I turn 'jacobian','off' and I use F = F.*(~image(F)) (if I don't the fit is not as good)

With jacobian on i get the error

Error using eig
Input to EIG must not contain NaN or Inf.

are you using the same function?

When I use nlinfit it works with the jacobian on (but when I change the way the jacobian is written eg[ ; ; ;]' or [ , , ,] or [ ; ; ;] or [ , , ,]' it doesn't make any difference (i understand that some of those are the same)) With nlinfit i also need to use F = F.*(~image(F)) again to get a good fit...

You must be doing something else!/I must be doing something wrong!

Otherwise though I think this looks like it and I'd like to thank you very much for helping me out so much. I don't know if it was the most productive use of my time but I have enjoyed it and i'm much better with Matlab and am learning more about data analysis, which is good.

Maybe you could send me your m file?

I'm getting excellent fits, no condition problems with the Jacobian, and no imaginary numbers with the scaled Id. I sent you my calling routine lines for lsqcurvefit. (I'm sure nlinfit would work as well but I wanted to see what lsqcurvefit did considering that it was having more problems with the unscaled Id.) I consider the problem solved.

I described my new approach as fully as I could. The only thing I left out is that I am now using the generic:

x0 = rand(4,1); 

and am still getting excellent fits with scaled Id and no errors.

My full ‘RcFun.m’ m-file script is:

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

function [F,J] = RcFun(a,Vg)
% Adam Parry’s Organic FET Function
% DATE: 2012 07 23
%
%
%
% F = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1)
% 	      w.r.t:             K,     V,     b,     R
%         respectively       a(1),  a(2),  a(3),  a(4)
% x0 = [                 8.6E-10   ;1.7   ;0.8   ;5E6];

if a(2) > Vg
    F = 0;
    return
end

F = a(1)./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1));

Jfcn = @(a,Vg) [1./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)) - (a(1).*a(4))./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2;  -(a(1).*(a(3) + 1))./((Vg - a(2)).^(a(3) + 2).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2);  (a(1).*log(Vg - a(2)))./((Vg - a(2)).^(a(3) + 1).*(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2);  -a(1).^2./(a(1).*a(4) + 1./(Vg - a(2)).^(a(3) + 1)).^2]';

if nargout > 1
    J = Jfcn(a,Vg);
    if any(size(J) == 1)
        J = Jfcn(a,Vg');
    elseif ~any(size(J) == 1)
        return
    end
end

% =========================  END: RcFun.m  =========================

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

It's been interesting for me, too. I've learned a considerable amount about organic FETs, as well as data scaling as a solution to situations where the parameters assume extreme values. However I fail to see how it wasn't the most productive use of your time if you want to learn how to do parameter estimation in MATLAB and if this is an important part of your research.

If you're interested in learning more, I suggest an older but still very relevant text:

Beck, J.V.;Arnold, Kenneth J. Parameter Estimation in Engineering and Science (ISBN: 9780471061182).

It's pre-MATLAB but will tell you everything you need to know about nonlinear parameter estimation. There is also a lot on the Internet.

Thanks for all the help and I do think it has been worth while. Certainly quite fun. Just that I've been spending pretty much my entire time on Matlabb..

Just one thing. I have got it to work, results are exactly the same as yours, but only when I use data that starts at Vg>17. It seems that

if a(2) > Vg
  F = 0;
  return
end

doesn't stop my function giving complex numbers. Are you using data that goes from 0-60 in your fit?

It's my pleasure to help! I remember spending more time than I wanted to when I was first learning MATLAB (and being a bit frustrated in the process) because I also had other course work to do at the time. The benefit is that being proficient in MATLAB makes nearly everything else you do in engineering and maths significantly easier. It's definitely worth the effort.

I'm using:

Vg = Data(19:end,1);
Id = Data(19:end,2);

LinScl = median(Id);
IdS = Id/LinScl;

I'm avoiding the region of the discontinuity because the function doesn't account for behaviour in that region. (I get real parameter estimates for all values of Vg >= 18V but not for any less than that.) This is also entirely legitimate because there is no reason to believe the parameters are complex, and the function is not defined for (Vg < V) anyway. It is not legitimate to force complex parameters to be real simply because we want them to be real! Also, it is never appropriate to extrapolate outside the region of fit in any regression. So while my choice of data allows us to estimate ‘V’ [defined as a(4) = 17.29V in my latest attempt] as being outside the region of fit (that in my code starts with Vg = 18V), it does not allow us to estimate the function for Vg < 18V.