在python中用一个奇怪的函数拟合数据

Question

我正在用一个奇怪的函数拟合数据，我必须解一个方程才能得到一个项。 有代码：

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


y=np.array([6.9, 7.4, 8.2, 8.7, 9.2, 9.8, 10.4, 11.0, 11.8, 12.5, 13.2, 13.6, 14.7, 16.0, 17.7, 17.8, 18.8, 20.6, 22.5, 24.9, 26.2, 26.6, 29.3, 32.0, 33.2, 33.6, 36.2, 38.2, 39.3, 40.9, 41.8, 43.6, 44.8, 45.0, 45.2, 43.7, 37.1, 30.3, 30.6, 30.1, 27.7, 25.9, 25.5, 24.1, 22.9, 21.4, 19.8, 18.1, 16.4, 16.0, 15.0, 14.2, 13.0, 12.1, 11.0, 10.3, 10.0, 9.9, 8.7, 7.9])
s=np.array([0.36300625, 0.36905625, 0.37515625, 0.38130625, 0.38750625, 0.39375625, 0.40005625, 0.40640625, 0.41280625, 0.41925625, 0.42575625, 0.43230625, 0.43890625, 0.44555625, 0.45225625, 0.45900625, 0.46580625, 0.47265625, 0.47955625, 0.48650625, 0.49350625, 0.50055625, 0.50765625, 0.51480625, 0.52200625, 0.52925625, 0.53655625, 0.54390625, 0.55130625, 0.55875625, 0.56625625, 0.57380625, 0.58140625, 0.58905625, 0.59675625, 0.60450625, 0.61230625, 0.62015625, 0.62805625, 0.63600625, 0.64400625, 0.65205625, 0.66015625, 0.66830625, 0.67650625, 0.68475625, 0.69305625, 0.70140625, 0.70980625, 0.71825625, 0.72675625, 0.73530625, 0.74390625, 0.75255625, 0.76125625, 0.77000625, 0.77880625, 0.78765625, 0.79655625, 0.80550625])
err_y=np.array([0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6, 0.6, 0.6, 0.6, 0.6, 0.7, 0.7, 0.7, 0.7, 0.8, 0.8, 0.9, 0.9, 0.9, 0.9, 0.9, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.9, 0.8, 0.8, 0.8, 0.8, 0.7, 0.8, 0.7, 0.7, 0.7, 0.7, 0.6, 0.6, 0.6, 0.6, 0.6, 0.6, 0.5, 0.5, 0.5, 0.5, 0.5, 0.4, 0.4])

m_pi=139.57/1000
m_omega=782.65/1000
gamma_omega=8.49/1000


def F_s(s,a,b,c,alpha,kappa):
    u=(1-4*m_pi**2/s)**(1/2)
    g=-1/np.pi*u*np.log((1+u)/(1-u))+1j*u

    def f(x):
        x = float(x)
        u_1=(1-4*m_pi**2/x)**(1/2)
        return [
                a*x**2+b*x+c+(x-4*m_pi**2)/(4*np.pi)*u_1*np.log((1+u_1)/(u_1-1))
                ]
    s_p= fsolve(f, -2)[0].item()

    A=(c+m_pi**2*(-2/np.pi))*(s_p-s)
    B=s_p*(a*s**2+b*s+c-(s-4*m_pi**2)*g/4)
    P_s=1+alpha*s+kappa*s/(m_omega**2-s-1j*m_omega*gamma_omega)
    return abs(P_s*A/B)**2


popt, pcov = curve_fit(F_s, s, y,p0=(-1.43,0.24,0.22,0.083,0.0018),maxfev=50000)



def chi_square():
    chi_sq=0
    for i in range(len(err_y)):
        chi_sq=chi_sq+ (F_s(s[i],*popt)-y[i])**2/err_y[i]**2
    return chi_sq


yvals=F_s(s,*popt)


print('parameters:',popt,'\n','chi-square:',chi_square(),'\n','dof:',len(y),'-',len(popt),'\n','chi-square/dof:',chi_square()/(len(y)-len(popt)))

