I tested the problem with: 1. Numeric 23.1 under python 2.3.2 2. numarray 0.8 (I made a copy of the Scientific package where all calls to Numeric were replaced to numarray), under python 2.3.2 There results where about the same -- high coefficients for the 5th order polynomials. I would expect reliable fit for a high order polynomials only under very special circumstances, so this is not a big surprise. My advice is: * Make sure that this is a bug and not a result of a numerical instability. If you can trace it down and point to a bug, then report it. The numarray package is very usable and is under a very active and rapid development, thus bugs are being fixed fast. * Look for a solution in the scipy package: It is generally better then Scientific. * Polynomials fit is relatively very simple --- you may write one of you own in less then a one day work. Since, as I said, the problem is, in many cases, unstable, you'll have the chance to implement more stable linear-equation solvers. Nadav. On Fri, 2003-10-31 at 15:19, Rob W.W. Hooft wrote:
I am using Polynomial.py from Scientific Python 2.1, together with Numeric 17.1.2. This has always served me well, but now we are busy upgrading our software, and I am currently porting some code to Scientific Python 2.4.1, Numeric 22.0. Suddenly I do no longer manage to get proper 2D polynomial fits over 4x4th order. At 5x5 the coefficients that come back from LinearAlgebra.linear_least_squares have exploded. In the old setup, I easily managed 9x9th order if I needed to, but most of the time I'd stop at 6x6th order. Would anyone have any idea how this difference can come about? I managed to work around this for the moment by using the equivalent code in the fitPolynomial routine that uses LinearAlgebra.generalized_inverse (and it doesn't even have problems with the same data at 8x8), but this definitely feels not right! I can't remember reading anything like this here before.
Together with Konrad Hinsen, I came to the conclusion that the problem is not in Scientific Python, so it must be the underlying LinearAlgebra code that changed between releases 17 and 22.
I hacked up a simplified example. Not sure whether it is the most simple case, but this resembles what I have in my code, and I'm quite sure it worked with Numeric 17.x, but currently it is horrible over order (4,4):
-------------------------------------- import Numeric
def func(x,y): return x+0.1*x**2+0.01*x**4+0.002*x**6+0.03*x*y+0.001*x**4*y**5
x=[] y=[] z=[] for dx in Numeric.arange(0,1,0.01): for dy in Numeric.arange(0,1,0.01): x.append(dx) y.append(dy) z.append(func(dx,dy))
from Scientific.Functions import Polynomial data=Numeric.transpose([x,y]) z=Numeric.array(z) for i in range(10): print data[i],z[i] pol=Polynomial.fitPolynomial((4,4),data,z) print pol.coeff ------------------------------------ for 4,4 this prints: [[ 1.84845529e-05 -7.60502772e-13 2.71314749e-12 -3.66731796e-12 1.66977148e-12] [ 9.99422967e-01 3.00000000e-02 -3.26346097e-11 4.42406519e-11 -2.01549767e-11] [ 1.03899464e-01 -3.19668064e-11 1.14721790e-10 -1.55489826e-10 7.08425891e-11] [ -9.40275000e-03 4.28456838e-11 -1.53705205e-10 2.08279772e-10 -9.48840470e-11] [ 1.80352695e-02 -1.10999843e-04 8.00662570e-04 -2.17266676e-03 2.47500004e-03]]
for 5,5:
[[ -2.25705839e+03 6.69051337e+02 -6.60470163e+03 6.66572425e+03 -8.67897022e+02 1.83974866e+03] [ -2.58646837e+02 -2.46554689e+03 1.15965805e+03 7.01089888e+03 -2.11395436e+03 2.10884815e+03] [ 3.93307499e+03 4.34484805e+02 -4.84080392e+03 5.90375330e+03 1.16798049e+03 -4.14163933e+03] [ 1.62814750e+03 2.08717457e+03 1.15870693e+03 -3.37838057e+03 3.49821689e+03 5.80572585e+03] [ 4.54127557e+02 -1.56645524e+03 4.58997025e+00 1.69772635e+03 -1.37751039e+03 -7.59726558e+02] [ 2.37878239e+03 9.43032094e+02 8.58518644e+02 -8.35846339e+03 -5.55845668e+02 1.87502761e+03]]
Which is clearly wrong.
I appreciate any help!
Regards,
Rob