我试图用一个函数拟合测量数据。然而,即使拟合看起来很好,优化参数(R
和C
)也接近理论值,但这些值的标准差(通过np.sqrt(np.diag(covariance))
计算,如documentation中所述)大约是实际值的10^6倍。我在这里做错了什么?
我已经尝试过使用C * 10^6(因为我认为C = 10 ^-8对于curve_fit
来说太小了),但结果是一样的。
绘图:
输出:
=========================
Parameters:
R = 1005.7461880858972
sigR = 1066772639.8019556
C = 8.825249658942551e-09
sigC = 0.009360746289516244
=========================
代码(降为重要):
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# constants
R = 1000.0 # Ohm
C = 1E-8 # Farad
# [Reading data. f is read, G is calculated]
# Logarithmic Scale & Scatter plot
plt.xscale('log')
plt.scatter(f,G,s=4,color='black')
# curve fitting
## fit function. R and C should be optimized
def opt (fopt, ropt, copt):
return abs( ropt / (ropt + (1 / (1j * fopt * 2 * np.pi * copt ))))
## fit curve: opt is the fit function, f and G are the points the function should fit, p0 are the starting values.
## the fitted R/C values and the covariances are saved in the parms array and the covariance matrix
parms, covariance = curve_fit(opt, f, G, p0=[R,C])
### Print the fitted values for R and C
print('=========================')
print('Parameters:')
print('R = ' + str(parms[0]))
print('sigR = ' + str(np.sqrt(np.diag(covariance))[0]))
print('C = ' + str(parms[1]))
print('sigC = ' + str(np.sqrt(np.diag(covariance))[1]))
print('=========================')
## draw the line: generate an array of values (xopt) and plug them into the opt function
xopt = np.arange(10, 1E6, 10)
gopt = [opt(i, *parms) for i in xopt]
## plot the values calculated before
plt.plot(xopt, gopt, linewidth=0.5, color='black')
代码(完整):
# Import libraries: numpy for advanced math stuff, matplotlib.pyplot for creating fancy plots
import numpy as np
import matplotlib.pyplot as plt
import statistics as stat
from scipy.optimize import curve_fit
# Latex implementation: Process all text with latex
plt.rc('text', usetex=True)
plt.rc('text.latex', preamble=r'\usepackage{siunitx}')
# variables
FILE = '../data/CR.csv'
# constants
R = 1000.0 # Ohm
C = 1E-8 # Farad
# function to decomment a file. it goes through every line and and splits it according to #. Everything before will be passed to var raw.
def decomment(file):
for row in file:
raw = row.split('#')[0].strip()
if raw: yield raw
# Read file and save the data in the 'time' and 'voltage' lists
with open (FILE) as dat:
gen = decomment(dat)
lines = list(gen)
index = [line.split(' ')[0] for line in lines] # index of measurement
UCH1 = [float(line.split(' ')[1]) for line in lines] # Voltage of Channel 1 (reference) [ticks]
UCH2 = [float(line.split(' ')[2]) for line in lines] # Voltage of Channel 2 ('Ausgangsspannung') [ticks]
TS = [float(line.split(' ')[3]) for line in lines] # Time Shift [ticks]
f = [float(line.split(' ')[4]) for line in lines] # Frequency [Hz]
VCH1 = [float(line.split(' ')[5]) for line in lines] # Volts/div of Channel 1 [V]
VCH2 = [float(line.split(' ')[6]) for line in lines] # Volts/div of Channel 2 [V]
Time = [float(line.split(' ')[7]) for line in lines] # Time/div [us]
# scatter plot
UE = np.array(UCH1) * np.array(VCH1)
UA = np.array(UCH2) * np.array(VCH2)
G = UA / UE
# Axis Limits
plt.xlim(10,1E6)
plt.ylim(-0.05,1.05)
# Logarithmic Scale & Scatter plot
plt.xscale('log')
plt.scatter(f,G,s=4,color='black')
# curve fitting
## fit function. R and C should be optimized
def opt (fopt, ropt, copt):
return abs( ropt / (ropt + (1 / (1j * fopt * 2 * np.pi * copt ))))
## fit curve: opt is the fit function, f and G are the points the function should fit, p0 are the starting values.
## the fitted R/C values and the covariances are saved in the parms array and the covariance matrix
parms, covariance = curve_fit(opt, f, G, p0=[R,C])
### Print the fitted values for R and C
print('=========================')
print('Parameters:')
print('R = ' + str(parms[0]))
print('sigR = ' + str(np.sqrt(np.diag(covariance))[0]))
print('C = ' + str(parms[1]))
print('sigC = ' + str(np.sqrt(np.diag(covariance))[1]))
print('=========================')
## draw the line: generate an array of values (xopt) and plug them into the opt function
xopt = np.arange(10, 1E6, 10)
gopt = [opt(i, *parms) for i in xopt]
## plot the values calculated before
plt.plot(xopt, gopt, linewidth=0.5, color='black')
# find fg
## Horizontal line at G = 1 / sqrt(2)
plt.plot(xopt, 1 / np.sqrt(2) + 0 * xopt, linewidth=0.5, linestyle='--', color='black')
## Find the corresponding f value
y_fg = 1
for i in range(len(gopt)):
if abs(gopt[i] - 1/np.sqrt(2)) < y_fg:
y_fg = abs(gopt[i] - 1/np.sqrt(2))
fg_ind = i
fg = xopt[fg_ind]
## Print fg
print('fg = ' + str(fg) + ' Hz')
## Vertical line at f = fg
yval = np.arange(-0.5,1.5,0.01)
plt.plot(fg + yval * 0, yval , linewidth=1, linestyle='--')
## Print the value on the graph
plt.text(fg * 1.1, 0.3, r'$f_\text{g} = \SI{' + str(fg) + r'}{\hertz}$')
plt.text(100, 1 / np.sqrt(2) + 0.03, r'$1 / \sqrt{2}$')
# write fg to file for other scripts
file = open('fg.txt', 'w')
file.write(str(fg))
file.close()
# Axis labels
plt.xlabel(r'$f$ / \si{\hertz}', fontsize=16)
plt.ylabel(r'$|G(f)| = \frac{U_\text{A}}{U_\text{E}}$', fontsize=16)
# Grid
plt.grid(color='gray',which='both',linestyle=':',linewidth=0.1)
# make margins nice and print the plot to a pdf file
plt.tight_layout()
plt.savefig('../plots/CR.pdf')
1条答案
按热度按时间ql3eal8s1#
我自己也遇到过同样的问题,发现问题的根源在于使用了一个冗余的fit参数,在您的例子中,结果证明您的fit函数可以重写(假设ropt和copt是非零数字,您只需从分子和分母中删除
ropt
):啊,但是现在很明显,唯一重要的是
ropt * copt
,所以它们的最佳值不是唯一定义的!解决方案是重写fit,使用单个参数:当然,这意味着你的数据不能告诉你
R
和C
的值,只能告诉你它们的乘积RC
。