plt.errorbar(s**(1/2), y, yerr=err_y, fmt='.', color='black', ecolor='black', elinewidth=2, capsize=0,label='BESIII_data')
plt.plot(s**(1/2), yvals, 'r',label='curve_fit values')

我很合身，但有两个关于这个合身的警告。

E:/03_07/ParPhy/fit_code/TEST_9.py:37：RuntimeWarning：日志中遇到无效值 ax2+bx+c+(x-4*m_pi**2)/(4*np.pi)*u_1*np.log((1+u_1)/( u_1-1))

E:\anaconda\lib\site-packages\scipy\optimize\minpack.py:161: RuntimeWarning: 迭代没有取得好的进展，由 过去十次迭代的改进。 warnings.warn(msg, RuntimeWarning)

curve_fit 有时会向我的函数输入一些“坏”参数，因此我的方程没有根。实际上 x(s_p) 必须是

当我改变 fsolve 的初始猜测时，输出参数也会改变，太奇怪了！

任何帮助将不胜感激！

def f(x):
    x = float(x)
    u_1=(1-4*m_pi**2/x)**(1/2)
    return a*x**2+b*x+c+(x-4*m_pi**2)/(4*np.pi)*u_1*np.log((1+u_1)/(u_1-1))

# s_p= fsolve(f, -2)[0].item()
l1, l2 = -1000000000, -1e-10
if f(l1) * f(l2) >= 0:   #brentq will raise a value error if endpoints do not have the same sign
    s_p = l2
    # print(f(l2), a, b, c, alpha, kappa)
else:
    s_p = brentq(f, l1, l2)

我设置 l1 -100000000 ，并使用这些输出 a,b,c 在 [-2000,0] 处绘制 f(x)

import matplotlib.pyplot as plt
import numpy as np

m_pi=139.57/1000
m_omega=781.94/1000
gamma_omega=8.43/1000

a=-2.62498839e-04
b=-1.39411562e+00
c=7.08577488e-01


def test_chao(x):
    u=(1-4*m_pi**2/x)**(1/2)
    return a*x**2+b*x+c+(x-4*m_pi**2)/(4*np.pi)*u*np.log((1+u)/(u-1))

x=np.linspace(-2000,0)
yval=test_chao(x)
plot1=plt.plot(x, yval, '*',label='original values')

plt.xlabel('s axis')
plt.ylabel('y axis')
plt.legend(loc=3)
plt.title('curve_fit')
plt.show()

我发现我的方程的根不在@Mstaino 设置的 [-10,0] 的域中。

看来我必须将 l1 设置得足够大？

Answer 1

curve_fit 本质上是最小二乘拟合。当你有很多次使用这些算法时，问题来自于缺乏边界。在您的情况下，我认为问题是由fsolve 引起的，而这又是由于您的curve_fit 中缺少边界引起的，这导致某些f(x) 无法解决。

我设法通过消除 f(x) 中的列表并将 fsolve 更改为具有适当限制的 brentq 来消除警告，因为 f 中的 x 域显然是负面的。

from scipy.optimize import fsolve, brentq

#...rest of code

def f(x):
    x = float(x)
    u_1=(1-4*m_pi**2/x)**(1/2)
    return a*x**2+b*x+c+(x-4*m_pi**2)/(4*np.pi)*u_1*np.log((1+u_1)/(u_1-1))

# s_p= fsolve(f, -2)[0].item()
l1, l2 = -10, -1e-10
if f(l1) * f(l2) >= 0:   #brentq will raise a value error if endpoints do not have the same sign
    s_p = l2
    # print(f(l2), a, b, c, alpha, kappa)
else:
    s_p = brentq(f, l1, l2)

如果您取消注释print 行，您将看到curve_fit 尝试的一些参数值，这些值给出了无法解决的f。如果您可以绑定这些，则有助于确保良好的匹配，作为放置良好“初始值”的替代方案（或更好的补充）