Chi-squared minimization for a model with 4 parameters

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
Report Content as Inappropriate

Chi-squared minimization for a model with 4 parameters

This post has NOT been accepted by the mailing list yet.
I've been trying to fit a model with 4 parameters to my data for a while now, by chi-squared minimization. I've tried curve_fit, minimize, least_squares, I've even tried lmfit without success.  I've tried different minimizers and methods, bounded and unbounded ones, as well as rescalling my data or modifying the range of the fit. I still can't get a decent fit (as easily seen on plotted charts).
I am sure that my model is correct, as well as my constants and data. I've been shown that it is possible to solve this problem using Excel's Solver (on the exact same data). It's hard for me to imagine that it's better that scipy.optimize. Do I have any errors in my code? Do you have any suggestions? I can provide more data if neccesary, or look for an open access article with the model. I'm providing test data and the chart with results as attachments.
Here's my code:
import matplotlib.pyplot as pyplot
import numpy
import pandas
from scipy.optimize import minimize

# Constants

# sample length/width, unitless
sampleDimension = 6
# electron charge in [C]
echarge = 1.6021766208 * 1e-19
# vacuum permittivity in [F/m]
epsilonZero = 8.854187817 * 1e-12 * 1e-2  # *(1e-2)  to make F/cm from F/m
# relative electric permittivity of the substrate: 3.9 SIO2
epsilon = 3.9
# oxidant thickness in [m]
tox = 285 * 1e-9 * 1e2  # *1e2 to make [cm]
#tox = 285
# capacitance of the oxidant on the gate, per unit of surface area
cox = epsilon * epsilonZero / tox


# Reads data from a file and cuts it off at a specific place
resistance = numpy.array([])
voltages = numpy.array([])
resistance_error = numpy.array([])
currents = numpy.array([])
rescale = 1

with open("testData.txt") as dataFile:
    data = pandas.read_csv(dataFile, sep='\t', decimal=',').values
    voltages = numpy.array(data[:, 0])
    currents = numpy.array(data[:, 1])
    currentsErr = numpy.array(data[:, 2])
    resistance1 = numpy.array([a / b for a, b in zip(abs(voltages), currents)])
    resistance2 = numpy.array([a / b for a, b in zip(resistance1, currents)])
    resistance_error = numpy.array([a * b for a, b in zip(resistance2, currentsErr)])

    resistance = numpy.array([a / b for a, b in zip(abs(voltages), currents)])
    #rescale = resistance.max()

    resistance = resistance / rescale
    resistance_error = resistance_error / rescale

def modelFunction(parameters):
    mobility, rcontact, n0, vdirac = parameters
    model = 2 * rcontact + (sampleDimension / (numpy.sqrt(n0 ** 2 + (cox * (voltages - vdirac)) ** 2) * mobility))
    return model/rescale

def chisqFuction(initial_parameters):
    theory = modelFunction(initial_parameters)
    chisq = numpy.sum(((resistance - theory) / resistance_error) ** 2)
    return chisq

# set initial parameters - mobility, rcontact, n0 * echarge, vdirac
initial_parameters = numpy.array([3e3, 1e5, 1e12*echarge, 65])

# set bounds for fitting parameters
bnds = ((1e2, 2e5), (1e1, 1e8), (1e8*echarge, 1e14*echarge), (55, 75))

# fitting to data using leastsq method, using minimum tolerance and displaying
result = minimize(chisqFuction, initial_parameters, bounds=bnds, method='L-BFGS-B', options={'disp': True})
best_parameters = result.x
print("initial parameters:", initial_parameters)
print("parameters:", best_parameters)
scatter = pyplot.scatter(voltages, resistance*rescale)
initFitLine, = pyplot.plot(voltages, modelFunction(initial_parameters), 'k--')
bestFitLine, = pyplot.plot(voltages, modelFunction(best_parameters), 'r-')
pyplot.xlabel('Gate voltage [ V ]')
pyplot.ylabel('Resistance [ M\u2126 ]')
pyplot.legend([scatter, bestFitLine, initFitLine], ['Data', 'Best Fit', 'Init fit'], loc='upper left')

And here's output from minimize:
initial parameters: [  3.00000000e+03   1.00000000e+05   1.60217662e-07   6.50000000e+01]
parameters: [  1.21380486e+02   1.98859683e+05   1.60217662e-11   5.50000000e+01]

           * * *

Machine precision = 2.220D-16
N =            4     M =           10

At X0         0 variables are exactly at the bounds

At iterate    0    f=  4.18981D+08    |proj g|=  9.29832D+02

At iterate    1    f=  4.18397D+08    |proj g|=  9.23872D+02
  ys=-4.324E+05  -gs= 4.574E+05 BFGS update SKIPPED

At iterate    2    f=  3.82981D+08    |proj g|=  6.02007D+02
  ys=-6.973E+08  -gs= 6.416E+06 BFGS update SKIPPED

At iterate    3    f=  3.75899D+08    |proj g|=  5.18560D+02

At iterate    4    f=  3.55938D+08    |proj g|=  2.40624D+04

At iterate    5    f=  3.55843D+08    |proj g|=  7.15256D+01

At iterate    6    f=  3.55727D+08    |proj g|=  3.57628D+01

At iterate    7    f=  3.55726D+08    |proj g|=  1.20997D+03

At iterate    8    f=  3.55726D+08    |proj g|=  3.57628D+01

At iterate    9    f=  3.55725D+08    |proj g|=  2.38419D+01

At iterate   10    f=  3.55725D+08    |proj g|=  2.38419D+01

At iterate   11    f=  3.55725D+08    |proj g|=  4.17233D+01

           * * *

Tit   = total number of iterations
Tnf   = total number of function evaluations
Tnint = total number of segments explored during Cauchy searches
Skip  = number of BFGS updates skipped
Nact  = number of active bounds at final generalized Cauchy point
Projg = norm of the final projected gradient
F     = final function value

           * * *

   N    Tit     Tnf  Tnint  Skip  Nact     Projg        F
    4     11     28     15     2     2   4.172D+01   3.557D+08
  F =   355725362.52186227    


Cauchy                time 0.000E+00 seconds.
Subspace minimization time 0.000E+00 seconds.
Line search           time 0.000E+00 seconds.

Total User time 0.000E+00 seconds.

Nonpositive definiteness in Cholesky factorization in formt;
   refresh the lbfgs memory and restart the iteration.

Nonpositive definiteness in Cholesky factorization in formt;
   refresh the lbfgs memory and restart the iteration.
  ascent direction in projection gd =    1752509.1627548405    

Bad direction in the line search;
   refresh the lbfgs memory and restart the iteration.

Process finished with exit code 0

Ps. If the code is hard to read, I'll gladly reformat, just give me some pointers on how to do this.

Yours truly,
Adam Ł.