scipy.optimize.curve_fit得到巨大的标准差/协方差

klsxnrf1  于 2022-12-23  发布在  其他
关注(0)|答案(1)|浏览(165)

我试图用一个函数拟合测量数据。然而,即使拟合看起来很好,优化参数(RC)也接近理论值,但这些值的标准差(通过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')
ql3eal8s

ql3eal8s1#

我自己也遇到过同样的问题,发现问题的根源在于使用了一个冗余的fit参数,在您的例子中,结果证明您的fit函数可以重写(假设ropt和copt是非零数字,您只需从分子和分母中删除ropt):

def opt (fopt, ropt, copt):
    return abs( 1 / (1 + (1 / (1j * fopt * 2 * np.pi * copt * ropt))))

啊,但是现在很明显,唯一重要的是ropt * copt,所以它们的最佳值不是唯一定义的!解决方案是重写fit,使用单个参数:

def opt (fopt, rcopt):
    return abs( 1 / (1 + (1 / (1j * fopt * 2 * np.pi * rcopt))))

当然,这意味着你的数据不能告诉你RC的值,只能告诉你它们的乘积RC

相关问题