Search code examples
pythonscipycurve-fittingscientific-computingfunction-fitting

Curve fitting with python error


I'm trying to fit my data to (cos(x))^n. The vale of n in theory is 2, but my data should give me around 1.7. When I define my fitting function and I try curve_fit, I get an error

def f(x,a,b,c):
   return a+b*np.power(np.cos(x),c)

param, extras = curve_fit(f, x, y)

This is my data

x   y               error
90  3.3888756187    1.8408898986
60  2.7662844365    1.6632150903
45  2.137309503     1.4619540017
30  1.5256883339    1.2351875703
0   1.4665463518    1.2110104672

The error looks like this:

/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in power after removing the cwd from sys.path.

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py:690: OptimizeWarning: Covariance of the parameters could not be estimated
category=OptimizeWarning)


Solution

  • The problem is that cos(x) can become negative and then cos(x) ^ n can be undefined. Illustration:

    np.cos(90)
    -0.44807361612917013
    

    and e.g.

    np.cos(90) ** 1.7
    nan
    

    That causes the two error messages you receive.

    It works fine, if you modify your model, e.g. to a + b * np.cos(c * x + d). Then the plot looks as follows:

    enter image description here

    The code can be found below with some inline comments:

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.optimize import curve_fit
    
    
    def f(x, a, b, c, d):
    
        return a + b * np.cos(c * x + d)
    
    # your data
    xdata = [90, 60, 45, 30, 0]
    ydata = [3.3888756187, 2.7662844365, 2.137309503, 1.5256883339, 1.4665463518]
    
    # plot data
    plt.plot(xdata, ydata, 'bo', label='data')
    
    # fit the data
    popt, pcov = curve_fit(f, xdata, ydata, p0=[3., .5, 0.1, 10.])
    
    # plot the result
    xdata_new = np.linspace(0, 100, 200)
    plt.plot(xdata_new, f(xdata_new, *popt), 'r-', label='fit')
    plt.legend(loc='best')
    plt.show()