Non Linear Fit help?

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

Please post your code, objective function, and some data. (Adding it by ‘Edit’ to your original post would be best.)

What do you mean by the residuals being ‘a little funny’? What was the MSE? If Origin Pro gave you covariances or confidence limits on the parameters, please post them, too.

Where do the imaginary numbers appear — data, parameter estimates, or somewhere else?

I'm not sure i really know how to add code correctly. How would I add data?

On origin the MSE seems to be about x10^-19,

Reduced Chi squared 7.29777E-19

Adj, R-square 0.99985

The imaginary numbers come from the function as one of the parameters appears as a power which I think leads to roots of negative numbers. There may be a way to prevent this but I haven't looked into it too deeply yet.

You added your code correctly. To add some of your data (that seems to be an [n x 2] array), append it (by editing) to the code you already posted. I suggest for clarity something like:

DATA:

12.475  0.156
14.330  1.271

etc.

Please format you data as ‘code’ as I did here. That makes it easier to read.

Also, I see ‘myfun’ but I don't see ‘RcFun’. It would help if you posted it. Is there a particular reason you're taking the absolute value of ‘ydata’? Please post the original ‘xdata’ and ‘ydata’ rather than ‘Vg’ and ‘Id’.

I think in the data I gave above it dows not matter that I have taken the absolute, I don't think that it really matters. I changed myfun to Rcfun just copied the wrong file.

Just curious, but, what are the good values of K,V,b,R that you obtained through origin?

Ah. Sorry. You had it down below and I missed it. k = 8.6E-10 V = 17.3 b = 0.618 R = 2.3E6

Answer 1

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.

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 negated the function:

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

