Pfthb

the link of data from dropboxbadfittingI tried use the curve_fit to fit the data with my pre_defined function in python, but the result was far to perfect. The code is simple and shown as below. I have no idea what's wrong.

Since I am new to python, are there any other optimization or fitting methods which are suitable for my case with predefined function?

Thanks in advance!

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, r1, r2, r3,l,c):
 w=2*math.pi*x
 m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
 n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
 y= (m**2+n**2)**.5
 return y

def readdata(filename):
 x = filename.readlines()
 x = list(map(lambda s: s.strip(), x))
 x = list(map(float, x))
 return x

 # test data
f_x= open(r'C:UsersadmDesktopsimpletryfre.txt')
xdata = readdata(f_x)

f_y= open(r'C:UsersadmDesktopsimpletryimpedance.txt')
ydata = readdata(f_y)

xdata = np.array(xdata)
ydata = np.array(ydata)
plt.semilogx(xdata, ydata, 'b-', label='data')

popt, pcov = curve_fit(func, xdata, ydata, bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf)))
plt.semilogx(xdata, func(xdata, *popt), 'r-', label='fitted curve') 

print(popt)
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

as you guessed, this is a LCR circuit model. now I am trying to fit two curves with the same parameters like

def func1(x, r1, r2, r3,l,c):
w=2*math.pi*x
m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
return m

def func2(x, r1, r2, r3,l,c):
w=2*math.pi*x
n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
return n

is it possible to use the curve_fitting to optimize the parameters?

Here are my results using scipy's differential_evolution genetic algorithm module to generate the initial parameter estimates for curve_fit, along with a simple "brick wall" in the function to ensure all parameters are positive. Scipy's implementation of Differential Evolution uses the Latin Hypercube algorithm to ensure a thorough search of parameter space, which requires bounds within which to search - in this example, those bounds are taken from the data maximum and minimum values. My results:

RMSE: 7.415

R-squared: 0.999995

r1 = 1.16614005e+00

r2 = 2.00000664e+05

r3 = 1.54718886e+01

l = 1.94473531e+04

c = 4.32515535e+05

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

def func(x, r1, r2, r3,l,c):
 # "brick wall" ensuring all parameters are positive
 if r1 < 0.0 or r2 < 0.0 or r3 < 0.0 or l < 0.0 or c < 0.0:
 return 1.0E10 # large value gives large error, curve_fit hits a brick wall

 w=2*numpy.pi*x
 m=r1+(r2*l*w)/(r2**2+l**2*w**2)+r3/(1+r3*c**2*w**2)
 n=(r2**2*l*w)/(r2**2+l**2*w**2)-r3**3*c*w/(1+r3*c**2*w**2)
 y= (m**2+n**2)**.5
 return y


def readdata(filename):
 x = filename.readlines()
 x = list(map(lambda s: s.strip(), x))
 x = list(map(float, x))
 return x

 # test data
f_x= open('/home/zunzun/temp/data/fre.txt')
xData = readdata(f_x)

f_y= open('/home/zunzun/temp/data/impedance.txt')
yData = readdata(f_y)

xData = numpy.array(xData)
yData = numpy.array(yData)


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
 warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
 val = func(xData, *parameterTuple)
 return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
 # min and max used for bounds
 maxX = max(xData)
 minX = min(xData)
 maxY = max(yData)
 minY = min(yData)
 minBound = min(minX, minY)
 maxBound = max(maxX, maxY)
 parameterBounds = 
 parameterBounds.append([minBound, maxBound]) # search bounds for r1
 parameterBounds.append([minBound, maxBound]) # search bounds for r2
 parameterBounds.append([minBound, maxBound]) # search bounds for r3
 parameterBounds.append([minBound, maxBound]) # search bounds for l
 parameterBounds.append([minBound, maxBound]) # search bounds for c

 # "seed" the numpy random number generator for repeatable results
 result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
 return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# now call curve_fit without passing bounds from the genetic algorithm,
# just in case the best fit parameters are aoutside those bounds
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Fitted parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
 f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
 axes = f.add_subplot(111)

 # first the raw data as a scatter plot
 plt.semilogx(xData, yData, 'D')

 # create data for the fitted equation plot
 yModel = func(xData, *fittedParameters)

 # now the model as a line plot
 plt.semilogx(xData, yModel) 

 axes.set_xlabel('X Data') # X axis data label
 axes.set_ylabel('Y Data') # Y axis data label

 plt.show()
 plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)

To have your Least Squares regression make sense, you'll have to at least supply initial parameters which make sense.

As all parameters are by default initiated to the value 1, the biggest influence on the initial regression will be resistor r1 which adds a constant to the mix.

Most probably you'll end up in something like the following configuration:

popt
Out[241]: 
array([1.66581563e+03, 2.43663552e+02, 1.13019744e+00, 1.20233767e+00,
 5.04984535e-04])

Which will output a neat-looking flat line, due to m = something big + ~0 + ~0 ; n=~0 - ~0, so y = r1.

However, if you initialize your parameters somewhat differently,

popt, pcov = curve_fit(func, xdata.flatten(), ydata.flatten(), p0=[0.1,1e5,1000,1000,0.2],
 bounds=((0, 0, 0, 0, 0), (np.inf, np.inf, np.inf, np.inf, np.inf)))

You will get a better looking fit,

popt
Out[244]: 
array([1.14947146e+00, 4.12512324e+05, 1.36182466e+02, 8.29771756e+04,
 1.77593448e+03])

((fitted-ydata.flatten())**2).mean()
Out[257]: 0.6099524982664816
#RMSE hence 0.78

P.s. My data starts at the second data point, due to a conversion error with pd.read_clipboard where the first row became headers instead of data. Shouldn't change the overall picture though.