181
182
def polyfit(x, y, deg, rcond=None, full=False):
182
183
"""Least squares polynomial fit.
186
x -- vector of sample points
187
y -- vector or 2D array of values to fit
188
deg -- degree of the fitting polynomial
192
rcond -- relative condition number of the fit (default len(x)*eps)
193
full -- return full diagnostic output (default False)
197
full == False -- coefficients
198
full == True -- coefficients, residuals, rank, singular values, rcond.
202
RankWarning -- if rank is reduced and not full output
204
185
Do a best fit polynomial of degree 'deg' of 'x' to 'y'. Return value is a
205
186
vector of polynomial coefficients [pk ... p1 p0]. Eg, for n=2
207
p2*x0^2 + p1*x0 + p0 = y1
208
p2*x1^2 + p1*x1 + p0 = y1
209
p2*x2^2 + p1*x2 + p0 = y2
211
p2*xk^2 + p1*xk + p0 = yk
214
Method: if X is a the Vandermonde Matrix computed from x (see
188
p2*x0^2 + p1*x0 + p0 = y1
189
p2*x1^2 + p1*x1 + p0 = y1
190
p2*x2^2 + p1*x2 + p0 = y2
192
p2*xk^2 + p1*xk + p0 = yk
197
1D vector of sample points.
199
1D vector or 2D array of values to fit. The values should run down the
200
columes in the 2D case.
202
Degree of the fitting polynomial
203
rcond: {None, float}, optional
204
Relative condition number of the fit. Singular values smaller than this
205
relative to the largest singular value will be ignored. The defaul value
206
is len(x)*eps, where eps is the relative precision of the float type,
207
about 2e-16 in most cases.
208
full : {False, boolean}, optional
209
Switch determining nature of return value. When it is False just the
210
coefficients are returned, when True diagnostic information from the
211
singular value decomposition is also returned.
215
coefficients, [residuals, rank, singular_values, rcond] : variable
216
When full=False, only the coefficients are returned, running down the
217
appropriate colume when y is a 2D array. When full=True, the rank of the
218
scaled Vandermonde matrix, it's effective rank in light of the rcond
219
value, its singular values, and the specified value of rcond are also
224
RankWarning : if rank is reduced and not full output
225
The warnings can be turned off by:
226
>>> import numpy as np
228
>>> warnings.simplefilter('ignore',np.RankWarning)
233
polyval : computes polynomial values.
237
If X is a the Vandermonde Matrix computed from x (see
215
238
http://mathworld.wolfram.com/VandermondeMatrix.html), then the
216
239
polynomial least squares solution is given by the 'p' in
220
where X is a len(x) x N+1 matrix, p is a N+1 length vector, and y
221
is a len(x) x 1 vector
243
where X.shape is a matrix of dimensions (len(x), deg + 1), p is a vector of
244
dimensions (deg + 1, 1), and y is a vector of dimensions (len(x), 1).
223
246
This equation can be solved as
225
p = (XT*X)^-1 * XT * y
248
p = (XT*X)^-1 * XT * y
227
250
where XT is the transpose of X and -1 denotes the inverse. However, this
228
251
method is susceptible to rounding errors and generally the singular value
229
decomposition is preferred and that is the method used here. The singular
230
value method takes a paramenter, 'rcond', which sets a limit on the
231
relative size of the smallest singular value to be used in solving the
232
equation. This may result in lowering the rank of the Vandermonde matrix,
233
in which case a RankWarning is issued. If polyfit issues a RankWarning, try
234
a fit of lower degree or replace x by x - x.mean(), both of which will
252
decomposition of the matrix X is preferred and that is what is done here.
253
The singular value method takes a paramenter, 'rcond', which sets a limit on
254
the relative size of the smallest singular value to be used in solving the
255
equation. This may result in lowering the rank of the Vandermonde matrix, in
256
which case a RankWarning is issued. If polyfit issues a RankWarning, try a
257
fit of lower degree or replace x by x - x.mean(), both of which will
235
258
generally improve the condition number. The routine already normalizes the
236
259
vector x by its maximum absolute value to help in this regard. The rcond
237
parameter may also be set to a value smaller than its default, but this may
238
result in bad fits. The current default value of rcond is len(x)*eps, where
260
parameter can be set to a value smaller than its default, but the resulting
261
fit may be spurious. The current default value of rcond is len(x)*eps, where
239
262
eps is the relative precision of the floating type being used, generally
240
263
around 1e-7 and 2e-16 for IEEE single and double precision respectively.
241
264
This value of rcond is fairly conservative but works pretty well when x -
242
265
x.mean() is used in place of x.
244
The warnings can be turned off by:
248
>>> warnings.simplefilter('ignore',numpy.RankWarning)
250
268
DISCLAIMER: Power series fits are full of pitfalls for the unwary once the
251
269
degree of the fit becomes large or the interval of sample points is badly
252
centered. The basic problem is that the powers x**n are generally a poor
253
basis for the functions on the sample interval with the result that the
254
Vandermonde matrix is ill conditioned and computation of the polynomial
255
values is sensitive to coefficient error. The quality of the resulting fit
256
should be checked against the data whenever the condition number is large,
257
as the quality of polynomial fits *can not* be taken for granted. If all
258
you want to do is draw a smooth curve through the y values and polyfit is
259
not doing the job, try centering the sample range or look into
260
scipy.interpolate, which includes some nice spline fitting functions that
270
centered. The problem is that the powers x**n are generally a poor basis for
271
the polynomial functions on the sample interval, resulting in a Vandermonde
272
matrix is ill conditioned and coefficients sensitive to rounding erros. The
273
computation of the polynomial values will also sensitive to rounding errors.
274
Consequently, the quality of the polynomial fit should be checked against
275
the data whenever the condition number is large. The quality of polynomial
276
fits *can not* be taken for granted. If all you want to do is draw a smooth
277
curve through the y values and polyfit is not doing the job, try centering
278
the sample range or look into scipy.interpolate, which includes some nice
279
spline fitting functions that may be of use.
263
281
For more info, see
264
282
http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html,
265
283
but note that the k's and n's in the superscripts and subscripts
266
284
on that page. The linear algebra is correct, however.
271
287
order = int(deg) + 1
272
288
x = NX.asarray(x) + 0.0
318
339
def polyval(p, x):
319
"""Evaluate the polynomial p at x. If x is a polynomial then composition.
323
If p is of length N, this function returns the value:
324
p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
326
x can be a sequence and p(x) will be returned for all elements of x.
327
or x can be another polynomial and the composite polynomial p(x) will be
330
Notice: This can produce inaccurate results for polynomials with
331
significant variability. Use carefully.
340
"""Evaluate the polynomial p at x.
342
If p is of length N, this function returns the value:
344
p[0]*(x**N-1) + p[1]*(x**N-2) + ... + p[N-2]*x + p[N-1]
346
If x is a sequence then p(x) will be returned for all elements of x. If x is
347
another polynomial then the composite polynomial p(x) will be returned.
351
p : {array_like, poly1d}
352
1D array of polynomial coefficients from highest degree to zero or an
354
x : {array_like, poly1d}
355
A number, a 1D array of numbers, or an instance of poly1d.
359
values : {array, poly1d}
360
If either p or x is an instance of poly1d, then an instance of poly1d is
361
returned, otherwise a 1D array is returned. In the case where x is a
362
poly1d, the result is the composition of the two polynomials, i.e.,
363
substitution is used.
367
Horners method is used to evaluate the polynomial. Even so, for polynomial
368
if high degree the values may be inaccurate due to rounding errors. Use
333
372
p = NX.asarray(p)
334
373
if isinstance(x, poly1d):
added
removed
numpy/lib/polynomial.py
