17.2.1  Least-squares polynomial approximation

The fitpoly command is used for replacing tabular data or a function by a polynomial.

Smoothing data with polynomial regression.

fitpoly can be used for finding a polynomial that fits a given tabular data.

• fitpoly takes one mandatory argument and two optional arguments:
  • data, a two-column matrix or a sequence of two lists, representing the x- and y-values of data points (x₀,y₀),(x₁,y₁),…,(x_m,y_m).
  • Optionally, x[=a..b] or a..b, which is a variable x and optionally its domain [a,b] or a segment [a,b] (by default, unset).
  • Optionally, a sequence of options, each of which may be one of:
    • limit=N, a nonnegative integer, which is the maximal degree for the resulting polynomial (by default, N=15).
    • degree=d, where d is a nonnegative integer, which is the desired degree of the resulting polynomial (by default, unset).
    • threshold=tol, a positive real number, which is used for finding an appropriate fitting degree, as described below (by default, threshold=0.01, ignored when d is set).
    • length=l, a positive integer, which is used for finding an appropriate fitting degree, as described below (by default, length=5, ignored when d is set).
• If d is set, then fitpoly returns the polynomial of degree min{d,N} which best fits the data in the least-squares sense.
• If d is not set, then a polynomial p_n of modest degree n≤ N but a good error suppression (thus representing the trend of the data) is chosen using tol and l such that raising the degree does not make a significant improvement to the approximation. Precisely, n is the first nonnegative integer such that

  stddev[σn2,σn+12,…,σn+l−12] mean[σ02,σ12,…,σn+l−12]

  ≤tol,

  where σ_k² is the sum of squared residuals ∑_i=0^m(y_i−p_k(x_i))² and p_k is the best-fitting polynomial of degree k. Ideally, n is the smallest degree for which σ_n²≈σ_n+1²≈σ_n+2²≈⋯, meaning that higher-order polynomials would overfit on the noise.
• Lowering the parameter length, as well as increasing threshold, would make the algorithm more greedy, effectively lowering the degree of the output polynomial. Note that when degree is set, these two parameters are ignored.
• If no variable is specified, then fitpoly returns the list of polynomial coefficients.
• If a segment [a,b] is given, then only the data points (x_k,y_k) for which a≤ x_k≤ b are considered.

Example

We find a polynomial which best approximates data infected by noise. Data is obtained by evaluating the polynomial f(x)=1/100x⁴−1/5x³+6x at 100 points between x=1 and x=20.8. The noise is generated by a random, normally-distributed variable.

We obtain the list of polynomial coefficients, starting with the leading coefficient. Note that the polynomial of degree 4 is returned (it has five coefficients), which is, in this case, obviously optimal for data smoothing. To visualize the fit, enter:

When approximating only the data for 1≤ x≤ 10, you obtain a polynomial of 9th degree. This is called overfitting. Namely, data curvature at the restriction is small, which makes the noise more prominent.

To make the approximating polynomial less sensitive to noise, you can increase the threshold value tol.

(Alternatively, you could set the parameter length to a smaller value, e.g. 3.) The command lines below show the curves p1 and p2 against data points (the overfitted one is drawn in grey).

Approximating functions by polynomials.

fitpoly can be used for approximating a continuous function f:[a,b]→ℝ by a polynomial p_n of certain degree n which minimizes the error in the sense of ℓ² norm.

• fitpoly takes two mandatory arguments and one optional argument:
  • f, an univariate symbolic expression representing the function to be approximated.
  • An interval a..b, which is the domain [a,b] of f.
  • Optionally, opts, a sequence of options each of which may be one of:
    • limit=N, as before.
    • degree=d, as before.
    • ε or threshold=ε, which is a positive real number such that for the resulting polynomial p_n the following holds:

      ||f−pn||2=

      ∫

      b a

      (f(x)−pn(x))2dx≤ε     (1)

      (by default, ε=10⁻⁸). This parameter is ignored when d is set.
• fitpoly(f,a..b ⟨,opts ⟩) finds p_n with smallest n≤ N such that (1) holds. If such n does not exist, then p_N is returned.
• If the parameter d is set, then p_d which minimizes ||f−p_d|| among all polynomials of degree d is returned.
• Polynomial approximation is fast and robust even for large d and N resp. for small ε.

Examples

The mean absolute error of the above approximation can be estimated as follows: