Python minimization leastsq with columns

Question

Please excuse me but I started recently using Python...and code. I have some troubles with optimization in python... I would like to determine coefficients k0,k1...of a polynom defined in residual_x with data of vector "T" for each column of 2D array "test" (as example here with 3 columns), is it possible? and store results as 2D array.

As you can see, I do for only a list of data.

Thanks a lot for help!!!

import numpy as np
import scipy
import pylab 
from scipy.optimize import leastsq

T = np.arange(5,56,5)
T = np.reshape(T,(11,1))

test = np.array([[  3051.11,   2984.85,   3059.17],
       [  3510.78,   3442.43,   3520.7 ],
       [  4045.91,   3975.03,   4058.15],
       [  4646.37,   4575.01,   4662.29],
       [  5322.75,   5249.33,   5342.1 ],
       [  6102.73,   6025.72,   6127.86],
       [  6985.96,   6906.81,   7018.22],
       [  7979.81,   7901.04,   8021.  ],
       [  9107.18,   9021.98,   9156.44],
       [ 10364.26,  10277.02,  10423.1 ],
       [ 11776.65,  11682.76,  11843.18]])

k0 = 100.
k1 = 100.
k2 = 1.
k3 = 1.
k4 = 1.
k5 = 10.

def residual_x(vars, T, donnees):
    k0 = vars[0]
    k1 = vars[1]
    k2 = vars[2]
    k3 = vars[3]
    k4 = vars[4]
    k5 = vars[5]
    modele = T**5*k0+T**4*k1+T**3*k2+T**2*k3+T*k4+k5
    return (donnees-modele)

vars = [k0, k1, k2, k3, k4, k5]
out_x = leastsq(residual_x, vars, args=(T, test),epsfcn=0.001)

Answer 1

I think you're trying to fit three different 5th-degree polynomials with a common set of domain points, {5,10,15,...,55}, to three different range sets (the columns of test).

Your original code had indentation errors: from k0 = 100. and down was incorrectly indented; that was probably a cut-and-paste error.

More importantly, you need to fit each polynomial separately, so you should call scipy.optimize.leastsq for each column. This is easiest to do with a for loop (see code). You also need to collect the coefficients for each polynomial; you could do this in a 3x6 array; I've used a list in the example below.

EDIT: I was focused on getting your code working when I answered this, so I completely forgot that fitting polynomial coefficients doesn't need a non-linear solver. You should take a look at numpy.polyfit, which is a convenience wrapper around a least-squares linear solver for polynomial coefficients. You can find out more in the numpy docs or at Wikipedia, and I've added sample usage of numpy.polyfit to the example below.

import numpy as np
from scipy.optimize import leastsq

T = np.arange(5,56,5)

test = np.array([[  3051.11,   2984.85,   3059.17],
                 [  3510.78,   3442.43,   3520.7 ],
                 [  4045.91,   3975.03,   4058.15],
                 [  4646.37,   4575.01,   4662.29],
                 [  5322.75,   5249.33,   5342.1 ],
                 [  6102.73,   6025.72,   6127.86],
                 [  6985.96,   6906.81,   7018.22],
                 [  7979.81,   7901.04,   8021.  ],
                 [  9107.18,   9021.98,   9156.44],
                 [ 10364.26,  10277.02,  10423.1 ],
                 [ 11776.65,  11682.76,  11843.18]])


k0 = 100.
k1 = 100.
k2 = 1.
k3 = 1.
k4 = 1.
k5 = 10.

def residual_x(vars, T, donnees):
    k0 = vars[0]
    k1 = vars[1]
    k2 = vars[2]
    k3 = vars[3]
    k4 = vars[4]
    k5 = vars[5]
    modele = T**5*k0+T**4*k1+T**3*k2+T**2*k3+T*k4+k5
    return donnees-modele

vars = [k0, k1, k2, k3, k4, k5]
coeffs=[]

for i in range(test.shape[1]):
#EDIT: previously appended the complete output of leastsq to coeffs list
#      this is wrong; leastsq returns two outputs (i.e., a tuple)
#      and we're only really interested in the coefficients
    thiscoeffs,_=leastsq(residual_x, vars, args=(T, test[:,i]),epsfcn=0.001)
    coeffs.append(thiscoeffs)

for i in range(test.shape[1]):
    print 'poly[%d] coeffs: %s' % (i,str(coeffs[i]))

# try first example with numpy.polyfit
coeffs2 = np.polyfit(T,test[:,0],5)
print 'np.polyfit for case 0: %s',str(coeffs2)

#EDIT: this added in second edit
# if you want coeffs as an array, you could convert the list like this
coeffs_array=np.asarray(coeffs)
# or you could make coeffs a numpy array in the first place

This gives [EDIT: coeffs2 solution added]

poly[0] coeffs: (array([ -9.70769231e-07,   1.56254079e-04,   5.27449301e-03,
         1.03258957e+00,   7.62019210e+01,   2.64248394e+03]), 3)
poly[1] coeffs: (array([ -1.37025641e-06,   2.10074592e-04,   2.49477855e-03,
         1.09748240e+00,   7.51708855e+01,   2.58007182e+03]), 2)
poly[2] coeffs: (array([ -1.09435897e-06,   1.58983683e-04,   5.97712121e-03,
         1.01704522e+00,   7.67132035e+01,   2.64824061e+03]), 3)
np.polyfit for case 0: %s [ -9.70769231e-07   1.56254079e-04   5.27449301e-03   1.03258957e+00
   7.62019210e+01   2.64248394e+03]