3

I'm new to python (and programming in general) and want to make a polynomial fit using curve_fit, where the order of the polynomials (or the number of fit parameters) is variable.

I made this code which is working for a fixed number of 3 parameters a,b,c

# fit function
def fit_func(x, a,b,c):
    p = np.polyval([a,b,c], x)
    return p

# do the fitting
popt, pcov = curve_fit(fit_func, x_data, y_data)

But now I'd like to have my fit function to only depend on a number N of parameters instead of a,b,c,....

I'm guessing that's not a very hard thing to do, but because of my limited knowledge I can't get it work.

I've already looked at this question, but I wasn't able to apply it to my problem.

scavi
  • 57
  • 7
  • I recommend looking at the numpy documentation https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html which has example code at the bottom of the page. – James Phillips Jun 06 '18 at 19:16
  • 1
    But that's the documentation for `polyfit`. I need to use `curve_fit` because I have to expand the code for other types of functions later on. – scavi Jun 06 '18 at 19:21

2 Answers2

7

You can define the function to be fit to your data like this:

def fit_func(x, *coeffs):
    y = np.polyval(coeffs, x)
    return y

Then, when you call curve_fit, set the argument p0 to the initial guess of the polynomial coefficients. For example, this plot is generated by the script that follows.

plot

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


# Generate a sample input dataset for the demonstration.
x = np.arange(12)
y = np.cos(0.4*x)


def fit_func(x, *coeffs):
    y = np.polyval(coeffs, x)
    return y


fit_results = []
for n in range(2, 6):
    # The initial guess of the parameters to be found by curve_fit.
    # Warning: in general, an array of ones might not be a good enough
    # guess for `curve_fit`, but in this example, it works.
    p0 = np.ones(n)

    popt, pcov = curve_fit(fit_func, x, y, p0=p0)
    # XXX Should check pcov here, but in this example, curve_fit converges.

    fit_results.append(popt)


plt.plot(x, y, 'k.', label='data')

xx = np.linspace(x.min(), x.max(), 100)
for p in fit_results:
    yy = fit_func(xx, *p)
    plt.plot(xx, yy, alpha=0.6, label='n = %d' % len(p))

plt.legend(framealpha=1, shadow=True)
plt.grid(True)
plt.xlabel('x')
plt.show()
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
0

The parameters of polyval specify p is an array of coefficients from the highest to lowest. With x being a number or array of numbers to evaluate the polynomial at. It says, the following.

If p is of length N, this function returns the value:

p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]

def fit_func(p,x):
    z = np.polyval(p,x)
    return z

e.g.

t= np.array([3,4,5,3])
y = fit_func(t,5)
503

which is if you do the math here is right.