Pass tuple as input argument for scipy.optimize.curve_fit

Question

I have the following code:

import numpy as np
from scipy.optimize import curve_fit


def func(x, p): return p[0] + p[1] + x


popt, pcov = curve_fit(func, np.arange(10), np.arange(10), p0=(0, 0)) 

It will raise TypeError: func() takes exactly 2 arguments (3 given). Well, that sounds fair - curve_fit unpact the (0, 0) to be two scalar inputs. So I tried this:

popt, pcov = curve_fit(func, np.arange(10), np.arange(10), p0=((0, 0),))

Again, it said: ValueError: object too deep for desired array

If I left it as default (not specifying p0):

popt, pcov = curve_fit(func, np.arange(10), np.arange(10))  

It will raise IndexError: invalid index to scalar variable. Obviously, it only gave the function a scalar for p.

I can make def func(x, p1, p2): return p1 + p2 + x to get it working, but with more complicated situations the code is going to look verbose and messy. I'd really love it if there's a cleaner solution to this problem.

Thanks!

Answer 1

This is an old thread now, but I also just ran into this issue. Building on Emile Maras' solution, but expanding the function to to return either the nth order polynomial fitting function for curve_fit, or the y values based on fit results. This facilitates plotting and residual calculations. Here is an example that fits data to progressively higher order polynomials and plots the results and residuals.

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

def funToFit(x):
    return 0.5 + 2*x -3*x**2 + 0.2*x**3 + 0.1*x**4

x=np.arange(30)
y=funToFit(x)+np.random.normal(0,5000,x.size)

def polyfun(order=0,x=np.arange(0),P=np.arange(0)):
    if P.size>0 and x.size>0:
        y=0
        for i in range(P.size):
            y+=P[i]*np.power(x,i)
        return y
    elif order>0:
        def fitfun(x,*a):
            y=0
            for i in range(order+1):
                y+=a[i]*np.power(x,i)
            return y
        return fitfun
    else:
        raise Exception("Either order or x and P must be provided")

plt.figure("fits")
plt.plot(x,y,color="black")
for i in range(4):
    order = i+1
    [fit,covar] = curve_fit(polyfun(order=order),x,y,p0=(1,)*(order+1))
    yfit = polyfun(x=x,P=fit)
    res = yfit-y
    plt.figure("fits")  
    plt.plot(x,yfit)
    plt.figure("res")
    plt.plot(x,res)

Answer 2

Problem

When using curve_fit you must explicitly say the number of fit parameters. Doing something like:

def f(x, *p):
    return sum( [p[i]*x**i for i in range(len(p))] )

would be great, since it would be a general nth-order polynomial fitting function, but unfortunately, in my SciPy 0.12.0, it raises:

ValueError: Unable to determine number of fit parameters.

Solution

So you should do:

def f_1(x, p0, p1):
    return p0 + p1*x

def f_2(x, p0, p1, p2):
    return p0 + p1*x + p2*x**2

and so forth...

Then you can call using the p0 argument:

curve_fit(f_1, xdata, ydata, p0=(0,0))

Answer 3

You can define functions that return other functions (see Passing additional arguments using scipy.optimize.curve_fit? )

Working example :

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

def funToFit(x):
    return 0.5+2*x-3*x*x+0.2*x*x*x+0.1*x*x*x*x


xx=[random.uniform(1,5) for i in range(30)]
yy=[funToFit(xx[i])+random.uniform(-1,1) for i in range(len(xx))]


a=np.zeros(5)
def make_func(numarg):
    def func(x,*a):
        ng=numarg
        v=0
        for i in range(ng):
            v+=a[i]*np.power(x,i)
        return v
    return func

leastsq, covar = curve_fit(make_func(len(a)),xx,yy,tuple(a))
print leastsq
def fFited(x):
    v=0
    for i in range(len(leastsq)):
        v+=leastsq[i]*np.power(x,i)
    return v


xfine=np.linspace(1,5,200)
plt.plot(xx,yy,".")
plt.plot(xfine,fFited(xfine))
plt.show()

Answer 4

scipy.optimize.curve_fit

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, **kw)

Use non-linear least squares to fit a function, f, to data.

Assumes ydata = f(xdata, *params) + eps

Explaining the idea

The function to be fitted should take only scalars (not: *p0). Remember that the result of the fit depends on the initialization parameters.

Working example

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

def func(x, a0, a1):
    return a0 + a1 * x

x, y = np.arange(10), np.arange(10) + np.random.randn(10)/10
popt, pcov = curve_fit(func, x, y, p0=(1, 1))

# Plot the results
plt.title('Fit parameters:\n a0=%.2e a1=%.2e' % (popt[0], popt[1]))
# Data
plt.plot(x, y, 'rx')
# Fitted function
x_fine = np.linspace(x[0], x[-1], 100)
plt.plot(x_fine, func(x_fine, popt[0], popt[1]), 'b-')
plt.savefig('Linear_fit.png')
plt.show()

Answer 5

Not sure if this is cleaner, but at least it is easier now to add more parameters to the fitting function. Maybe one could even make an even better solution out of this.

import numpy as np
from scipy.optimize import curve_fit


def func(x, p): return p[0] + p[1] * x

def func2(*args):
    return func(args[0],args[1:])

popt, pcov = curve_fit(func2, np.arange(10), np.arange(10), p0=(0, 0))
print popt,pcov

EDIT: This works for me

import numpy as np
from scipy.optimize import curve_fit

def func(x, *p): return p[0] + p[1] * x

popt, pcov = curve_fit(func, np.arange(10), np.arange(10), p0=(0, 0))
print popt,pcov