Error values in motofit

Share icon

I have been using motofit to do curve-fitting for a while now. Some functions will yield +/- values, for instance:

Genetic Optimisation Successful
  Fitting: YCumData to MOTO_NewGlblFitFunc
  V_fitIters = 420; V_Chisq = 96.7877; V_npnts= 162; V_nterms= 4; V_nheld= 0; V_logBayes = 65.7031
    w[0]    =   3001.7   +/-   0.249107
    w[1]    =   9.91829e-07   +/-   8.10371e-10
    w[2]    =   0.9985   +/-   0.000139441
    w[3]    =   199.641   +/-   0.164732
Genetic Optimisation Successful
  Fitting: YCumData to MOTO_NewGlblFitFunc
  V_fitIters = 420; V_Chisq = 96.7877; V_npnts= 162; V_nterms= 4; V_nheld= 0; V_logBayes = 65.7031
    w[0]    =   3001.7   +/-   0.249107
    w[1]    =   9.91829e-07   +/-   8.10371e-10
    w[2]    =   0.9985   +/-   0.000139441
    w[3]    =   199.641   +/-   0.164732
  _______________________________
  Global Fit converged normally.

while other functions do not, for instance:

Genetic Optimisation Successful
  Fitting: YCumData to MOTO_NewGlblFitFunc
  V_fitIters = 287; V_Chisq = 212.499; V_npnts= 162; V_nterms= 10; V_nheld= 0; V_logBayes = 0
    w[0]    =   28.0851   +/-   0
    w[1]    =   5.48056   +/-   0
    w[2]    =   3593.68   +/-   0
    w[3]    =   18178.4   +/-   0
    w[4]    =   3.41564   +/-   0
    w[5]    =   9.834e-07   +/-   0
    w[6]    =   1.65546e-10   +/-   0
    w[7]    =   1.35591e-09   +/-   0
    w[8]    =   9.91592e-12   +/-   0
    w[9]    =   200.196   +/-   0
  _______________________________
  Global Fit converged normally.


Here is the function for the function giving +/- values
Function ZrealQandlesU(w,frequency) : FitFunc
    Wave w
    Variable frequency

    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Equation:
    //CurveFitDialog/ f(frequency) = -1*(2*Cdl*frequency*3.14159)/(4*Cdl^2*frequency^2*3.14159^2+1/(RP^2))
    //CurveFitDialog/ End of Equation
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ frequency
    //CurveFitDialog/ Coefficients 4
    //CurveFitDialog/ w[0] = RP
    //CurveFitDialog/ w[1] = Qdl
    //CurveFitDialog/ w[2] = ndl
    //CurveFitDialog/ w[3] = RS

    Variable RP,Qdl,RS,ndl

    RP=w[0]
    Qdl=w[1]
    ndl=w[2]
    RS=w[3]

    return RS+ real(1/(1/RP +Qdl*(2*Pi*frequency*sqrt(-1))^ndl))
End


And the function that does not:
Function ZrealTlporosity(w,frequency) : FitFunc
    Wave w
    Variable frequency

    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Equation:
    //CurveFitDialog/ f(frequency) = transmission line
    //CurveFitDialog/ End of Equation
    //CurveFitDialog/ Independent Variables 1
    //CurveFitDialog/ frequency
    //CurveFitDialog/ Coefficients 10
    //CurveFitDialog/ w[0] = R1
    //CurveFitDialog/ w[1] = R2
    //CurveFitDialog/ w[2] = Ra
    //CurveFitDialog/ w[3] = Rb
    //CurveFitDialog/ w[4] = R0
    //CurveFitDialog/ w[5] = Ca
    //CurveFitDialog/ w[6] = Cb
    //CurveFitDialog/ w[7] = C0
    //CurveFitDialog/ w[8] = L1
    //CurveFitDialog/ w[9] = Rsb
   
    Variable R1,R2,Ra,Rb,R0,Ca,Cb,C0,Rsb,x1,x2,L1
    Variable/C pr,Za,Zb,lambda,Cl,Sl,Ztline
   
    R1=w[0]
    R2=w[1]
    Ra=w[2]
    Rb=w[3]
    R0=w[4]
    Ca=w[5]
    Cb=w[6]
    C0=w[7]
    L1=w[8]
    Rsb=w[9]
   
    pr = 1/(1/R0 - (sqrt(-1)*2*pi*frequency*C0))
    x1 = R1
    x2 = R2
    Za = 1/(1/Ra - sqrt(-1)*2*pi*frequency*Ca)
    Zb = 1/(1/Rb - sqrt(-1)*2*pi*frequency*Cb)

    lambda = (pr/(x1 + x2))^(1/2)
    Cl = Cosh(L1/lambda)
    Sl = Sinh(L1/lambda)
       
       return Rsb+real(1/(x1 + x2)*(lambda*(x1 + x2)*Sl + (Za + Zb)*Cl + (Za*Zb*Sl)/(lambda*(x1 + x2)))^(-1)*(L1*x1*x2*(x1 + x2)*Sl + x1*(lambda*x1*Sl + L1*x2*Cl)*Za + x2*(lambda*x2*Sl + L1*x1*Cl)*Zb + (2*X1*x2 + (x1^2 + x2^2)*Cl + (L1/lambda)*x1*x2*Sl)*(Za*Zb)/(x1 + x2)))
