python - 曲线拟合方程python

Question

我正在尝试使用具有 6 个组件的 rlc 电路的导纳方程拟合数据。我正在按照给定he[fit] 1 re 的示例并插入我的方程式。该方程是使用 Mathcad 简化的 6 分量电路的导纳的实部。在附图中，x 轴是欧米茄 (w=2*pi*f)，y 是导纳，以毫西门子为单位。

该程序运行，但尽管有良好的试用功能，但它并没有进行拟合。我感谢任何帮助，为什么适合是一条直线。我还附上了一个高斯拟合示例。

这就是我尝试拟合方程式时得到的结果。数据是左侧峰值较小的数据，试验函数是虚线。拟合是一条直线

from numpy import sqrt, pi, exp, linspace, loadtxt
from lmfit import  Model

import matplotlib.pyplot as plt

data = loadtxt("C:/Users/susu/circuit_eq_real5.dat")
x = data[:, 0]
y = data[:, 1]

def circuit(x,C0,Cm,Lm,Rm,R0,Rs):


    return ((C0**2*Cm**2*Lm**2*R0*x**4)+(Rs*C0**2*Cm**2*Lm**2*x**4)+(C0**2*Cm**2*R0**2*Rm*x**2)+(Rs*C0**2*Cm**2*R0**2*x**2)+(C0**2*Cm**2*R0*Rm**2*x**2)+(2*Rs*C0**2*Cm**2*R0*Rm*x**2)+(Rs*C0**2*Cm**2*Rm**2*x**2)-(2*C0**2*Cm*Lm*R0*x**2)-(2*Rs*C0**2*Cm*Lm*x**2)+(C0**2*R0)+(Rs*C0**2)-(2*Rs*C0*Cm**2*Lm*x**2)+(2*Rs*C0*Cm)+(Cm**2*Rm)+(Rs*Cm**2))/((C0**2*Cm**2*Lm**2*x**4)+(C0**2*Cm**2*R0**2*x**2)+(2*C0**2*Cm**2*R0*Rm*x**2)+(C0**2*Cm**2*Rm**2*x**2)-(2*C0**2*Cm*Lm*x**2)+(C0**2)-(2*C0*Cm**2*Lm*x**2)+(2*C0*Cm)+(Cm**2))

gmodel = Model(circuit)
result = gmodel.fit(y, x=x, C0=1.0408*10**(-12), Cm=5.953*10**(-14), 
Lm=1.475*10**(-7), Rm=1.571, R0=2.44088, Rs=0.42)

print(result.fit_report())

plt.plot(x, y,         'bo')
plt.plot(x, result.init_fit, 'k--')
plt.plot(x, result.best_fit, 'r-')
plt.show()

以下是合身报告

 [[Fit Statistics]]
# fitting method   = leastsq
# function evals   = 14005
# data points      = 237
# variables        = 6
chi-square         = 32134074.5
reduced chi-square = 139108.548
Akaike info crit   = 2812.71607
Bayesian info crit = 2833.52443
[[Variables]]
C0: -7.5344e-15 +/- 6.3081e-09 (83723736.65%) (init = 1.0408e-12)
Cm: -8.9529e-13 +/- 1.4518e-06 (162164237.47%) (init = 5.953e-14)
Lm:  2.4263e-06 +/- 1.94051104 (79978205.20%) (init = 1.475e-07)
Rm: -557.974399 +/- 1.3689e+09 (245334051.75%) (init = 1.571)
R0: -5178.53517 +/- 6.7885e+08 (13108904.45%) (init = 2.44088)
Rs:  2697.67659 +/- 7.3197e+08 (27133477.70%) (init = 0.42)
[[Correlations]] (unreported correlations are < 0.100)
C(R0, Rs) = -1.003
C(Rm, Rs) = -0.987
C(Rm, R0) =  0.973
C(C0, Lm) =  0.952
C(C0, Cm) = -0.502
C(Cm, R0) = -0.483
C(Cm, Rs) =  0.453
C(Cm, Rm) = -0.388
C(Cm, Lm) = -0.349
C(C0, R0) =  0.310
C(C0, Rs) = -0.248
C(C0, Rm) =  0.148

