fit a parabola over 2D data -- edited

Question

I have searched for an example code 2-D data fitting using parabola and/or hyperbola, I was not successful to run any.

I have used scipy.optimize.leastsq for my fitting exercise:

NEW modified the code as follows:

def hyprsinc_errors(pararr, t,x,datarr):
    x=np.array(x)
    pararr[4] = np.abs(pararr[4])
    outarr = np.zeros((np.size(t),np.size(x)),float)
    for ix in x:
        it = int(np.sqrt(pararr[0]*(ix-pararr[1])**2)+pararr[2])
        if it < max(t)-int(pararr[4])-1:
            for iit in range(-int(pararr[4])+it,it+int(pararr[4])+1,1):
                outarr[iit,ix] =(-1)**int(pararr[4])*pararr[0]/pararr[3]
           #end for it
            outarr[it,ix]=pararr[3]
   #end for ix
    output =  ((outarr-datarr).ravel()).sum()
    print(type(output),output)
    return np.float(output)
#
import numpy as np
import matplotlib as plt
from scipy.optimize import curve_fit
#
datarr = np.array([[ 4,  0,  1,  0,  2,  3,  1,  5,  2,  0],\
                  [ 2,  0,  0,  2,  1,  0,  5,  5,  3,  5],\
                  [ 4,  2,  0,  2,  0,  1,  5,  4,  3,  4],\
                  [ 2,  0,  1,  3,  5,  2,  3,  5,  3,  3],\
                  [ 5,  3,  3,  4, 12, 12,  5,  0,  2,  3],\
                  [ 2,  0,  5, 12, 12, 11, 13,  0,  4,  3],\
                  [ 5,  3, 12, 11,  2,  2, 10, 15,  2,  3],\
                  [ 1, 15, 11,  3,  4,  0,  0, 11, 10,  3],\
                  [14, 12,  1,  1,  2,  5,  3,  2, 12, 14],\
                  [10,  3,  4,  4,  1,  4,  0,  5,  4, 10]])
#
T = np.linspace(0, 9, 10)
X = np.linspace(0, 9, 10)


hyprsinc_errors((T,X),datarr, 1,4,4,10,0)
optimized_result = leastsq(hyprsinc_errors,x0=np.array([1.,5,5,10,0]),args=(T,X,datarr))
print("opt_result = ", optimized_result[0])

I receive new error:

In [92]: p1,success = leastsq(hyprsinc_errors,x0=np.array([1.,5,5,10,0]),args=(T,X,datarr))

(, -388.0) <-- print in the last line of error function!

> - TypeError                                 Traceback (most recent call last)
> /nfs/rvl/sip/gs/nobackup3/holland/interferometry16/multiples/python_build/intel/python/<ipython-input-92-5f7acb18d23f>
> in <module>()
> ----> 1 p1,success = leastsq(hyprsinc_errors,x0=np.array([1.,5,5,10,0]),args=(T,X,datarr))
> 
> /apps/sss/epd/7.2.2/lib/python2.7/site-packages/scipy/optimize/minpack.py
> in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol,
> gtol, maxfev, epsfcn, factor, diag)
>     276     m = _check_func('leastsq', 'func', func, x0, args, n)[0]
>     277     if n > m:
> --> 278         raise TypeError('Improper input: N=%s must not exceed M=%s' % (n,m))
>     279     if Dfun is None:
>     280         if (maxfev == 0):
> 
> TypeError: Improper input: N=5 must not exceed M=1

---

Answer 1

curve_fit is more suited towards fitting 2-dimensional functions. The error you see is because you're passing a 2-dimensional array as the ydata parameter, which is meant to be an m-length array (i.e. one dimensional). Fitting the parameters along a single slice of the function will yield incorrect results most of the time.

I suggest you use least_squares instead of curve_fit. Bear in mind that it is a bit lower level: it requires you to manually calculate the errors and to provide a guess on the parameters:

from scipy.optimize import least_squares
XX,YY = np.meshgrid(X,Y)
prbola_errors = lambda args: (prbola(XX,YY, *args) - dataarr).ravel()
optimize_result = least_squares(prbola_errors, (0., 1., 0., 1.))

Here's a demonstration of this: https://gist.github.com/FranciscoDA/378b2223957d2b0e201350b0e66aec84

Edit: About the updated question, some things needed to be fixed:

• leastsq() expects the passed function to return a 1d-array with the calculated errors
• The number of arguments passed as the x0 argument of leastsq() didn't match the number of arguments expected by the function (here I added a dummy variable initialized to 0)
• The function tries to index an array using non-int scalars. Here I've casted them to int where needed

Updated code:

import numpy as np
import matplotlib as plt
from scipy.optimize import curve_fit, leastsq

def hyprsinc_errors(pararr, t,x,datarr):
    x=np.array(x)
    pararr[4] = np.abs(pararr[4])
    outarr = np.zeros((np.size(t),np.size(x)),float)
    for ix in x:
        it = int(np.sqrt(pararr[0]*(ix-pararr[1])**2)+pararr[2])
        if it < max(t)-int(pararr[4])-1:
            for iit in range(-int(pararr[4])+it,it+int(pararr[4])+1,1):
                outarr[int(iit),int(ix)] =(-1)**int(pararr[4])*pararr[0]/pararr[3]
            outarr[int(it),int(ix)]=pararr[3]
    output = (outarr-datarr).ravel()
    return output
#

#
datarr = np.array([[ 4,  0,  1,  0,  2,  3,  1,  5,  2,  0],\
                  [ 2,  0,  0,  2,  1,  0,  5,  5,  3,  5],\
                  [ 4,  2,  0,  2,  0,  1,  5,  4,  3,  4],\
                  [ 2,  0,  1,  3,  5,  2,  3,  5,  3,  3],\
                  [ 5,  3,  3,  4, 12, 12,  5,  0,  2,  3],\
                  [ 2,  0,  5, 12, 12, 11, 13,  0,  4,  3],\
                  [ 5,  3, 12, 11,  2,  2, 10, 15,  2,  3],\
                  [ 1, 15, 11,  3,  4,  0,  0, 11, 10,  3],\
                  [14, 12,  1,  1,  2,  5,  3,  2, 12, 14],\
                  [10,  3,  4,  4,  1,  4,  0,  5,  4, 10]])
#
T = np.linspace(0, 9, 10)
X = np.linspace(0, 9, 10)

optimized_result = leastsq(hyprsinc_errors,x0=np.array([1.,5,5,10,0,0]),args=(T,X,datarr))
print("opt_result = ", optimized_result[0])

Sample output: opt_result = [ 1. 5.00000006 5. 10.42857143 0. 0. ]