as well as negating ‘Vg’. (I copied it here so you can check it to be sure that's what you intended. The extra parentheses in it are likely not necessary, but I want to be sure it's doing what I want it to do. At this point, I'm taking no chances!)

The parameter estimates with those changes were:

beta_est =

     282.3452e-012 -335.6849e-012i
     -10.2542e+000 -  1.6440e+000i
     444.5718e-003 -125.7975e-003i
     845.8825e+003                

with ‘abs(beta_est)’ for information purposes only:

absbeta_est =

     438.6378e-012
      10.3851e+000
     462.0272e-003
     845.8825e+003

With a decent-looking fit of real(Id_new) as well as abs(Id_new).

I did a search on MOSFET models this afternoon (GMT-7) in order to understand a bit more precisely what you're doing. Unfortunately, I didn't learn much. I'd appreciate your providing me with an open-source PDF or website that can explain what you're doing if that's an option. I'm sensitive to your doing research and that you would likely not want to divulge too much to competing investigators.

I could send you a couple of PDF's over dropbox if you like? The fitting procedure and model is not confidential in anyway so that should be fine.

Could You post the whole of the code that you are using for me to have a look at?

Thanks

I would appreciate the PDFs. My circuit analysis and synthesis knowledge is relatively basic, but they could help me understand what you're doing.

This is the code snippet I've used:

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

        Vg = -Data(:,1);
        Id = Data(:,2);

        % x0 = [8.6E-10;1.7;0.8;5E6];
        x0 = rand(4,1).*[1E-10;  1;  1;  1E+6];
        % x0 = rand(4,1).*1E-1;
        % options = optimset('Display','iter', 'TolFun',1E-100,'TolX',1E-100,'MaxIter',5000, 'MaxFunEvals', 2500);    
        options = optimset('Display','final', 'TolFun',1E-100,'TolX',1E-100,'MaxIter',5000, 'MaxFunEvals', 2500);    
        [beta,resid,J,COVB,mse] = nlinfit(Vg,Id,RcFun,[x0],options);
        % Lobnd = [];
        % Upbnd = ones(4,1)*realmax;
        % Upbnd(2) = 0;
        % [beta,resnrm,resid,xitflg,outpt,lmda,J] = lsqcurvefit(RcFun,x0,Vg,Id,Lobnd,Upbnd,options);
        beta_est = beta 
        absbeta_est = abs(beta)
        se = sqrt(diag(COVB/size(Data,1)));
        if isreal(beta)
            ci = nlparci(beta,resid,'covar',COVB);
            % ci = nlparci(beta,resid,'jacobian',J);
        end

        Id_new = RcFun(beta,Vg);

        figure(4096)
        plot(Vg,real(Id_new),'r',Vg,Id,'o')
        hold on
        plot(Vg,imag(Id_new),'--r',Vg,Id,'o')
        plot(Vg,real(resid),'+')
        plot(Vg,imag(resid),'x')
        hold off
        grid

        figure(4097)
        plot(Vg,abs(Id_new),'r',Vg,Id,'o')
        hold on
        plot(Vg,abs(resid),'+')
        hold off
        grid

I didn't include the code and data I copied from your original post. I decided to keep in the lines I used for various diagnostic options, including my call to ‘lsqcurvefit’ and a variation on the objective function I tried.

Answer 2

https://dl.dropbox.com/u/11653970/Extracting%20Parameters%20from%20the%20Current-Voltage%20Characteristics%20of%20Organic%20Field-Effect%20Transistors%20-%202004%20-%20Horowitz%20et%20al.pdf

https://dl.dropbox.com/u/11653970/Temperature%20and%20gate%20voltage%20dependence%20of%20hole%20mobility%20in%20polycrystalline%20oligothiophene%20thin%20film%20transistors%20-%202000%20-%20Horowitz%2C%20Hajlaoui.pdf

Couple of links to pdfs

Also, the fits that I produce from nlinfit and lsqcurvefit still are no where near as good as on origin, do you get good fits, could you post a pic?

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?

I'm intrigued to know how origin decides to cut off the plot below Vg-V>Vd without it being specified. And manages it in only 18 or so iterations. In fact origin is able to get almost exactly the same parameters even setting the initial guesses quite far out....

so i tried a few things that maybe i shouldn't

I changed nlinfit.m so that the nlrobustfit function included an if else statement making it do a normal weighted robustfit if x-beta(2)>1 but would call a function in statrobustfit.m that set w = 0 when x-beta(2)<=1.

Where x is Vg and beta(2) is V

But that didn't really seem to do anything

Also I'm still getting the problem with beta(4) (which is R) not changing?

Do you think there is an easier way to implement the weighted fit?

I gave it some thought too overnight. (However I'm reluctant to change nlinfit or the others.) I thought about the weighting vector, and am considering:

WgtV = sqrt(Id);

thus forcing the routine to ignore the region where (Id = 0). I've done this successfully in the past with other functions and data, but in those situations using the usual inverse-variance weighting vectors, (clearly inappropriate here since there are no variances on the measurements).

The other thing I thought of is that since I have the Global Optimization Toolbox and a fair amount of experience with genetic algorithms (GA), I'm going to give that a go. We already have a good start on the fitness function (the objective function you derived), and setting up a weighted metric [to avoid the (Vg < V) and (Id = 0) problems]. An Euclidean distance measure is probably the easiest, although we can get creative and use something more robust if necessary. The problem is that GAs usually take a while to converge, but have the significant benefit that they are essentially immune to the response surface because they completely ignore it.

If we're not happy with the GA results, we'll give GlobalSearch a go, even though that may take longer than GA to converge. It is the most likely to produce the best solution. With luck, we'll have its results today or tomorrow.

These are problems I'll have to solve for my own applications in the not distant future, so I'm interested in learning about their solutions.

Answer 3

This problem is clearly difficult. It has multiple minima, and a strong dependence on initial conditions. Instead of NLINFIT, I tried two different approaches which both give me good fits for values of Vg > 18.

In both cases, since it seems this equation can't really model the Vg < 18 region, I weighted that region to be zero. Also, just to scale things better, I am working with the log10's of K and R.

Option 1: Pattern Search Algorithm from the Global Optimization Toolbox

This gives a solution that is almost identical to the one that you found from your other software.

G = @(X) (((10.^(X(1)).*(Vg-X(2)).^(X(3)+1)).^(-1))+10.^(X(4))).^(-1);
E = @(X) ( sum(((real(G(X)) - Id)).^2 .* (Vg >= 18)) );
opts = psoptimset;
opts.Display = 'iter';
opts.MaxIter = realmax;
opts.MaxFunEvals = realmax;
X0 = [log10(8.6E-10);1.7;0.8;log10(5E6)];
X = patternsearch(E,X0,[],[],[],[],[],[],[],opts);

% Your previous answer
K = 8.6E-10; V = 17.3; b = 0.618; R = 2.3E6;
F0 = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1);

% Comparing results: They are virtually identical
plot(Vg, Id);hold all; plot(Vg, G(X)); plot(Vg,F0);

Option 2: Simply using FMINSEARCH from various starting points

Again, identical to the previous results.

G = @(X) (((10.^(X(1)).*(Vg-X(2)).^(X(3)+1)).^(-1))+10.^(X(4))).^(-1);
E = @(X) ( sum(((real(G(X)) - Id)).^2 .* (Vg >= 18)) ); 
opts = optimset('fminsearch');
opts.MaxIter = realmax;
opts.MaxFunEvals = realmax;
options.TolFun = realmin;
options.TolX = realmin;
X0 = [log10(8.6E-10);1.7;0.8;log10(5E6)];

Emin = Inf;
for n = 1:100
    disp(n)
    [Xn, err(n)] = fminsearch(E,X0.*(1+0.01*randn(size(X0))),opts);
    if err(n) < Emin
        Emin = err(n);
        X = Xn;
    end
end
% Your previous answer
K = 8.6E-10; V = 17.3; b = 0.618; R = 2.3E6;
F0 = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1);

