Python curve fitting for a function with input vectors(변수가 여러개인 함수의 fitting)

변수가 여러개인 함수 예를 들어 f(d,t; c_0,c_1...) 인 경우에는
어떻게 fitting을 할까? 여기서, d,t는 변수이고, c_0,c_1..등은
fitting parameter 이다.

Suppose we have data which corresponds to

y=f[d,t]

Something like (d,t)=(0,0) gives y=0, (0,1) gives y=1....

We want to fit these data with certain fitting function
with parameters c1,c2,c3...

fit[d,t; c1,c2,c3]

Even in this case, we can use curve_fit
(Following example copied from http://stackoverflow.com/a/27096056/2775514

from scipy.optimize import curve_fit
import scipy

def fn(x, a, b, c):
    return a + b*x[0] + c*x[1]

# y(x0,x1) data:
#    x0=0 1 2
# ___________
# x1=0 |0 1 2
# x1=1 |1 2 3
# x1=2 |2 3 4

x = scipy.array([[0,1,2,0,1,2,0,1,2,],[0,0,0,1,1,1,2,2,2]])
y = scipy.array([0,1,2,1,2,3,2,3,4])
popt, pcov = curve_fit(fn, x, y)
print popt

In this example,
x[:, i-th], y[i-th] corresponds to  i-th data set.

result returns popt -> [0,1,1].

One can use different form : 

x = np.array([ [1.0, 1.0 ],[1.0,1.5],[2.0,0.5],[2.1,1.6],[0.5,0.5],[0.8,1.2] ]  )

y = 0.5*x[:,0] + 2.0 *x[:,1]**2

def fitf(x,*paras):

    return paras[0]*x[:,0]+paras[1]*x[:,1]**2

xx = np.array( [x[:,0],x[:,1]])

def testf(x,*para):

    return para[0]*x[0]+x[1]**2*para[1]

print( curve_fit(fitf,x,y,p0=[1,1]) )

print( curve_fit(testf,xx,y,p0=[1,1]) )

Note that curve_fit calls the function as 

ydata = f(xdata, *params) +eps 

In other words, when xdata is an (k,M) array with M-data, 

the function should return array ydata with M-data. 

if function only works work for a scalar f(x), it cannot be used directly with curve_fit. 

On the other hand, leastsq(cost, p0) assumes the 

cost function  returns a number or array of length of p0

#--------------------------------------------------------------------------------------------

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, **kw)[source]

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

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

Parameters:	f : callable The model function, f(x, ...). It must take the independent variable as the first argument and the parameters to fit as separate remaining arguments. xdata : An M-length sequence or an (k,M)-shaped array for functions with k predictors. The independent variable where the data is measured. ydata : M-length sequence The dependent data — nominally f(xdata, ...) p0 : None, scalar, or N-length sequence Initial guess for the parameters. If None, then the initial values will all be 1 (if the number of parameters for the function can be determined using introspection, otherwise a ValueError is raised). sigma : None or M-length sequence, optional If not None, the uncertainties in the ydata array. These are used as weights in the least-squares problem i.e. minimising np.sum( ((f(xdata, *popt) -ydata) / sigma)**2 ) If None, the uncertainties are assumed to be 1. absolute_sigma : bool, optional If False, sigma denotes relative weights of the data points. The returned covariance matrix pcov is based on estimated errors in the data, and is not affected by the overall magnitude of the values in sigma. Only the relative magnitudes of the sigma values matter. If True, sigma describes one standard deviation errors of the input data points. The estimated covariance in pcov is based on these values. check_finite : bool, optional If True, check that the input arrays do not contain nans of infs, and raise a ValueError if they do. Setting this parameter to False may silently produce nonsensical results if the input arrays do contain nans. Default is True.
Returns:	popt : array Optimal values for the parameters so that the sum of the squared error of f(xdata, *popt) - ydata is minimized pcov : 2d array The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov)). How the sigma parameter affects the estimated covariance depends on absolute_sigma argument, as described above.
Raises:	OptimizeWarning if covariance of the parameters can not be estimated. ValueError if ydata and xdata contain NaNs.

python fitting / optimization

주어진 xdata와 ydata에 대해,
Least square method로 data를 fitting 하는 방법.

결론은 , scipy의 optimize module을 이용한다.

import scipy.optimize as optimization

def func(x,para1,para2):
   return para1*sin(para2*x)

popt, pcov = optimization.curve_fit(func,xdata, ydata, p0=(1.0,1.0))

p0는 initial guess이고

popt는 parameter fitting 의 결과를, pcov는 covariance 계산 결과를 보여준다.

least square method를 사용하므로,
cost function chi(para...) 함수를 최소로 만드는 parameter를 찾는다.

pcov는 (parameter 갯수)*(parameter 갯수) matrix로,
cost function의 parameter 값의 변화에 대한 변화값을 나타낸다.
pcov(i,j)= d^2[chi(c1,c2..)]/dc_i dc_j

몇가지 다른 fitting 방법이 있다.

1. scipy.optimize.curve_fit : 

   입력: f(x,para) , xdata, ydata, yerr(optional)

    popt, pcov = optimization.curve_fit(f,xdata, ydata, p0=(1.0,1.0))

   curve_fit은 least square 방법을 통해  sum ( ydata - f(xdata,para))^2/yerr^2  

   을 minimize 시킨다. 

   parameter의 error estimation은 sqrt(diag(pcov)) 를 이용한다.

2. scipy.optimize.leastsq :

   입력: errFunc(para,x,y) , xdata, ydata 

    popt, hess_inv, infodict, errmsg, success = optimize.leastsq(errFunc, pstart, args=(xdata, ydata), full_output=1)   

    입력하는 함수는 errFunc(para,xdata,ydata) = ydata - f(xdata,para) 

    이고,  sum ( errFunc )^2 를 minimize 시킨다. 

    parameter의 error estimation은 hess_inv 와 residual variance 를 이용한다. 

    sqrt( diag( hess_inv * resVariance)) with 

    resVariance = sum( errFunc(para,xdata,ydata)^2) /(len(ydata)-len(para)) 

3. scipy.optimize.minimize :

    입력: minFunc(para,x,y) , xdata, ydata 

    

    result = optimize.minimize(lambda p,x,y: np.sum( errFunc(p,x,y)**2 ), pstart, args=(xdata, ydata))

    이 경우 입력함수 

    minFunc(para,xdata,ydata) = sum ( ydata- f(xdata,para))^2 이고, 

    minFunc을 minimize 시킨다.  parameter의 error estimation은 leastsq 와 같이

    result.hess_inv 와 residual variance 을 이용한다. 

    그러나, minimization의 방법에 따라 hess_inv가 제공되지 않는 경우가 있다.