End


Any idea as to why this is?
Dear mtaylor,
I've looked into this and I've got some good and bad news.
THe bad news first. The code is determining the fit uncertainty by estimating the covariance matrix using a gradient technique (similar to the Levenberg Marquardt). This involves determining the Hessian (partial derivative of Chi2 w.r.t each of the parameters) followed by matrix inversion.
The code is having a problem inverting the matrix, probably due to numeric instability caused by the presence of some very small parameters (1e-12). I could probably condition the matrix better by scaling those values up, but it's not a high priority fix for me as you still get good fit parameters and it's not the prime output of the genetic optimisation. Following this analysis with a LM analysis may give you the fit uncertainties you desire, in this case.

THe good news is that you can speed your fits up a lot by choosing a population size multiplier of 10, and a tolerance of 0.01. The values you had in there were much too small and were making the fit take way too long.

Please note that the genetic optimisation itself does not use gradients.

December 16, 2011 at 12:48 am - Permalink

Andy,

Curiously, running the levenberg method yielded this result:

 Coefficient values ± one standard deviation
    K0  =28.085 ± 7.13e+05
    K1  =5.4806 ± 5.89e+05
    K2  =3593.7 ± 1.77e+04
    K3  =18178 ± 4.25e+05
    K4  =3.4156 ± 4.65e+05
    K5  =9.834e-07 ± nan
    K6  =1.6555e-10 ± nan
    K7  =1.3559e-09 ± 171
    K8  =9.9159e-12 ± 0.00232
    K9  =200.2 ± 0.194
  Hmm... Global Fit stopped for an unknown reason.


Which doesn't look so useful ... I think that despite reporting values, that the inversion is again failing, probably due to the nature of my function.

I appreciated your advice about the fit variables. I did not realize that the population size was a multiplier! I have been adjusting popsize based on the number of coefficients that I was using. For some functions, this meant using popsize=400. Clearly, for 20 coefficients, that would actually generate a population of 8000, which slows things down a bit. I am going to read your documentation before asking any further questions in order to respect your time. :)

Thanks again

Matt

December 16, 2011 at 01:47 pm - Permalink

Matt-

The fit result you posted looks like you have "Identifiability" problems. That's regression buzz talk for a situation where two or more of your fit coefficients are nearly (or actually) linearly dependent. That means that if you were to perturb one of the affected coefficients and hold it, then do a fit you would get a fit that is nearly as good as the original solution. In that case, the fitting algorithm can't distinguish the different solutions and the chi-square surface around the solution is very flat. That leads to those very large error bars.

Unfortunately, it's not something that is easily fixed. It's a mathematical problem, not a software or algorithm problem. Sometimes you can relieve the situation by adding more data points to cover some area in the coefficient space where the coefficients have different effects. Sometimes you can reformulate the model to reduce the number of coefficients.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

December 20, 2011 at 10:20 am - Permalink

Checkout the
Moto_montecarlo(fn, w, yy, xx, ee, holdstring, Iters,[cursA, cursB, outf]) function in the geneticoptimisation.ipf file. This may be able to assist in evaluating parameter uncertainties (but not in a global fit).

December 20, 2011 at 02:33 pm - Permalink

Thank you for the information and advice. The function in this case is being fit to data that is for a different system; I am using this data while I am troubleshooting my functions. It is possible that a better situation will emerge when I use appropriate data, however, I am surprised with this result as certain values should go to zero since it is a linear combination of the function that created the data (Randles Cell) and another (transmission line).

December 21, 2011 at 11:41 am - Permalink

When doing a global fit, is there any way to add switches (like /DUMP, or /POL)? At this point, I am editing the procedure file to do this...

Matt

December 21, 2011 at 12:10 pm - Permalink

Another question about global fit when using motofot: Does the epsilon value get taken into account? I know it can be useful to tune granularity when using a gradient method.

December 21, 2011 at 12:12 pm - Permalink

With Gencurvefit there are no such thing as epsilon values, differential evolution is a stochastic process not gradient based. You can't add in flags apart from editing the procedure file.

December 21, 2011 at 01:04 pm - Permalink