非常感谢 M Newville 和 Mikuszefski 以及其他人的见解和反馈。我同意我放在那里的东西可能会让程序一团糟。从 Python 代码中可以看出我并不精通 Python 或编程。

Mikuszefsky，感谢您发布 rlc 示例代码。你的方法简洁有趣。我不知道 Python 会直接进行复杂的拟合。我会尝试您的方法，看看是否可以进行拟合。我想同时拟合 Y（导纳）的实部和虚部。我肯定会被困在某个地方，并会在这里发布我的进度。最好的，苏苏

Answer 1

这是一种清理具有并联和串联连接的 RLC 电路的方法。这避免了这条超长的线路和难以检查的功能。它还避免使用 Matlab 或类似程序，因为它直接计算电路。当然，它可以很容易地扩展 OP 的电路。正如 M Newville 所指出的，简单的拟合失败了。另一方面，如果将单位缩放为自然单位，即使没有初始参数，它也可以工作。请注意，结果仅通过比例因子是正确的。需要知道至少一个组件的值。

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

def r_l( w, l ):
    return 1j * w * l

def r_c( w, c ):
    return 1. / ( 1j * w * c )

def parallel( a, b ):
    return( 1. / ( 1./ a + 1. / b ) )

def series( a, b ):
    return a + b

# simple rlc band pass filter (to be extended)
def rlc_band( w , r, l, c ):
    lc = parallel( r_c( w , c ), r_l( w, l ) )
    return lc / series( r, lc )

def rlc_band_real( w , r, l, c ):
    return rlc_band( w , r, l, c ).real

def rlc_band_real_milli_nano( w , r, l, c ):
    return rlc_band_real( w , r, 1e-6 * l, 1e-9 * c ).real

wList = np.logspace( 5, 7, 25 )
wFullList = np.logspace( 5, 7, 500 )
rComplexList = np.fromiter( ( rlc_band(w, 12, 1.3e-5, 1e-7 ) for w in wList ), np.complex )
rList = np.fromiter( ( r.real for r in rComplexList ), np.float )
pList = np.fromiter( ( np.angle( r ) for r in rComplexList ), np.float )

fit1, pcov = curve_fit( rlc_band_real, wList, rList )
print fit1
print "does not work"

fit2, pcov = curve_fit( rlc_band_real_milli_nano, wList, rList )
print fit2
print "works, but is not unique (scaling is possible)"
print 12, fit2[1] * 12 / fit2[0], fit2[2] * fit2[0] / 12.

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )

ax.plot( wList, rList , ls='', marker='o', label='data')
#~ ax.plot( wList, pList )
ax.plot( wFullList, [ rlc_band_real( w , *fit1 ) for w in wFullList ], label='naive fit')
ax.plot( wFullList, [ rlc_band_real_milli_nano( w , *fit2 ) for w in wFullList ], label='scaled units')
ax.set_xscale('log')
ax.legend( loc=0 )

plt.show()

提供：

>> /...minpack.py:785: OptimizeWarning: Covariance of the parameters could not be estimated category=OptimizeWarning)
>> [1. 1. 1.]
>> does not work
>> [ 98.869924   107.10908434  12.13715912]
>> works, but is not unique (scaling is possible)
>> 12 13.0 100.0

Answer 2

Providing a real link to a text file of the actual data you are using and/or a real plot of what you are actually seeing would be most helpful. Also, please provide an accurate and complete description of the results including the text of what is actually printed out by the print(result.fit_report()). Basically, ask yourself how you might try to help someone who asked such a question, and provide as much information as you can.

No one (including you) is ever going to be able to spell-check the implementation of your function. You will need thorough and robust testing of this function in order to convince anyone (including you, I hope) that it is doing what you think it should do. You should provide the results of those tests before worrying about why it is not working as a fitting function. You should definitely consider refactoring that mess of an equation into more manageable and readable pieces.

That said, I also strongly recommend that you do not work in units of Farads and Henrys but picoFarads or nanoFarads and microHenrys. That will make the values much closer to 1 (say, order 1e-6 to 1e+6), which will make it much easier for the fit to do its job.