Curve Fitting In Scipy With 3d Data And Parameters

I am working on fitting a 3d distribution function in scipy. I have a numpy array with counts in x- and y-bins, and I am trying to fit that to a rather complicated 3-d distributio

Solution 1:

So what leastsq does is try to:

"Minimize the sum of squares of a set of equations" -scipy docs

as it says it's minimizing a set of functions and therefore doesn't actually take any x or y data inputs in the easiest manner if you look at the arguments here so you can do it as you like and pass a residual function however, it's significantly easier to just use curve_fit which does it for you :) and creates the necessary equations

For fitting you should use: curve_fit if you are ok with the generic residual they use which is actually the function you pass itself res = leastsq(func, p0, args=args, full_output=1, **kw) if you look in the code here.

e.g. If I fit the rosenbrock function in 2d and guess the y-parameter:

from scipy.optimize import curve_fit
from itertools import imap
import numpy as np
# use only an even number of arguments
def rosen2d(x,a):
    return (1-x)**2 + 100*(a - (x**2))**2
#generate some random data slightly off

datax = np.array([.01*x for x in range(-10,10)])
datay = 2.3
dataz = np.array(map(lambda x: rosen2d(x,datay), datax))
optimalparams, covmatrix = curve_fit(rosen2d, datax, dataz)
print 'opt:',optimalparams

fitting the colville function in 4d:

from scipy.optimize import curve_fit
import numpy as np

# 4 dimensional colville function
# definition from http://www.sfu.ca/~ssurjano/colville.html
def colville(x,x3,x4):
    x1,x2 = x[:,0],x[:,1]
    return 100*(x1**2 - x2)**2 + (x1-1)**2 + (x3-1)**2 + \
            90*(x3**2 - x4)**2 + \
            10.1*((x2 - 1)**2 + (x4 - 1)**2) + \
            19.8*(x2 - 1)*(x4 - 1)
#generate some random data slightly off

datax = np.array([[x,x] for x in range(-10,10)])
#add gaussian noise
datax+= np.random.rand(*datax.shape)
#set 2 of the 4 parameters to constants
x3 = 3.5
x4 = 4.5
#calculate the function
dataz = colville(datax, x3, x4)
#fit the function
optimalparams, covmatrix = curve_fit(colville, datax, dataz)
print 'opt:',optimalparams

Using a custom residual function:

from scipy.optimize import leastsq
import numpy as np

# 4 dimensional colville function
# definition from http://www.sfu.ca/~ssurjano/colville.html
def colville(x,x3,x4):
    x1,x2 = x[:,0],x[:,1]
    return 100*(x1**2 - x2)**2 + (x1-1)**2 + (x3-1)**2 + \
            90*(x3**2 - x4)**2 + \
            10.1*((x2 - 1)**2 + (x4 - 1)**2) + \
            19.8*(x2 - 1)*(x4 - 1)
#generate some random data slightly off


datax = np.array([[x,x] for x in range(-10,10)])
#add gaussian noise
datax+= np.random.rand(*datax.shape)
#set 2 of the 4 parameters to constants
x3 = 3.5
x4 = 4.5

def residual(p, x, y):
    return y - colville(x,*p)
#calculate the function
dataz = colville(datax, x3, x4)
#guess some initial parameter values
p0 = [0,0]
#calculate a minimization of the residual
optimalparams = leastsq(residual, p0, args=(datax, dataz))[0]
print 'opt:',optimalparams

Edit: you used both the position and the keyword arg for args: if you look at the docs you'll see it uses position 3, but also can be used as a keyword argument. You used both which means the function is as expected, confused.