% Comparing results: They are virtually identical
plot(Vg, Id);hold all; plot(Vg, G(X)); plot(Vg,F0);

Although I don't have global optimisation toolbox and don't full understand whats going on, it looks like it could be on the right track.

What I would like to know is.

1. Is it possible to change the code so it can work for different data that will have a different V and so may work at Vg < 18. I.e can we ad a statement that confines Vg-V < Vd (Vd = 1 in the data I have provided but could be different for other data)

2. Being useless at statistics etc. How do I go about getting errors on the parameters and on the fit?

Kind Regards

Could you explain the

 Emin = Inf;
 for n = 1:100
    disp(n)
    [Xn, err(n)] = fminsearch(E,X0.*(1+0.01*randn(size(X0))),opts);
    if err(n) < Emin
        Emin = err(n);
        X = Xn;
    end
 end

Are you using MATLAB R2012a?

Answer 4

Ah! Sorry sorry sorry. There's actually a MUCH easier way to make all of this work. The important point is to set the model equal to zero whenever it goes imaginary. I bet your other software probably was doing this automatically; MATLAB does not. See how I modified rcfun below to include this line:

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

This sets F to be identically zero whenever it has an imaginary component.

Anyways, copy the entire function below into a file and run it from the command line:

>> [BETA,Rout,J,COVB,MSE] = domyfit

and it will work perfectly (and it should work for other data in general). The standard errors on your coefficients are given by diag(COVB).

function [BETA,Rout,J,COVB,MSE] = domyfit
A = struct;
A.data = [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];

xdata = A.data(:,1);
ydata = A.data(:,2);
Vg = xdata;
Id = abs(ydata);

X0 = [log10(8.6E-10);1.7;0.8;log10(5E6)];
[BETA,Rout,J,COVB,MSE] = nlinfit(Vg,Id,@rcfun,X0);
K = 8.6E-10; V = 17.3; b = 0.618; R = 2.3E6;
F0 = (((K.*(Vg-V).^(b+1)).^(-1))+R).^(-1);
plot(Vg, Id);hold all; plot(Vg,rcfun(BETA,Vg)); plot(Vg,F0);

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

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?????

