Fitting data with a strange function in python

Question

I am fitting data with a strange function,I have to solve a equation to get one term. there is the code:

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

I've got a good fit ,but there are two warning about this fit.

E:/03_07/ParPhy/fit_code/TEST_9.py:37: RuntimeWarning: invalid value encountered in log 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: The iteration is not making good progress, as measured by the improvement from the last ten iterations. warnings.warn(msg, RuntimeWarning)

The curve_fit input some 'bad'parameters to my function sometimes , so there is no root of my equation. In fact x(s_p) must be <0.

and when I change the initial guess of fsolve ,the output parameters change too,that's so weird!

Any help will be great appreciated!

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)

I set l1 -100000000 , and I use these output a,b,c to plot the f(x) at[-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()

What I find is that the root of my equation is not in the domain of [-10,0] which @Mstaino sets.

It seems like that I have to set l1 large enough?

Answer 1

curve_fit is essentially a least squares fit. The problem when you have these algorithms a lot of times arises from the lack of boundaries. In your case, I think the problem is caused by fsolve, which in turn is caused by the lack of boundaries in your curve_fit, which cause some of the f(x) to be unsolvable.

I managed to eliminate the warnings by eliminating the list in f(x) and changing fsolve to brentq with appropriate limits, since your x domain in f is clearly negative.

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)

If you uncomment the print line you will see some values of parameters which are attempted by curve_fit which give an f that cannot be solved. If you can bound those, you help ensure a good matching, as an alternative (or better yet, an addition) to putting a good "initial value"