[Numpy-svn] r3338 - trunk/numpy/lib
numpy-svn at scipy.org
numpy-svn at scipy.org
Sat Oct 14 16:34:17 CDT 2006
Author: charris
Date: 2006-10-14 16:34:16 -0500 (Sat, 14 Oct 2006)
New Revision: 3338
Modified:
trunk/numpy/lib/polynomial.py
Log:
Scale the x vector in polyfit to improve condition of Vandermonde matrix.
Throw error on rank reduction in polyfit.
Error should probably be replace with a warning at some point.
Add argument checks in polyfit.
Modified: trunk/numpy/lib/polynomial.py
===================================================================
--- trunk/numpy/lib/polynomial.py 2006-10-14 21:22:07 UTC (rev 3337)
+++ trunk/numpy/lib/polynomial.py 2006-10-14 21:34:16 UTC (rev 3338)
@@ -7,9 +7,10 @@
'polyfit']
import re
+import warnings
import numpy.core.numeric as NX
-from numpy.core import isscalar
+from numpy.core import isscalar, abs
from numpy.lib.twodim_base import diag, vander
from numpy.lib.shape_base import hstack, atleast_1d
from numpy.lib.function_base import trim_zeros, sort_complex
@@ -169,11 +170,11 @@
val = poly1d(val)
return val
-def polyfit(x, y, N, rcond=-1):
+def polyfit(x, y, deg, rcond=-1):
"""
- Do a best fit polynomial of degree N of y to x. Return value is a
- vector of polynomial coefficients [pk ... p1 p0]. Eg, for N=2
+ Do a best fit polynomial of degree 'deg' of 'x' to 'y'. Return value is a
+ vector of polynomial coefficients [pk ... p1 p0]. Eg, for n=2
p2*x0^2 + p1*x0 + p0 = y1
p2*x1^2 + p1*x1 + p0 = y1
@@ -201,21 +202,22 @@
value method takes a paramenter, 'rcond', which sets a limit on the
relative size of the smallest singular value to be used in solving the
equation. The default value of rcond is the double precision machine
- precision as the actual solution is carried out in double precision. If you
- are simply interested in a polynomial line drawn through the data points
- and *not* in a true power series expansion about zero, then the best bet is
- to scale the x sample points to the interval [0,1] as the problem will
- probably be much better posed.
+ precision as the actual solution is carried out in double precision. The
+ vector x is also normalized by its maximum absolute value before the fit
+ to improve the condition of the Vandermonde matrix.
WARNING: Power series fits are full of pitfalls for the unwary once the
- degree of the fit get up around 4 or 5 and the interval of sample points
- gets large. Computation of the polynomial values are sensitive to
- coefficient errors and the Vandermonde matrix is ill conditioned. The knobs
- available to tune the fit are degree and rcond. The rcond knob is a bit
- flaky and it can be useful to use values of rcond less than the machine
- precision, 1e-24 for instance, but the quality of the resulting fit needs
- to be checked against the data. The quality of polynomial fits *can not* be
- taken for granted.
+ degree of the fit becomes large or the interval of sample points is badly
+ uncentered. The basic problem is that the powers x**n are not generally a
+ good basis for the functions on the sample interval with the result that
+ the Vandermonde matrix is ill conditioned and computation of the polynomial
+ values is sensitive to coefficient error. The quality of the resulting fit
+ should be checked against the data whenever the condition number is large,
+ as the quality of polynomial fits *can not* be taken for granted. If all
+ you want to do is draw a smooth curve through the y values and polyfit is
+ not doing the job, try centering the sample range or look into
+ scipy.interpolate, which includes some nice spline fitting functions that
+ may be of use.
For more info, see
http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html,
@@ -225,10 +227,38 @@
See also polyval
"""
+ order = deg + 1.
x = NX.asarray(x) + 0.0
y = NX.asarray(y) + 0.0
- v = vander(x, N+1)
+
+ # check arguments.
+ if deg < 0 :
+ raise ValueError, "expected deg >= 0"
+ if x.ndim != 1 or x.size == 0:
+ raise TypeError, "expected non-empty vector for x"
+ if y.ndim < 1 or y.ndim > 2 :
+ raise TypeError, "expected 1D or 2D array for y"
+ if x.shape[0] != y.shape[0] :
+ raise TypeError, "expected x and y to have same length"
+
+ # scale x to improve condition number
+ scale = abs(x).max()
+ if scale != 0 :
+ x /= scale
+
+ # solve least squares equation for powers of x
+ v = vander(x, order)
c, resids, rank, s = _lstsq(v, y, rcond)
+
+ # warn on rank reduction, which indicates an ill conditioned matrix
+ if rank != order :
+ # fixme -- want a warning here, not an exception
+ raise RuntimeError, "degree reduced to %d" % (rank - 1)
+
+ # scale returned coefficients
+ if scale != 0 :
+ c /= vander([scale], order)[0]
+
return c
More information about the Numpy-svn
mailing list