I'm going to continue experimenting with GlobalSearch, since your problem is likely the best test for it that I've encountered. The Origin fit didn't give much oscillation in the residuals if I remember correctly.

As for the fit being good enough to use — that's definitely your call! I'm going to see what GlobalSearch comes up with. That's likely as close to a definitive set of parameter estimates as is possible. (I'm not happy with the fits the curve-fitting routines have come up with so far.) They may be the set that Origin estimates, but at least we'll know.

I suggest continuing to use RobustWgtFun and the Jacobian. If you use the Jacobian with ‘lsqcurvefit’, insert this line before you calculate the Jacobian in RcFun:

Vg = Vg';

The MATLAB convention is column-major order, and nlinfit passes column vectors (as I would expect) to its objective functions, so I have no idea why lsqcurvefit passes row vectors. (I figured this out myself by experimenting with it.)

I've successfully estimated equivalent circuit models (to physiological problems) with pF capacitances, mH inductances, and MOhm resistances with MATLAB curve-fitting routines without having problems (although it usually took several guesses of initial parameters to get them to fit), so I doubt parameter magnitude ranges are the problem. My data oscillated so were not as challenging to fit as yours are.

I'm quite enjoying this.

What is it you mean by experimeting with GlobalSearch, i assume (by looking into it abit) that it finds global minima in a problem?

Could you use MultiStart?

Also (to Teja in particular) how do you go from log10 parameters to estimating the error in the parameter itself using COVB?? I just cannot work it out

SUCCESS!

After all that, it came down to the initial conditions, although I also made some changes to ‘RcFun’ to make it work seamlessly with both ‘nlinfit’ and ‘lsqcurvefit’.

The new initial conditions:

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

and the new objective function:

