Подгонка данных со странной функцией в 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: в журнале обнаружено недопустимое значение a x 2 + b x + 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 \ optimpack \ minpack.py: 161: RuntimeWarning: Итерацияне добился хорошего прогресса, что измеряется улучшением за последние десять итераций.warnings.warn (msg, RuntimeWarning)

Curve_fit иногда вводит некоторые «плохие» параметры в мою функцию, поэтому корень моего уравнения отсутствует.На самом деле x (s_p) должно быть <0. </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, чтобы построить f (x) на [-2000,0]

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()

Я обнаружил, что корень моего уравнения не находится в области [-10,0], которую устанавливает @Mstaino.

Кажется, я должен установить l1 достаточно большим?

Answer 1

curve_fit по сути подходит для наименьших квадратов.Проблема, когда у вас есть эти алгоритмы много раз возникает из-за отсутствия границ.В вашем случае, я думаю, что проблема вызвана fsolve, что, в свою очередь, вызвано отсутствием границ в вашем Curve_fit, что приводит к неразрешимости некоторых из f(x).

Мне удалосьудалите предупреждения, исключив список в f (x) и изменив fsolve на brentq с соответствующими ограничениями, поскольку ваш домен x в f явно отрицательный.

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, который не может быть решен.Если вы можете связать их, вы поможете обеспечить хорошее соответствие в качестве альтернативы (или, что еще лучше, дополнения), чтобы поставить хорошее «начальное значение»