0

这是我的代码，这是我的数据，这是代码的输出。我尝试在 x 轴上添加一个值，认为也许值这么少可以解释为零。我不知道 true_divide 可能是什么，我无法解释这个除以零错误，因为我的数据中没有一个零，检查了我所有的 2500 个数据点。希望你们中的一些人可以提供一些澄清。提前致谢。

import pandas as pd
import matplotlib.pyplot as plt
from iminuit import cost, Minuit
import numpy as np

frame = pd.read_excel('/Users/lorenzotecchia/Desktop/Analisi Laboratorio/Analisi     dati/Quinta Esperienza/500Hz/F0000CH2.xlsx', 'F0000CH2')
data =  pd.read_excel('/Users/lorenzotecchia/Desktop/Analisi Laboratorio/Analisi     dati/Quinta Esperienza/500Hz/F0000CH1.xlsx', 'F0000CH1')
# tempi_500Hz = pd.DataFrame(frame,columns=['x'])
# Vout_500Hz = pd.DataFrame(frame,columns=['y'])
tempi_500Hz = pd.DataFrame(frame,columns=['x1'])
Vout_500Hz = pd.DataFrame(frame,columns=['y1']) 
# Vin_500Hz = pd.DataFrame(data,columns=['y'])

def fit_esponenziale(x, α, β):
    return α * (1 - np.exp(-x / β))

plt.xlabel('ω(Hz)')
plt.ylabel('Attenuazioni')
plt.title('Fit Parabolico')
plt.scatter(tempi_500Hz, Vout_500Hz)
least_squares = cost.LeastSquares(tempi_500Hz, Vout_500Hz, np.sqrt(Vout_500Hz), fit_esponenziale)
m = Minuit(least_squares, α=0, β=0)
m.migrad()
m.hesse()
plt.errorbar(tempi_500Hz, Vout_500Hz, fmt="o", label="data")
plt.plot(tempi_500Hz, fit_esponenziale(tempi_500Hz, *m.values), label="fit")
fit_info = [
f"$\\chi^2$ / $n_\\mathrm{{dof}}$ = {m.fval:.1f} / {len(tempi_500Hz) - m.nfit}",]

for p, v, e in zip(m.parameters, m.values, m.errors):
    fit_info.append(f"{p} = ${v:.3f} \\pm {e:.3f}$")

plt.legend()
plt.show()

输入

输出和数据示例

4

2 回答 2

0

这是一个工作Minuit与curve_fit示例。我对函数进行了缩放，使得指数的衰减约为 1（对于非线性拟合通常是一个好主意）。最终，这两种方法都给出了非常相似的结果。

注意：无论错误是否有意义，我都将其保持打开状态。等于零的起始值绝对是个坏主意。

import pandas as pd
import matplotlib.pyplot as plt
from iminuit import cost, Minuit
from scipy.optimize import curve_fit
import numpy as np

def fit_exp(x, a, b, c):
    return a * (1 - np.exp(- 1000 * b * x) ) + c

nn = 170
xl = np.linspace( 0, 0.001, nn )
yl = fit_exp( xl, 15, 5.3, -8.1 ) + np.random.normal( size=nn, scale=0.05 )

#######################
### Minuit
#######################
least_squares = cost.LeastSquares(xl, yl, np.sqrt( np.abs( yl ) ), fit_exp )
print(least_squares)
m = Minuit(least_squares, a=1, b=5, c=-7)
print( "grad: ")
print( m.migrad() ) ### needs to be called to get fit values
print( m.values )### gives slightly different output
print("Hesse:")
print( m.hesse() )

#######################
### curve_fit
#######################

opt, cov = curve_fit(
    fit_exp, xl, yl, sigma=np.sqrt( np.abs( yl ) ),
    absolute_sigma=True
)
print( " curve_fit: ")
print( opt )
print( " covariance ")
print( cov )

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1 )
ax.plot(xl, yl, marker ='+', ls='' )

ax.plot(xl, fit_exp(xl, *m.values), ls="--")
ax.plot(xl, fit_exp(xl, *opt), ls=":")


plt.show()
于 2021-12-06T06:51:36.297 回答
0

实际上有一种方法可以完全线性拟合。参见例如这里

import matplotlib.pyplot as plt
import numpy as np

from scipy.integrate import cumtrapz



def fit_exp(x, a, b, c):
    return a * (1 - np.exp( -b * x) ) + c

nn = 170
xl = np.linspace( 0, 0.001, nn )
yl = fit_exp( xl, 15, 5300, -8.1 ) + np.random.normal( size=nn, scale=0.05 )

"""
with y = a( 1- exp(-bx) ) + c
we have Y = int y = -1/b y + d x + h ....try it out or see below
so we get a linear equation for b (actually 1/b ) to optimize
this goes as:
"""

Yl = cumtrapz( yl, xl, initial=0 )
ST = [xl, yl, np.ones( nn ) ]
S = np.transpose( ST )

eta = np.dot( ST, Yl )
A = np.dot( ST, S )

sol = np.linalg.solve( A, eta )
bFit = -1/sol[1]
print( bFit )
"""
now we can do a normal linear fit
"""

ST = [ fit_exp(xl, 1, bFit, 0), np.ones( nn ) ]
S = np.transpose( ST )
A = np.dot( ST, S )
eta = np.dot( ST, yl )
aFit, cFit = np.linalg.solve( A, eta )

print( aFit, cFit )
print(aFit + cFit, sol[0] ) ### consistency check

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1 )
ax.plot(xl, yl, marker ='+', ls='' )
## at best a sufficient fit, at worst a good start for a non-linear fit
ax.plot(xl, fit_exp( xl, aFit, bFit, cFit ) ) 

plt.show()

"""
a ( 1 - exp(-b x)) + c = a + c - a exp(-b x) = d - a exp( -b x )
int y = d x + a/b exp( -b x ) + g
    = d x +a/b exp( -b x ) + a/b - a/b + c/b - c/b + g
    = d x - 1/b ( a - a exp( -b x ) + c ) + c/b + a/b + g
    = d x + k y + h
with k = -1/b and h = g + c/b + a/b.
d and h are fitted but not used, but as a+c = d we can check 
for consistency
"""
于 2021-12-06T12:29:59.430 回答