function [F,J] = RdFun(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  =========================

produces the parameters:

	Parameter Estimates:
		K =   3.254448861356919E-10    0.000000000000000E+00j
		V =   1.559046341560319E+01    0.000000000000000E+00j
		b =   9.096959820500091E-01    0.000000000000000E+00j
		R =   2.828879053496115E+06    0.000000000000000E+00j

	Parameter Estimates:
		|K| =   3.254448861356919E-10
		|V| =   1.559046341560319E+01
		|b| =   9.096959820500091E-01
		|R| =   2.828879053496115E+06

and it only took ‘nlinfit’ 2.075 seconds to converge.

However they do not appear to be unique because these produce and equivalently close fit:

	Parameter Estimates:
		K =   6.692378493392830E-10    0.000000000000000E+00j
		V =   1.685567986182001E+01    0.000000000000000E+00j
		b =   6.960550385747625E-01    0.000000000000000E+00j
		R =   2.476857099346900E+06    0.000000000000000E+00j

	Parameter Estimates:
		|K| =   6.692378493392830E-10
		|V| =   1.685567986182001E+01
		|b| =   6.960550385747625E-01
		|R| =   2.476857099346900E+06

as do these:

	Parameter Estimates:
		K =   3.054487171664253E-10    0.000000000000000E+00j
		V =   1.547909799526304E+01    0.000000000000000E+00j
		b =   9.280342803415647E-01    0.000000000000000E+00j
		R =   2.854955008915018E+06    0.000000000000000E+00j

	Parameter Estimates:
		|K| =   3.054487171664253E-10
		|V| =   1.547909799526304E+01
		|b| =   9.280342803415647E-01
		|R| =   2.854955008915018E+06

that produced these 95% confidence intervals:

     137.6613e-012   473.2361e-012
      14.4356e+000    16.5226e+000
     772.1086e-003     1.0840e+000
       2.6426e+006     3.0673e+006

and this covariance matrix:

       8.2089e-021    49.8600e-012    -7.6134e-012    -9.6597e-006
      49.8600e-012   317.5260e-003   -45.6541e-003   -54.6273e+003
      -7.6134e-012   -45.6541e-003     7.0893e-003     9.1732e+003
      -9.6597e-006   -54.6273e+003     9.1732e+003    13.1525e+009

although sometimes it still had a bit of a problem converging. The residuals were reasonably flat with no significant trend, at least in my opinion. I guess that's the result of the ‘rand’ component of the initial parameter guess vector. I'm still going to work on the GlobalSearch, but I feel reasonably confident about these, considering the fit and the residuals. Unfortunately, ‘lsqcurvefit’ consistently failed to find a fit with the same initial parameter guesses. I'll let you know what GlobalSearch comes up with. It may take a few days for me to learn about and successfully run it.

You might want to experiment a bit with it yourself to see what the various parameter estimates look like. That will give you a better idea than even the confidence intervals will about the region of fit.

Answer 5

"How do you go from log10 parameters to estimating the error in the parameter itself using COVB??"

It's very easy, you just use the chain rule.

My change of coordinates was y = 10^x

So then dy/dx = 10^x * log(10)

Let A = Covariance in regular coordinates (this is what you want) and B = Covariance in log coordinates (this is what you have)

Then you can find A as

A = 10.^(B) * log(10)

In hindsight, I think this was all a bad idea, and you should just do it without any rescaling, so you don't have to worry about these kinds of things.

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?

Do you mean ‘other method’ of device fabrication or parameter estimation?

I would be interested in other data to fit, yes.

I finally gave up on GlobalSearch, at least for now. I've since set a two-hour time limit on it. I would have done that earlier but I had no idea it would still be running after 18 hours!

Out of sheer desperation, I asked the Symbolic Math Toolbox to calculate a second-order Taylor series approximation of your function to see if that would make it easier to estimate the parameters, at least initially. The problem is that the Taylor series wants to calculate (-V).^b, guaranteeing a complex result. Here's the Taylor series, formatted by the Symbolic Math Toolbox as an anonymous function:

Taylr_Id = @(K,R,V,Vg,b)K.*(-V).^b.*1.0./(K.*R.*V.*(-V).^b-1.0).^2.*(-V+Vg+Vg.*b+K.*R.*V.^2.*(-V).^b)

that after this substitution: {K V b R}, {'a(1)' 'a(2)' 'a(3)' 'a(4)'} (and apparently some rearranging) becomes:

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

I doubt if it contributes anything of significance to the effort, but so long as I have it I'll send it along for you to experiment with.

Yeah there is another method for parameter estimatio, but it involves taking some other measurements, but it may only require fitting straight lines, which would be nice. I'll have to look into it.

In other news I still can't get rid of this error without using log10(parameters). aAnd I still can't figure out how to get the correct errors from them.

Warning: Rank deficient, rank = 2, tol = 2.661976e-08. > In nlinfit>LMfit at 351 In nlinfit>nlrobustfit at 532 In nlinfit at 215 In nlinfitRc at 22

I'm managing to get fits with your code but they never see to give me residuals that are random, and the fits are still a little off in comparison to origin. Am I doing this right?

The nice thing is that fitting straight lines only require two parameters. However you have four in your model, so I don't see what the advantage is in that unless you can ignore two of them. When I did this just now, I got estimates for R of 195.8839E+006 and for V of 19.9989E+000.

I believe the nature of your data and the nature of your model make fitting your data difficult. The model doesn't really account for the sharp discontinuity where Vg=V, or for the very small oscillations in the Id data, that are likely real although quantitatively unimportant. (When I looked at them on the model I estimated from the parameter sets I posted, the magnitude of the oscillations with the best parameter estimates was about 1% of the value of the function at that point.)

The problem is that taking the logs of the parameters and then reconsituting them in the model changes the model. It's easier to see if instead of using 10.^a(1) and 10.^a(4) you substitute exp(a(1)) and exp(a(4)). That's not the model you specified and it doesn't make any sense in the context of the papers you referred me to, so I am not using the parameter transformation.

A Taylor series approximation might still be the way to go initially to get a good start on the parameter estimates, but you will have to change your data and model to avoid the possiblity of ‘(-V).^p’. Changing ‘(Vg-V)’ to ‘(V-Vg)’ does this, and likely only requires changing the Vg data to -Vg. In the Taylor series I sent along, to make that transformation simply change the signs of ‘V’ (a(2)) and ‘Vg’. I suggest that only for the Taylor approximation in order to get a good initial set of stable parameter estimates.

The parameter estimates I posted earlier are among the best I got. You can plug them in directly and calculate the relevant statistics. There is a small oscillation in the residuals, but it is barely discernable. This set took less than three seconds to converge and produced an MSE = 1.3182E-018:

	Parameter Estimates:
		|K| =   2.627883069159599E-10
		|V| =   1.521501224663270E+01
		|b| =   9.713005068703799E-01
		|R| =   2.914483559512664E+06

I gave up on GlobalSearch because GlobalSearch either didn't find a global minimum or got caught in a loop that required me to either stop or completely exit MATLAB, and so about which it gave me no informaton. It seemed to have the most trouble fitting ‘R’ (it seemed to persisstently underestimate it), but didn't seem to have trouble with the other parameters.

I suggest you go with the parameter estimates that make the most sense and that you are the most comfortable with. I know nothing about Origin, how it decides the region to fit, or how it calculates its parameters, so I can't advise you on what parameter set is best.

I'll be glad to continue to provide what help I can.

Answer 6

Hi

Long time

To be honest. I think that I am just about understanding the statistical things that you are doing in order to get a good fit. But my main problem is actually writing the code that does the same job as yours is doing. I am still unable to get lsqcurvefit to work with Jacobian on and nlin fit still gives me Rank deficient warnings. The main problem that I have though is that no matter waht initial parameters I start with R does not change. You mentioned a similar thing you noticed with R above and I just wonder how I can overcome it?

I know this is not exactly what you want to do, but could you walk me through the code. I can show you what I've got so far if that helps?

Thank you

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.

Ok, thats fine, do you know anyway that Matlab could work out for itself when not to take into account data?

Like if I test a device that has a threshold voltage (V) around 0.5, i'm going to want most of the data up until 0 aren't I?

This is going to be a bit cheeky, but could you tell me if this jacobian is right, I'm getting an error telling me it is the wrong size

F = (a(1).*((VgX-a(2)).^(a(3)+1)));

J = [(VgX-a(2)).^(a(3)+1),...
  -(a(1).*((VgX-a(2)).^a(3)).*(a(3)+1)),...
  a(1).*log(VgX-a(2)).*((VgX-a(2)).^(a(3)+1))];

I don't know how MATLAB routines could decide for themselves what part of the data to use. I have no idea how Origin managed to fit your function as well as it did and only over the linear region of fit that produced real parameters, assuming you gave it and MATLAB the same function and data.

I don't know the data you're commenting on. Since I know only the data from the present problem, I suggest you start the fit with the greatest ‘Id’ data, take a range of those data, do the fit, and estimate ‘V’. Then take the data greater than the estimated ‘V’ and fit all of it. If your parameters are linear functions of your data (as ‘K’ and ‘R’ were here), and if your ‘Id’ data are again on the order of 1E-7, you can scale the ‘Id’ data and then ‘un-scale’ ‘K’ and ‘R’ by the same factor. Since ‘V’ and ‘b’ are not linear functions of either variable, you have to fit them as they are. So you can't scale ‘Vg’.

Your Jacobian looks fine to me. The Optimization and Statistics Toolboxes use different conventions — Statistics uses column-major, Optimization row-major — so to get your Jacobian to work with both, you have to transpose your ‘VgX’ vector to use the Optimization functions. If your data are column vectors, the Jacobian as listed should work with Statistics Toolbox functions, but you'll have to transpose the ‘VgX’ vector to make it work with Optimization Toolbox functions. Take the ‘RcFun.m’ I sent you, rename it to your present problem, then paste your ‘F’ and ‘J’ functions in it in place of the ones I provided. It should work fine. Note that my variable ‘Jcbn’ is an anonymous function, so be sure to paste your ‘J’ expression there and change the argument from ‘Vg’ to ‘VgX’.

That should work.