Подгонка данных со странной функцией в Python

Я подбираю данные с помощью странной функции, мне нужно решить уравнение, чтобы получить один член. вот код:

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 \ optimize \ minpack.py: 161: RuntimeWarning: итерация не продвигается хорошо, если судить по улучшению по сравнению с последними десятью итерациями. warnings.warn (msg, RuntimeWarning)

Curve_fit иногда вводит некоторые «плохие» параметры в мою функцию, поэтому в моем уравнении нет корня. Фактически x (s_p) должен быть ‹0.

и когда я изменяю первоначальное предположение 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 достаточно большим?


python-3.x scipy model-fitting curve-fitting
person Gene Y    schedule 11.03.2019    source источник
comment
Если вы знаете, что можете безопасно игнорировать предупреждения, вы можете импортировать предупреждения, а затем вызвать warnings.filterwarnings ('ignore') перед подгонкой кривой - это подавит все предупреждения. Затем вы можете вызвать warnings.filterwarnings («по умолчанию») после подбора кривой, чтобы снова восстановить предупреждения.   -  person James Phillips    schedule 11.03.2019

Ответы (1)


arrow_upward
1
arrow_downward

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

person Mstaino    schedule 11.03.2019
comment
Это отлично сработало для меня и, очевидно, лучше, чем мое (постыдно неубедительное) предложение отключить предупреждения. - person James Phillips; 11.03.2019
comment
Я отредактировал свой вопрос, и есть еще одна проблема, и мне интересно, почему вы устанавливаете s_p на l2 (такое маленькое число), когда f (l1) и f (l2) имеют одинаковый знак, я знаю, что curve_fit () хочет найти минимум хи-квадрат. Итак, мой вопрос: не удастся ли s_p = l2 подойти? @Mstaino - person Gene Y; 12.03.2019
comment
Я установил его произвольно, так как видел, что функция убывает с x и в основном положительна, и что для некоторых значений коэффициентов нет 0, а f выравнивается. Решение состоит в том, чтобы вообще избежать этих значений коэффициентов (но в этом случае вам нужно связать их в curve_fit), альтернативы позволяют небольшое значение x, чтобы f было как можно ближе к 0 (что я и сделал). Если вы удалите if, brentq вызовет ValueError, в отличие от fsolve, который вернет начальное число (-2), что также неверно, и выдаст вам предупреждающее сообщение о том, что хорошее решение не найдено. - person Mstaino; 12.03.2019